Commit 0c7fa12a authored by Médéric Boquien's avatar Médéric Boquien
Browse files

Rather than reinterpolating all the components when adding a new one, let's be...

Rather than reinterpolating all the components when adding a new one, let's be smart about that and compute the interpolation only for new wavelengths. At the same time, make use of memoisation not to repeat the same computation all the time when marging different wavelength grids.
parent 11534c3b
......@@ -35,7 +35,6 @@ from collections import OrderedDict
from . import utils
from .io.vo import save_sed_to_vo
from scipy.constants import c, parsec
from scipy.interpolate import interp1d
from ..data import Database
......@@ -196,26 +195,14 @@ class SED(object):
# grid, we interpolate everything on a common wavelength grid.
if (results_wavelengths.size != self.wavelength_grid.size or
not np.all(results_wavelengths == self.wavelength_grid)):
# Compute the new wavelength grid for the spectrum.
new_wavelength_grid = utils.best_grid(results_wavelengths,
self.wavelength_grid)
# Interpolate each luminosity component to the new wavelength
# grid setting everything outside the wavelength domain to 0.
new_luminosities = interp1d(self.wavelength_grid,
self.wavelength_grid, self.luminosities = \
utils.interpolate_lumin(self.wavelength_grid,
self.luminosities,
bounds_error=False,
assume_sorted=True,
fill_value=0.)(new_wavelength_grid)
# Interpolate the added luminosity array to the new wavelength
# grid
interp_lumin = np.interp(new_wavelength_grid,
results_wavelengths, results_lumin,
left=0., right=0.)
self.wavelength_grid = new_wavelength_grid
self.luminosities = np.vstack((new_luminosities, interp_lumin))
results_wavelengths,
results_lumin)
self.luminosity = self.luminosities.sum(0)
else:
self.luminosities = np.vstack((self.luminosities,
......
......@@ -5,6 +5,7 @@
import numpy as np
from scipy.constants import c, pi
from scipy.interpolate import interp1d
def lambda_to_nu(wavelength):
......@@ -41,67 +42,6 @@ def nu_to_lambda(frequency):
return 1.e-9 * c / frequency
def best_grid_memoise(f):
"""
Memoisation helper for the best_grid() function. Given that best_grid takes
numpy arrays that are not hashable as input, we cannot use the standard
memoisation functions. Here we use the array sizes, minimum and maximum
values as hashes.
Parameters
----------
f: function
Function to memoise.
Returns
-------
best_grid_helper: function
Meomoised best_grid function
"""
memo = {}
def best_grid_helper(x, y):
sx = x.size
minx = x[0]
maxx = x[-1]
sy = y.size
miny = y[0]
maxy = y[-1]
if (sx, minx, maxx, sy, miny, maxy) not in memo:
memo[(sx, minx, maxx, sy, miny, maxy)] = f(x, y)
return memo[(sx, minx, maxx, sy, miny, maxy)]
return best_grid_helper
@best_grid_memoise
def best_grid(wavelengths1, wavelengths2):
"""
Return the best wavelength grid to regrid to arrays
Considering the two wavelength grids passed in parameters, this function
compute the best new grid that will be used to regrid the two spectra
before combining them. We do not use np.unique as it is much slowe than
finding the unique elements by hand.
Parameters
----------
wavelengths1, wavelengths2: array of floats
The wavelength grids to be 'regridded'.
Returns
-------
new_grid: array of floats
Array containing all the wavelengths found in the input arrays.
"""
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:])
return wl[flag]
def luminosity_to_flux(luminosity, dist):
"""
Convert a luminosity (or luminosity density) to a flux (or flux density).
......@@ -228,3 +168,201 @@ def redshift_spectrum(wavelength, flux, redshift, is_fnu=False):
flux = lambda_flambda_to_fnu(wavelength, flux)
return wavelength, flux
def memoise_1var(f):
"""
Memoisation helper for a function taking 1 numpy array as input. As it is
not hashable, we cannot use the standard memoisation function. Here we
we use the array size, minimum, and maximum values as hashes.
Parameters
----------
f: function
Function to memoise.
Returns
-------
memoise_helper: function
Meomoised best_grid function
"""
memo = {}
def memoise_helper(x):
sx = x.size
minx = x[0]
maxx = x[-1]
if (sx, minx, maxx) not in memo:
memo[(sx, minx, maxx)] = f(x)
return memo[(sx, minx, maxx)]
return memoise_helper
def memoise_2var(f):
"""
Memoisation helper for a function taking 2 numpy arrays as input. As they
are not hashable, we cannot use the standard memoisation function. Here
we use the array sizes, minimum, and maximum values as hashes.
Parameters
----------
f: function
Function to memoise.
Returns
-------
memoise_helper: function
Meomoised best_grid function
"""
memo = {}
def memoise_helper(x, y):
sx = x.size
minx = x[0]
maxx = x[-1]
sy = y.size
miny = y[0]
maxy = y[-1]
if (sx, minx, maxx, sy, miny, maxy) not in memo:
memo[(sx, minx, maxx, sy, miny, maxy)] = f(x, y)
return memo[(sx, minx, maxx, sy, miny, maxy)]
return memoise_helper
@memoise_2var
def best_grid(wavelengths1, wavelengths2):
"""
Return the best wavelength grid to regrid to arrays
Considering the two wavelength grids passed in parameters, this function
compute the best new grid that will be used to regrid the two spectra
before combining them. We do not use np.unique as it is much slowe than
finding the unique elements by hand.
Parameters
----------
wavelengths1, wavelengths2: array of floats
The wavelength grids to be 'regridded'.
Returns
-------
new_grid: array of floats
Array containing all the wavelengths found in the input arrays.
"""
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:])
return wl[flag]
@memoise_2var
def diff_wl(wavelengths1, wavelengths2):
"""
Memoised version of numpy.setdiff1d. It is used to find the values that
are in wavelengths1 but not in wavelengths2.
Parameters
----------
wavelengths1, wavelengths2: array of floats
The wavelength grids to be compared
Returns
-------
setdiff1d: array
Array containing the wavelengths in wavelengths1 but not in
wavelengths2.
"""
return np.setdiff1d(wavelengths1, wavelengths2, assume_unique=True)
@memoise_1var
def argsort_wl(wavelengths):
"""
Memoised version of numpy.argsort. It is used to sort the wavelengths.
Parameters
----------
wavelengths: array
Array containing the wavelengths to be sorted.
Returns
-------
argsort: array
Array containing the sorting indices.
"""
return np.argsort(wavelengths,
kind='mergesort')
def interpolate_lumin(wl, lumin, wl_new, lumin_new):
"""
Procedure to interpolate the luminosity components given the wavelengths
grid and the luminosity of a new component.
Parameters
----------
wl: 1D array
Wavelengths grid of the luminosity components.
lumin: 2D array
Luminosity components. The first index gives a given component and the
second gives the wavelength.
wl_new: 1D array
Wavelengths grid of he new component.
lumin_new: 1D array
Luminosity of the new component.
Returns
-------
wl_best: 1D array
The best merged wavelengths grid
lumin_out: 2D array
New luminosities array resampled if necessary to the best merge
wavelengths grid, and including the new component.
"""
# First we remove from the wavelengths grid of the new component the
# wavelengths that are already present for the other components. We do
# that to avoid having the interpolate on a wavelength for which we
# already know the answer.
wl_unique = diff_wl(wl_new, wl)
# Output 2D array. We increase the number of components by one to include
# the new directly rather than vstack-ing later. We also increase the size
# to include possible new wavelengths
lumin_out = np.zeros((lumin.shape[0]+1, lumin.shape[1]+wl_unique.size))
# If the new component has non already existing wavelengths then we
# interpolate on these new wavelengths and only those ones. The first
# lumin.shape[0] elements are the same as the input luminosity components.
# We only carry out the interpolation of the new components. Finally, we
# reorder the output array with increasing wavelength.
if wl_unique.size > 0:
lumin_out[:-1, :lumin.shape[1]] = lumin
# We interpolate only on the wavelengths where the components are
# 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, :])
wl_best = np.concatenate((wl, wl_unique))
s = argsort_wl(wl_best)
wl_best = wl_best[s]
lumin_out = lumin_out[:, s]
else:
wl_best = wl
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.)
return (wl_best, lumin_out)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment