Commit 36a3b9ec authored by Médéric Boquien's avatar Médéric Boquien

We do not store the time grid in the SED anymore given that we assume it...

We do not store the time grid in the SED anymore given that we assume it starts at 0 Myr with steps of 1 Myr, we can easily reconstruct to save it if needed. It should save a little bit of memory and it should go a little bit faster.
parent 53b3a832
......@@ -34,6 +34,7 @@
### 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)
- The sfhfromfile module is now fully initialised when it is instantiated rather than doing so when processing the SED. This should be especially sensitive when processing different SED. (Médéric Boquien)
- We do not store the time grid in the SED anymore given that we assume it starts at 0 Myr with steps of 1 Myr, we can easily reconstruct to save it if needed. It should save a little bit of memory and it should go a little bit faster. (Médéric Boquien)
## 0.8.1 (2015-12-07)
### Fixed
......
......@@ -61,23 +61,16 @@ class BC03(object):
self.color_table = color_table
self.lumin_table = lumin_table
def convolve(self, sfh_time, sfh_sfr):
def convolve(self, sfh):
"""Convolve the SSP with a Star Formation History
Given a SFH (an time grid and the corresponding star formation rate
SFR), this method convolves the color table and the SSP luminosity
spectrum along the whole SFR.
The time grid of the SFH is expected to be ordered and must not run
beyong 14 Gyr (the maximum time in the database).
Given an SFH, this method convolves the color table and the SSP
luminosity spectrum.
Parameters
----------
sfh_time: array of floats
Time grid [Myr] of the star formation history. It must start at 0
and not run beyond 14 Gyr, with a step of 1 Myr.
sfh_sfr: array of floats
Star Formation Rates in Msun/yr at each time of the SFH time grid.
sfh: array of floats
Star Formation History in Msun/yr.
Returns
-------
......@@ -97,20 +90,15 @@ class BC03(object):
- "b_912": Amplitude of Lyman discontinuity
"""
# We make sure the time grid starts from 0. Otherwise the selection of
# the SSP will be wrong.
if sfh_time[0] != 0:
raise Exception("Age grid must start from 0. Exiting.")
# The convolution is just a matter of reverting the SFH and computing
# the sum of the data from the SSP one to one product. This is done
# using the dot product.
color_table = self.color_table[:, :sfh_time[-1] + 1]
lumin_table = self.lumin_table[:, :sfh_time[-1] + 1]
color_table = self.color_table[:, :sfh.size]
lumin_table = self.lumin_table[:, :sfh.size]
# The 1.e6 factor is because the SFH is in solar mass per year.
color_info = 1.e6 * np.dot(color_table, sfh_sfr[::-1])
luminosity = 1.e6 * np.dot(lumin_table, sfh_sfr[::-1])
color_info = 1.e6 * np.dot(color_table, sfh[::-1])
luminosity = 1.e6 * np.dot(lumin_table, sfh[::-1])
bc03_info = dict(zip(
["m_star", "m_gas", "n_ly", "b_4000", "b4_vn", "b4_sdss", "b_912"],
......
......@@ -70,23 +70,16 @@ class M2005(object):
self.mass_table = mass_table
self.spec_table = spec_table
def convolve(self, sfh_time, sfh_sfr):
def convolve(self, sfh):
"""Convolve the SSP with a Star Formation History
Given a SFH (an time grid and the corresponding star formation rate
SFR), this method convolves the mass distribution and the SSP spectrum
along the whole SFR.
The time grid of the SFH is expected to be ordered and must not run
beyond 13.7 Gyr (the maximum time for Maraston 2005 SSP).
Given an SFH, this method convolves the masses and the spectra of the
SSPs.
Parameters
----------
sfh_time: array of floats
Time grid [Myr)] of the star formation history. It must start at 0
and not run beyond 13.7 Gyr, with a step of 1 Myr.
sfh_sfr: array of floats
Star Formation Rates in Msun/yr at each time of the SFH time grid.
sfh: array of floats
Star Formation History in Msun/yr.
Returns
-------
......@@ -104,20 +97,15 @@ class M2005(object):
- spectra[1]: luminosity in W/nm
"""
# We make sure the time grid starts from 0. Otherwise the selection of
# the SSP will be wrong.
if sfh_time[0] != 0:
raise Exception("Age grid must start from 0. Exiting.")
# As both the SFH and the SSP (limited to the age of the SFH) data now
# share the same time grid, the convolution is just a matter of
# reverting one and computing the sum of the one to one product; this
# is done using the dot product.
mass_table = self.mass_table[:, :sfh_time[-1] + 1]
spec_table = self.spec_table[:, :sfh_time[-1] + 1]
mass_table = self.mass_table[:, :sfh.size]
spec_table = self.spec_table[:, :sfh.size]
# The 1.e6 factor is because the SFH is in solar mass per year.
masses = 1.e6 * np.dot(mass_table, sfh_sfr[::-1])
spectra = 1.e6 * np.dot(spec_table, sfh_sfr[::-1])
masses = 1.e6 * np.dot(mass_table, sfh[::-1])
spectra = 1.e6 * np.dot(spec_table, sfh[::-1])
return masses, spectra
......@@ -7,8 +7,7 @@
This class represents a Spectral Energy Distribution (SED) as used by pcigale.
Such SED is characterised by:
- sfh: a tuple (time [Myr], Star Formation Rate [Msun/yr]) representing the
Star Formation History of the galaxy.
- sfh: the Star Formation History of the galaxy.
- modules: a list of tuples (module name, parameter dictionary) containing all
the pcigale modules the SED 'went through'.
......@@ -84,15 +83,15 @@ class SED(object):
# it's needed.
self._sfh = value
if value:
sfh_time, sfh_sfr = value
if value is not None:
sfh_sfr = value
self._sfh = value
self.add_info("sfh.sfr", sfh_sfr[-1], True, force=True)
self.add_info("sfh.sfr10Myrs", np.mean(sfh_sfr[-10:]), True,
force=True)
self.add_info("sfh.sfr100Myrs", np.mean(sfh_sfr[-100:]), True,
force=True)
self.add_info("sfh.age", sfh_time[-1]+1, False, force=True)
self.add_info("sfh.age", sfh_sfr.size, False, force=True)
@property
def fnu(self):
......@@ -364,7 +363,7 @@ class SED(object):
"""
sed = SED()
if self._sfh is not None:
sed._sfh = (self._sfh[0], self._sfh[1])
sed._sfh = self._sfh
sed.modules = self.modules[:]
if self.wavelength_grid is not None:
sed.wavelength_grid = self.wavelength_grid.copy()
......
......@@ -6,6 +6,7 @@
from collections import OrderedDict
from astropy.table import Table, Column
import numpy as np
def save_sed_to_fits(sed, prefix, norm=1.):
"""
......@@ -38,6 +39,6 @@ def save_sed_to_fits(sed, prefix, norm=1.):
table.write("{}_best_model.fits".format(prefix))
table = Table(meta=info)
table["time"] = Column(sed.sfh[0], unit="Myr")
table["SFR"] = Column(norm * sed.sfh[1], unit="Msun/yr")
table["time"] = Column(np.arange(sed.sfh.size), unit="Myr")
table["SFR"] = Column(norm * sed.sfh, unit="Msun/yr")
table.write("{}_SFH.fits".format(prefix))
......@@ -4,7 +4,7 @@
# Author: Yannick Roehlly
from astropy.io.votable.tree import VOTableFile, Resource, Table, Field, Info
import numpy as np
def save_sed_to_vo(sed, filename, norm=1.):
"""
......@@ -68,9 +68,9 @@ def save_sed_to_vo(sed, filename, norm=1.):
Field(votable, name="SFR", datatype="double", unit="Msun/yr",
ucd="phys.SFR")
])
sfh_table.create_arrays(len(sed.sfh[0]))
sfh_table.array["time"] = sed.sfh[0]
sfh_table.array["SFR"] = norm * sed.sfh[1]
sfh_table.create_arrays(sed.sfh.size)
sfh_table.array["time"] = np.arange(sed.sfh.size)
sfh_table.array["SFR"] = norm * sed.sfh
# SED information to keywords
if sed.sfh is not None:
......
......@@ -73,23 +73,17 @@ class BC03(SedModule):
SED object.
"""
sfh_time, sfh_sfr = sed.sfh
ssp = self.ssp
# Age of the galaxy at each time of the SFH
sfh_age = sfh_time[::-1]
# First, we process the young population (age lower than the
# separation age.)
young_sfh = np.copy(sfh_sfr)
young_sfh[sfh_age > self.separation_age - 1] = 0
young_wave, young_lumin, young_info = ssp.convolve(sfh_time, young_sfh)
young_sfh = np.copy(sed.sfh)
young_sfh[:-self.separation_age] = 0.
young_wave, young_lumin, young_info = self.ssp.convolve(young_sfh)
# Then, we process the old population. If the SFH is shorter than the
# separation age then all the arrays will consist only of 0.
old_sfh = np.copy(sfh_sfr)
old_sfh[sfh_age <= self.separation_age - 1] = 0
old_wave, old_lumin, old_info = ssp.convolve(sfh_time, old_sfh)
old_sfh = np.copy(sed.sfh)
old_sfh[-self.separation_age:] = 0.
old_wave, old_lumin, old_info = self.ssp.convolve(old_sfh)
# We compute the Lyman continuum luminosity as it is important to
# compute the energy absorbed by the dust before ionising gas.
......
......@@ -85,23 +85,17 @@ class M2005(SedModule):
SED object.
"""
sfh_time, sfh_sfr = sed.sfh
ssp = self.ssp
# Age of the galaxy at each time of the SFH
sfh_age = sfh_time[::-1]
# First, we process the young population (age lower than the
# separation age.)
young_sfh = np.copy(sfh_sfr)
young_sfh[sfh_age > self.separation_age] = 0
young_masses, young_spectrum = ssp.convolve(sfh_time, young_sfh)
young_sfh = np.copy(sed.sfh)
young_sfh[:-self.separation_age] = 0.
young_masses, young_spectrum = self.ssp.convolve(sfh_time, young_sfh)
# Then, we process the old population. If the SFH is shorter than the
# separation age then all the arrays will consist only of 0.
old_sfh = np.copy(sfh_sfr)
old_sfh[sfh_age <= self.separation_age] = 0
old_masses, old_spectrum = ssp.convolve(sfh_time, old_sfh)
old_sfh = np.copy(sed.sfh)
old_sfh[-self.separation_age:] = 0.
old_masses, old_spectrum = self.ssp.convolve(sfh_time, old_sfh)
sed.add_module(self.name, self.parameters)
......
......@@ -76,11 +76,11 @@ class Sfh2Exp(SedModule):
normalise = bool(self.parameters["normalise"])
# Time grid and age. If needed, the age is rounded to the inferior Myr
self.time_grid = np.arange(age)
time_grid = np.arange(age)
time_grid_burst = np.arange(self.burst_age)
# SFR for each component
self.sfr = np.exp(-self.time_grid / self.tau_main)
self.sfr = np.exp(-time_grid / self.tau_main)
sfr_burst = np.exp(-time_grid_burst / self.tau_burst)
# Height of the late burst to have the desired produced mass fraction
......@@ -112,7 +112,7 @@ class Sfh2Exp(SedModule):
sed.add_module(self.name, self.parameters)
# Add the sfh and the output parameters to the SED.
sed.sfh = (self.time_grid, self.sfr)
sed.sfh = self.sfr
sed.add_info("sfh.integrated", self.sfr_integrated, True)
sed.add_info("sfh.tau_main", self.tau_main)
sed.add_info("sfh.tau_burst", self.tau_burst)
......
......@@ -63,7 +63,7 @@ class SfhBuat08(SedModule):
normalise = bool(self.parameters["normalise"])
# Time grid and age. If needed, the age is rounded to the inferior Myr
self.time_grid = np.arange(self.age)
time_grid = np.arange(self.age)
# Values from Buat et al. (2008) table 2
paper_velocities = np.array([80., 150., 220., 290., 360.])
......@@ -77,7 +77,7 @@ class SfhBuat08(SedModule):
c = np.interp(self.velocity, paper_velocities, paper_cs)
# Main SFR
t = (self.time_grid+1) / 1000 # The time is in Gyr in the formulae
t = (time_grid+1) / 1000 # The time is in Gyr in the formulae
self.sfr = 10.**(a + b * np.log10(t) + c * t**.5) / 1.e9
# Compute the galaxy mass and normalise the SFH to 1 solar mass
......@@ -97,7 +97,7 @@ class SfhBuat08(SedModule):
sed.add_module(self.name, self.parameters)
# Add the sfh and the output parameters to the SED.
sed.sfh = (self.time_grid, self.sfr)
sed.sfh = self.sfr
sed.add_info("sfh.integrated", self.sfr_integrated, True)
sed.add_info("sfh.velocity", self.velocity)
......
......@@ -62,9 +62,9 @@ class SfhQuench(SedModule):
"""
# Read the star formation history of the SED
time, sfr = sed.sfh
sfr = sed.sfh
if self.quenching_age > time[-1]:
if self.quenching_age > sfr.size:
raise Exception("[sfh_quenching] The quenching age is greater "
"than the galaxy age. Please fix your parameters.")
......@@ -74,10 +74,8 @@ class SfhQuench(SedModule):
# We make the computation only if the quenching age and the quenching
# factor are not 0.
if self.quenching_age > 0 and self.quenching_factor > 0.:
age_lapse = time[1] - time[0]
quenching_idx = np.int(self.quenching_age/age_lapse)
sfr[-quenching_idx:] = sfr[-quenching_idx] * (
1 - self.quenching_factor)
sfr[-self.quenching_age:] = sfr[-quenching_idx] * (
1. - self.quenching_factor)
# Compute the galaxy mass and normalise the SFH to 1 solar mass
# produced if asked to.
......@@ -86,7 +84,7 @@ class SfhQuench(SedModule):
sfr /= sfr_integrated
sfr_integrated = 1.
sed.sfh = (time, sfr)
sed.sfh = sfr
sed.add_info("sfh.integrated", sfr_integrated, True, force=True)
sed.add_module(self.name, self.parameters)
......
......@@ -59,8 +59,8 @@ class SFHDelayed(SedModule):
normalise = bool(self.parameters["normalise"])
# Time grid and SFR
self.time_grid = np.arange(age)
self.sfr = self.time_grid * np.exp(-self.time_grid / self.tau_main) / \
time_grid = np.arange(age)
self.sfr = time_grid * np.exp(-time_grid / self.tau_main) / \
self.tau_main**2
# Compute the galaxy mass and normalise the SFH to 1 solar mass
......@@ -83,7 +83,7 @@ class SFHDelayed(SedModule):
sed.add_module(self.name, self.parameters)
# Add the sfh and the output parameters to the SED.
sed.sfh = (self.time_grid, self.sfr)
sed.sfh = self.sfr
sed.add_info("sfh.integrated", self.sfr_integrated, True)
sed.add_info("sfh.tau_main", self.tau_main)
......
......@@ -65,16 +65,16 @@ class SfhFromFile(SedModule):
table = read_table(filename)
self.sfr = table.columns[self.sfr_column_number].data.astype(np.float)
self.time_grid = table.columns[0].data.astype(np.int)
if self.time_grid[0] != 0:
time_grid = table.columns[0].data.astype(np.int)
if time_grid[0] != 0:
raise Exception("The time grid must start from 0.")
if np.all(self.time_grid[1:]-self.time_grid[:-1] == 1) == False:
if np.all(time_grid[1:]-time_grid[:-1] == 1) == False:
raise Exception("The time step must be 1 Myr. Computed models will"
" be wrong.")
# We cut the SFH to the desired age.
self.sfr = self.sfr[self.time_grid <= age]
self.time_grid = self.time_grid[self.time_grid <= age]
self.sfr = self.sfr[time_grid <= age]
time_grid = time_grid[time_grid <= age]
# Compute the galaxy mass and normalise the SFH to 1 solar mass
# produced if asked to.
......@@ -94,7 +94,7 @@ class SfhFromFile(SedModule):
"""
sed.add_module(self.name, self.parameters)
sed.sfh = (self.time_grid, self.sfr)
sed.sfh = self.sfr
sed.add_info("sfh.integrated", self.sfr_integrated, True)
sed.add_info("sfh.index", self.sfr_column_number)
......
......@@ -75,16 +75,16 @@ class SfhPeriodic(SedModule):
sfr_A = float(self.parameters["sfr_A"])
normalise = bool(self.parameters["normalise"])
self.time_grid = np.arange(0, age)
self.sfr = np.zeros_like(self.time_grid, dtype=np.float)
time_grid = np.arange(0, age)
self.sfr = np.zeros_like(time_grid, dtype=np.float)
if self.type_bursts == 0:
burst = np.exp(-self.time_grid/self.tau_bursts)
burst = np.exp(-time_grid/self.tau_bursts)
elif self.type_bursts == 1:
burst = np.exp(-self.time_grid/self.tau_bursts) * \
self.time_grid/self.tau_bursts**2
burst = np.exp(-time_grid/self.tau_bursts) * \
time_grid/self.tau_bursts**2
elif self.type_bursts == 2:
burst = np.zeros_like(self.time_grid)
burst = np.zeros_like(time_grid)
burst[:self.tau_bursts+1] = 1.
else:
raise Exception("Burst type {} unknown.".format(self.type_bursts))
......@@ -116,7 +116,7 @@ class SfhPeriodic(SedModule):
sed.add_module(self.name, self.parameters)
sed.sfh = (self.time_grid, self.sfr)
sed.sfh = self.sfr
sed.add_info("sfh.integrated", self.sfr_integrated, True)
sed.add_info("sfh.type_bursts", self.type_bursts)
sed.add_info("sfh.delta_bursts", self.delta_bursts)
......
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