v1.0.1 - fitting spectra with variable resolution

bagpipes/models/__init__.py

@@ -5,6 +5,7 @@
 from .stellar_model import stellar
 from .dust_emission_model import dust_emission
 from .dust_attenuation_model import dust_attenuation
+from .agn_model import agn
 
 from .nebular_model import nebular
 from .igm_model import igm

bagpipes/models/agn_model.py

@@ -0,0 +1,57 @@
+from __future__ import print_function, division, absolute_import
+
+import numpy as np
+import spectres
+
+from .. import config
+
+
+class agn(object):
+    """ A basic rest-optical AGN continuum + broad line model.
+
+    Parameters
+    ----------
+
+    wavelengths : np.ndarray
+        1D array of wavelength values for the continuum.
+
+    param : dict
+        Parameters for the AGN model.
+    """
+
+    def __init__(self, wavelengths):
+        self.wavelengths = wavelengths
+
+    def update(self, param):
+        self.param = param
+
+        tau = 5000.
+        mask1 = (self.wavelengths < 5000.)
+        mask2 = np.invert(mask1)
+
+        agn_spec = np.zeros_like(self.wavelengths)
+
+        agn_spec[mask1] = self.wavelengths[mask1] ** param["alphalam"]
+        agn_spec[mask2] = self.wavelengths[mask2] ** param["betalam"]
+
+        agn_spec[mask1] /= agn_spec[mask1][-1]
+        agn_spec[mask2] /= agn_spec[mask2][0]
+
+        agn_spec *= param["f5100A"]
+        agn_spec /= agn_spec[np.argmin(np.abs(self.wavelengths - 5100.))]
+
+        agn_spec += self.gaussian_model(4861.35, param["sigma"],
+                                        param["hanorm"]/2.86)
+
+        agn_spec += self.gaussian_model(6562.80, param["sigma"],
+                                        param["hanorm"])
+
+        self.spectrum = agn_spec
+
+    def gaussian_model(self, central_wav, sigma_vel, norm):
+        x = self.wavelengths
+        sigma_aa = sigma_vel*central_wav/(3*10**5)
+        gauss = (norm/(sigma_aa*np.sqrt(2*np.pi)))
+        gauss *= np.exp(-0.5*(x-central_wav)**2/sigma_aa**2)
+
+        return gauss

bagpipes/models/chemical_enrichment_history.py

@@ -24,6 +24,70 @@ def __init__(self, model_comp, sfh_weights):
 
                 self.grid += self.grid_comp[comp]
 
+    def metallicity_bins(self, comp, sfh):
+        bin_edges = np.array(comp["bin_edges"])[::-1]*10**6
+        n_bins = len(bin_edges) - 1
+        ages = config.age_sampling
+        grid = np.zeros((self.zmet_vals.shape[0], sfh.shape[0]))
+
+        for i in range(1, n_bins+1):
+            zmet = comp["metallicity" + str(i)]
+
+            weights = np.zeros(self.zmet_vals.shape[0])
+
+            high_ind = self.zmet_vals[self.zmet_vals < zmet].shape[0]
+
+            if high_ind == self.zmet_vals.shape[0]:
+                weights[-1] = 1.
+
+            elif high_ind == 0:
+                weights[0] = 1.
+
+            else:
+                low_ind = high_ind - 1
+                width = (self.zmet_vals[high_ind] - self.zmet_vals[low_ind])
+                weights[high_ind] = (zmet - self.zmet_vals[low_ind])/width
+                weights[high_ind-1] = 1 - weights[high_ind]
+
+            mask = (ages < bin_edges[i-1]) & (ages > bin_edges[i])
+            grid[:, mask] = np.expand_dims(weights, axis=1)
+
+        return grid*sfh
+
+    def metallicity_bins_continuity(self, comp, sfh):
+        bin_edges = np.array(comp["bin_edges"])[::-1]*10**6
+        n_bins = len(bin_edges) - 1
+        ages = config.age_sampling
+        grid = np.zeros((self.zmet_vals.shape[0], sfh.shape[0]))
+        dzmet = [comp["dzmet" + str(i)] for i in range(1, n_bins)]
+        for i in range(1, n_bins+1):
+
+            zmet = comp["metallicity1"]
+            if i >= 2:
+                zmet = 10**(np.log10(comp["metallicity1"])
+                            + np.sum(dzmet[:i-1]))
+
+            weights = np.zeros(self.zmet_vals.shape[0])
+
+            high_ind = self.zmet_vals[self.zmet_vals < zmet].shape[0]
+
+            if high_ind == self.zmet_vals.shape[0]:
+                weights[-1] = 1.
+
+            elif high_ind == 0:
+                weights[0] = 1.
+
+            else:
+                low_ind = high_ind - 1
+                width = (self.zmet_vals[high_ind] - self.zmet_vals[low_ind])
+                weights[high_ind] = (zmet - self.zmet_vals[low_ind])/width
+                weights[high_ind-1] = 1 - weights[high_ind]
+
+            mask = (ages < bin_edges[i-1]) & (ages > bin_edges[i])
+            grid[:, mask] = np.expand_dims(weights, axis=1)
+
+        return grid*sfh
+
     def delta(self, comp, sfh):
         """ Delta function metallicity history. """
 

