Added FAST_SPEC_CONVOLVE

CHANGELOG

@@ -6,6 +6,7 @@ Features:
  - The spectrum no longer inherits the same template error function as the photometry
  - AUTO_SCALE=1 no longer automatically disables the template error for the spectrum
  - Added option SPEC_LSF_FILE to apply a line spread function to the models
+ - Added option FAST_SPEC_CONVOLVE for faster line spread function convolution
 
 Bug fixes:
  - Fixed crash when saving extra outputs and using an absolute path in CATALOG

README.md

@@ -491,6 +491,12 @@ NB: applying the velocity dispersion needs to be done for each SED library that
 ## Line spread function
 ```SPEC_LSF_FILE``` can be used to specify a file holding the line spread function of the instrument that was used to acquire the spectra. This file must contain the "sigma" of the Gaussian profile of the line spread function, expressed in Angstroms (observer frame), as a function of observed wavelength (also in Angstroms). It will be linearly interpolated at each observed wavelength of the models used in the fit, and the resulting value will be used to broaden the model by the specified value around that observed wavelength.
 
+This is a costly process since it must be done for each model. If your model grid has a narrow redshift range, an approximation can be enbaled with ```FAST_SPEC_CONVOLVE=1``` (default is ```0```), in which case the LSF will be applied to the rest-frame templates, pre-attenuation & IGM absorption, instead of observer-frame templates. This is not quite correct, since the LSF spectral broadening is a phenomenon happening in the instrument and not to the stellar population, but in practice this is a very good approximation provided that:
+ - the range of redshifts in the grid is small (e.g., max delta_z/z = 0.05)
+ - the LSF varies smoothly as a function of wavelength
+
+If using ```FAST_SPEC_CONVOLVE=1```, please make sure to test its impact for your specific use case before using the results.
+
 
 ## Continuum indices
 Traditionally, absorption line strengths are measured as "equivalent widths" (in wavelength unit). These are part of a large set of continuum features call "spectral indices", see for example [Balogh et al. (1999)](http://adsabs.harvard.edu/abs/1999ApJ...527...54B), which includes other commonly used quantities such as the Dn4000 index. These features have been used extensively in the past to characterize the SEDs of galaxies and infer their star formation histories and/or metal abundances.

example/cosmos-20115/fast.param

@@ -192,13 +192,19 @@ USE_LIR             = 1          # 0 / 1
 #   resolution of is finer than the LSF, and if binning the spectrum is
 #   not acceptable.
 #
+# o FAST_SPEC_CONVOLVE: Set to 1 to enable a speed optimisation to
+#   convolve the line spread function. This optimisation is only valid
+#   if the LSF varies smoothly, and the redshift range of the grid is
+#   small.
+#
 #-----------------------------------------------------------------------
 
 SPECTRUM            = 'spectrum'
 TEMP_ERR_SPEC_FILE  = ''
 AUTO_SCALE          = 1          # 0 / 1
 APPLY_VDISP         = 350        # km/s
 SPEC_LSF_FILE       = ''
+FAST_SPEC_CONVOLVE  = 0
 
 
 #--- OUTPUT INFORMATION  -----------------------------------------------

example/hdfn_fs99/fast.param

@@ -192,13 +192,19 @@ USE_LIR             = 0          # 0 / 1
 #   resolution of is finer than the LSF, and if binning the spectrum is
 #   not acceptable.
 #
+# o FAST_SPEC_CONVOLVE: Set to 1 to enable a speed optimisation to
+#   convolve the line spread function. This optimisation is only valid
+#   if the LSF varies smoothly, and the redshift range of the grid is
+#   small.
+#
 #-----------------------------------------------------------------------
 
 SPECTRUM            = ''
 TEMP_ERR_SPEC_FILE  = ''
 AUTO_SCALE          = 0          # 0 / 1
 APPLY_VDISP         = 0          # km/s
 SPEC_LSF_FILE       = ''
+FAST_SPEC_CONVOLVE  = 0
 
 
 #--- OUTPUT INFORMATION  -----------------------------------------------

src/fast++-gridder-custom.cpp

@@ -118,9 +118,7 @@ bool gridder_t::build_and_send_custom(fitter_t& fitter) {
         }
 
         // Apply velocity dispersion
-        if (is_finite(opts.apply_vdisp)) {
-            ssp.sed = convolve_rest(ssp.lambda, ssp.sed);
-        }
+        convolve_rest(ssp.lambda, ssp.sed);
 
         // Pre-compute dust law & IGM absorption (they don't change with SFH)
         vec1d dust_law = build_dust_law(ssp.lambda);
@@ -285,9 +283,7 @@ bool gridder_t::build_template_custom(uint_t iflat, vec1f& lam, vec1f& flux) con
         cached_library = filename;
 
         // Apply velocity dispersion
-        if (is_finite(opts.apply_vdisp)) {
-            ssp->sed = convolve_rest(ssp->lambda, ssp->sed);
-        }
+        convolve_rest(ssp->lambda, ssp->sed);
     }
 
     // Integrate SFH on local time grid
@@ -339,9 +335,7 @@ bool gridder_t::build_template_custom(uint_t iflat, vec1f& lam, vec1f& flux_youn
         cached_library = filename;
 
         // Apply velocity dispersion
-        if (is_finite(opts.apply_vdisp)) {
-            ssp->sed = convolve_rest(ssp->lambda, ssp->sed);
-        }
+        convolve_rest(ssp->lambda, ssp->sed);
     }
 
     // Integrate SFH for template

src/fast++-gridder-ised.cpp

@@ -307,9 +307,7 @@ bool gridder_t::build_and_send_ised(fitter_t& fitter) {
         }
 
         // Apply velocity dispersion
-        if (is_finite(opts.apply_vdisp)) {
-            ised.fluxes = convolve_rest(ised.lambda, ised.fluxes);
-        }
+        convolve_rest(ised.lambda, ised.fluxes);
 
         // Pre-compute dust law & IGM absorption (they don't change with SFH)
         vec1d dust_law = build_dust_law(ised.lambda);
@@ -433,9 +431,7 @@ bool gridder_t::build_template_ised(uint_t iflat, vec1f& lam, vec1f& flux) const
             cached_library = filename;
 
             // Apply velocity dispersion
-            if (is_finite(opts.apply_vdisp)) {
-                ised->fluxes = convolve_rest(ised->lambda, ised->fluxes);
-            }
+            convolve_rest(ised->lambda, ised->fluxes);
         }
     }
 

src/fast++-gridder.cpp

@@ -189,6 +189,13 @@ gridder_t::gridder_t(const options_t& opt, const input_state_t& inp, output_stat
         output_z = max(cz, 0.00001);
     }
 