bagpipes/models/igm_model.py

@@ -6,12 +6,12 @@
 
 
 def interp_discont(x, xp, fp, xdiscont, left=None, right=None):
-    """Interpolates separately on both sides of a discontinuity (not over it)"""
-    i  = np.searchsorted(x, xdiscont)
+    """Interpolates separately on both sides of a discontinuity, not over it"""
+    i = np.searchsorted(x, xdiscont)
     ip = np.searchsorted(xp, xdiscont)
     y1 = np.interp(x[:i], xp[:ip], fp[:ip], left=left)
     y2 = np.interp(x[i:], xp[ip:], fp[ip:], right=right)
-    y  = np.concatenate([y1, y2])
+    y = np.concatenate([y1, y2])
     return y
 
 
@@ -42,7 +42,7 @@ def _resample_in_wavelength(self):
                                         config.igm_wavelengths,
                                         config.raw_igm_grid[i, :], 1215.67,
                                         left=0., right=1.)
-                                   
+
         # Make sure the pixel containing Lya is always IGM attenuated
         lya_ind = np.abs(self.wavelengths - 1215.67).argmin()
         if self.wavelengths[lya_ind] > 1215.67:

bagpipes/models/model_galaxy.py

@@ -2,6 +2,7 @@
 
 import numpy as np
 import warnings
+import spectres
 
 from copy import deepcopy
 
@@ -15,6 +16,7 @@
 from .dust_attenuation_model import dust_attenuation
 from .nebular_model import nebular
 from .igm_model import igm
+from .agn_model import agn
 from .star_formation_history import star_formation_history
 from ..input.spectral_indices import measure_index
 
@@ -99,9 +101,14 @@ def __init__(self, model_components, filt_list=None, spec_wavs=None,
         self.nebular = False
         self.dust_atten = False
         self.dust_emission = False
+        self.agn = False
 
         if "nebular" in list(model_components):
-            self.nebular = nebular(self.wavelengths)
+            if "velshift" not in model_components["nebular"]:
+                model_components["nebular"]["velshift"] = 0.
+
+            self.nebular = nebular(self.wavelengths,
+                                   model_components["nebular"]["velshift"])
 
             if "metallicity" in list(model_components["nebular"]):
                 self.neb_sfh = star_formation_history(model_components)
@@ -111,6 +118,9 @@ def __init__(self, model_components, filt_list=None, spec_wavs=None,
             self.dust_atten = dust_attenuation(self.wavelengths,
                                                model_components["dust"])
 
+        if "agn" in list(model_components):
+            self.agn = agn(self.wavelengths)
+
         self.update(model_components)
 
     def _get_wavelength_sampling(self):
@@ -179,6 +189,29 @@ def _get_wavelength_sampling(self):
 
         return np.array(x)
 
+    def _get_R_curve_wav_sampling(self, oversample=4):
+        """ Calculate wavelength sampling for the model to be resampled
+        onto in order to apply variable spectral broadening. Only used
+        if a resolution curve is supplied in model_components. Has to
+        be re-run when a model is updated as it depends on redshift.
+
+        Parameters
+        ----------
+
+        oversample : float
+            Number of spectral samples per full width at half maximum.
+        """
+
+        R_curve = self.model_comp["R_curve"]
+        x = [0.95*self.spec_wavs[0]]
+
+        while x[-1] < 1.05*self.spec_wavs[-1]:
+            R_val = np.interp(x[-1], R_curve[:, 0], R_curve[:, 1])
+            dwav = x[-1]/R_val/oversample
+            x.append(x[-1] + dwav)
+
+        return np.array(x)
+
     def _get_index_spec_wavs(self, model_components):
         """ Generate an appropriate spec_wavs array for covering the
         spectral indices specified in index_list. """
@@ -234,13 +267,30 @@ def update(self, model_components):
 
         else:
             self._calculate_full_spectrum(model_components)
-            self._calculate_uvj_mags()
+
+        if self.spec_wavs is not None:
+            self._calculate_spectrum(model_components)
+
+        # Add any AGN component:
+        if self.agn:
+            self.agn.update(self.model_comp["agn"])
+            agn_spec = self.agn.spectrum
+            agn_spec *= self.igm.trans(self.model_comp["redshift"])
+
+            self.spectrum_full += agn_spec/(1. + self.model_comp["redshift"])
+
+            if self.spec_wavs is not None:
+                zplus1 = (self.model_comp["redshift"] + 1.)
+                agn_interp = np.interp(self.spec_wavs, self.wavelengths*zplus1,
+                                       agn_spec/zplus1, left=0, right=0)
+
+                self.spectrum[:, 1] += agn_interp
 
         if self.filt_list is not None:
             self._calculate_photometry(model_components["redshift"])
 
-        if self.spec_wavs is not None:
-            self._calculate_spectrum(model_components)
+        if not self.sfh.unphysical:
+            self._calculate_uvj_mags()
 
         # Deal with any spectral index calculations.
         if self.index_list is not None:
@@ -358,10 +408,11 @@ def _calculate_full_spectrum(self, model_comp):
 
         if self.dust_atten:
             self.spectrum_bc *= 3.826*10**33
-            
+
         em_lines *= 3.826*10**33
 
         self.line_fluxes = dict(zip(config.line_names, em_lines))
+
         self.spectrum_full = spectrum
 
     def _calculate_photometry(self, redshift, uvj=False):
@@ -413,9 +464,33 @@ def _calculate_spectrum(self, model_comp):
             spectrum = self.spectrum_full
             redshifted_wavs = zplusone*self.wavelengths
 
-        fluxes = np.interp(self.spec_wavs,
-                           redshifted_wavs,
-                           spectrum, left=0, right=0)
+        if "R_curve" in list(model_comp):
+            oversample = 4  # Number of samples per FWHM at resolution R
+            new_wavs = self._get_R_curve_wav_sampling(oversample=oversample)
+
+            # spectrum = np.interp(new_wavs, redshifted_wavs, spectrum)
+            spectrum = spectres.spectres(new_wavs, redshifted_wavs,
+                                         spectrum, fill=0)
+            redshifted_wavs = new_wavs
+
+            sigma_pix = oversample/2.35  # sigma width of kernel in pixels
+            k_size = 4*int(sigma_pix+1)
+            x_kernel_pix = np.arange(-k_size, k_size+1)
+
+            kernel = np.exp(-(x_kernel_pix**2)/(2*sigma_pix**2))
+            kernel /= np.trapz(kernel)  # Explicitly normalise kernel
+
+            # Disperse non-uniformly sampled spectrum
+            spectrum = np.convolve(spectrum, kernel, mode="valid")
+            redshifted_wavs = redshifted_wavs[k_size:-k_size]
+
+        # Converted to using spectres in response to issue with interp,
+        # see https://github.com/ACCarnall/bagpipes/issues/15
+        # fluxes = np.interp(self.spec_wavs, redshifted_wavs,
+        #                    spectrum, left=0, right=0)
+
+        fluxes = spectres.spectres(self.spec_wavs, redshifted_wavs,
+                                   spectrum, fill=0)
 
         if self.spec_units == "mujy":
             fluxes /= ((10**-29*2.9979*10**18/self.spec_wavs**2))

bagpipes/models/nebular_model.py

@@ -19,9 +19,9 @@ class nebular(object):
         1D array of wavelength values desired for the stellar models.
     """
 
-    def __init__(self, wavelengths):
+    def __init__(self, wavelengths, velshift):
         self.wavelengths = wavelengths
-
+        self.velshift = velshift
         self.combined_grid, self.line_grid = self._setup_grids()
 
     def _setup_grids(self):
@@ -57,7 +57,8 @@ def _setup_grids(self):
 
         # Add the nebular lines to the resampled nebular continuum grid.
         for i in range(config.line_wavs.shape[0]):
-            ind = np.abs(self.wavelengths - config.line_wavs[i]).argmin()
+            line_wav_shift = config.line_wavs[i]*(1+(self.velshift/(3*10**5)))
+            ind = np.abs(self.wavelengths - line_wav_shift).argmin()
             if ind != 0 and ind != self.wavelengths.shape[0]-1:
                 width = (self.wavelengths[ind+1] - self.wavelengths[ind-1])/2
                 comb_grid[ind, :, :, :] += line_grid[i, :, :, :]/width

bagpipes/models/star_formation_history.py

@@ -134,9 +134,11 @@ def _calculate_derived_quantities(self):
         self.mass_weighted_age = np.sum(self.sfh*self.age_widths*self.ages)
         self.mass_weighted_age /= np.sum(self.sfh*self.age_widths)
 
-        self.mass_weighted_metallicity = np.sum(self.live_frac_grid*self.ceh.grid, axis=1)
-        self.mass_weighted_metallicity /= np.sum(self.live_frac_grid*self.ceh.grid)
-        self.mass_weighted_metallicity = np.sum(self.mass_weighted_metallicity*config.metallicities)
+        self.mass_weighted_met = np.sum(self.live_frac_grid*self.ceh.grid,
+                                        axis=1)
+        self.mass_weighted_met /= np.sum(self.live_frac_grid*self.ceh.grid)
+        self.mass_weighted_met *= config.metallicities
+        self.mass_weighted_met = np.sum(self.mass_weighted_met)
 
         self.tform = self.age_of_universe - self.mass_weighted_age
 
@@ -169,11 +171,13 @@ def _resample_live_frac_grid(self):
                                                   raw_live_frac_grid[:, i])
 
     def massformed_at_redshift(self, redshift):
-        t_hubble_at_z = 10**9*np.interp(redshift, utils.z_array, utils.age_at_z)
+        t_hubble_at_z = np.interp(redshift, utils.z_array, utils.age_at_z)
+        t_hubble_at_z *= 10**9
 
         mass_assembly = np.cumsum(self.sfh[::-1]*self.age_widths[::-1])[::-1]
 
-        ind = np.argmin(np.abs(self.ages - (self.age_of_universe - t_hubble_at_z)))
+        indices = np.abs(self.ages - (self.age_of_universe - t_hubble_at_z))
+        ind = np.argmin(indices)
 
         return np.log10(mass_assembly[ind])
 

docs/index.rst

@@ -16,7 +16,7 @@ Star-formation history recovery from spectroscopy (see `Carnall et al. 2019b <ht
 
 Identification of z > 3 quiescent galaxies from photometry (see `Carnall et al. 2020 <https://arxiv.org/abs/2001.11975>`_)
 
-Bagpipes has been used in ~30 refereed publications as of Jan 2021. For more example use cases take a look at `papers that cite the bagpipes paper <https://ui.adsabs.harvard.edu/abs/2018MNRAS.480.4379C/citations>`_, e.g. `Carnall et al. (2019a) <https://arxiv.org/abs/1811.03635>`_, `Williams et al. (2019) <https://arxiv.org/abs/1905.11996>`_ and `Wild et al. (2020) <https://arxiv.org/abs/2001.09154>`_.
+Bagpipes has been used in ~130 papers as of April 2023. For more example use cases take a look at this `ADS library listing papers that use Bagpipes <https://ui.adsabs.harvard.edu/public-libraries/VOrR8ITjTTSYNXVYiQ1oag>`_.
 
 
 Source and installation