+    if (opts.fast_spec_convolve) {
+        double dz = (max(output_z) - min(output_z))/(1.0 + mean(output_z));
+        if (dz > 0.05) {
+            warning("delta_z/(1+z) is large (", dz, "); consider turning off FAST_SPEC_CONVOLVE");
+        }
+    }
+
     // Pre-compute distances (galaxev templates are in Lsun/A)
     {
         const double dist_Mpc_to_cgs = 3.0856e24; // [cm/Mpc]
@@ -702,10 +709,8 @@ void gridder_t::redshift_and_send(fitter_t& fitter, progress_t& pg,
             tpl_att_z_lam.safe[il] *= (1.0 + output_z.safe[iz]);
         }
 
-        if (!input.tpllsf_sigma.empty()) {
-            // Apply velocity dispersion
-            tpl_att_z_flux = convolve_obs(tpl_att_z_lam, tpl_att_z_flux);
-        }
+        // Apply LSF
+        convolve_obs(tpl_att_z_lam, tpl_att_z_flux);
 
         // Integrate
         for (uint_t il : range(input.lambda)) {
@@ -1021,10 +1026,8 @@ bool gridder_t::build_template_impl(uint_t iflat, bool nodust,
         lam.safe[il] *= (1.0 + z);
     }
 
-    if (!input.tpllsf_sigma.empty()) {
-        // Apply velocity dispersion
-        flux = convolve_obs(lam, flux);
-    }
+    // Apply LSF
+    convolve_obs(lam, flux);
 
     // Integrate
     iflux.resize(input.lambda.size());
@@ -1117,26 +1120,55 @@ vec2d gridder_t::convolve_function(const vec1d& lam, const vec2d& osed, uint_t n
     return sed;
 }
 
-vec1d gridder_t::convolve_obs(const vec1d& lam, const vec1d& osed) const {
+void gridder_t::convolve_obs(const vec1f& lam, vec1f& osed) const {
     // We can't use multi-threading if 'generators' is chosen as parallel mode, since
     // we will already be generating multiple models in different threads.
     uint_t nthread = (opts.parallel == parallel_choice::generators ? 1 : opts.n_thread);
 
-    return flatten(convolve_function(lam, reform(osed, 1, osed.size()), nthread, [&](double lobs) {
-        if (lobs < input.tpllsf_lam.front()) {
-            return input.tpllsf_sigma.front();
-        } else if (lobs > input.tpllsf_lam.back()) {
-            return input.tpllsf_sigma.back();
-        } else {
-            return interpolate(input.tpllsf_sigma, input.tpllsf_lam, lobs);
+    if (!input.tpllsf_sigma.empty() && !opts.fast_spec_convolve) {
+        osed = flatten(convolve_function(lam, reform(osed, 1, osed.size()), nthread, [&](double lobs) {
+            if (lobs < input.tpllsf_lam.front()) {
+                return input.tpllsf_sigma.front();
+            } else if (lobs > input.tpllsf_lam.back()) {
+                return input.tpllsf_sigma.back();
+            } else {
+                return interpolate(input.tpllsf_sigma, input.tpllsf_lam, lobs);
+            }
+        }));
+    }
+}
+
+void gridder_t::convolve_rest(const vec1d& lam, vec2d& osed) const {
+    if (!input.tpllsf_sigma.empty() && opts.fast_spec_convolve) {
+        // Approximation: convolve LSF + velocity dispersion in one go for speed
+        double vdisp = is_finite(opts.apply_vdisp) ? opts.apply_vdisp : 0.0;
+        double z = mean(output.grid[grid_id::z]);
+        osed = convolve_function(lam, osed, opts.n_thread, [&](double lrest) {
+            double lobs = lrest*(1.0 + z);
+            double lsf = 0.0;
+            if (lobs < input.tpllsf_lam.front()) {
+                lsf = input.tpllsf_sigma.front();
+            } else if (lobs > input.tpllsf_lam.back()) {
+                lsf = input.tpllsf_sigma.back();
+            } else {
+                lsf = interpolate(input.tpllsf_sigma, input.tpllsf_lam, lobs);
+            }
+
+            return std::sqrt(sqr(lsf/(1.0 + z)) + sqr(lrest*vdisp/2.99792e5));
+        });
+    } else {
+        if (is_finite(opts.apply_vdisp)) {
+            osed = convolve_function(lam, osed, opts.n_thread, [&](double lrest) {
+                return lrest*opts.apply_vdisp/2.99792e5;
+            });
         }
-    }));
+    }
 }
 
-vec2d gridder_t::convolve_rest(const vec1d& lam, const vec2d& osed) const {
-    return convolve_function(lam, osed, opts.n_thread, [&](double lrest) {
-        return lrest*opts.apply_vdisp/2.99792e5;
-    });
+void gridder_t::convolve_rest(const vec1f& lam, vec2f& osed) const {
+    vec2d dsed = osed;
+    convolve_rest(lam, dsed);
+    osed = dsed;
 }
 
 struct sed_id_pair {

src/fast++-read_input.cpp

@@ -97,6 +97,7 @@ bool read_params(options_t& opts, input_state_t& state, const std::string& filen
         PARSE_OPTION(temp_err_file)
         PARSE_OPTION(temp_err_spec_file)
         PARSE_OPTION(spec_lsf_file)
+        PARSE_OPTION(fast_spec_convolve)
         PARSE_OPTION(name_zphot)
         PARSE_OPTION(spectrum)
         PARSE_OPTION(auto_scale)

src/fast++.hpp

@@ -101,6 +101,7 @@ struct options_t {
     // Spectra
     std::string temp_err_spec_file;
     std::string spec_lsf_file;
+    bool fast_spec_convolve = false;
 
     // z-phot behavior
     bool force_zphot = false;
@@ -419,8 +420,9 @@ private :
 
     vec2d convolve_function(const vec1d& lam, const vec2d& osed, uint_t n_thread,
         std::function<double(double)> sigma_fun) const;
-    vec1d convolve_obs(const vec1d& lam, const vec1d& osed) const;
-    vec2d convolve_rest(const vec1d& lam, const vec2d& osed) const;
+    void convolve_obs(const vec1f& lam, vec1f& osed) const;
+    void convolve_rest(const vec1d& lam, vec2d& osed) const;
+    void convolve_rest(const vec1f& lam, vec2f& osed) const;
 
     mutable tinyexpr_wrapper sfh_expr;
     mutable tinyexpr_wrapper exclude_expr;