From d1adef88eda1fbd5c9a2109dcf4eed041746b0ac Mon Sep 17 00:00:00 2001 From: Gabriel Brammer Date: Wed, 25 Jan 2023 14:58:00 +0100 Subject: [PATCH] use sedpy for smoothing --- eazy/__init__.py | 18 +++++----- eazy/photoz.py | 4 ++- eazy/sps.py | 44 ++++++++++++++++++----- eazy/templates.py | 67 ++++++++++++++++++++++-------------- eazy/tests/test_templates.py | 4 +-- requirements.txt | 2 +- setup.cfg | 3 +- 7 files changed, 92 insertions(+), 50 deletions(-) diff --git a/eazy/__init__.py b/eazy/__init__.py index 304190c..4228ec1 100644 --- a/eazy/__init__.py +++ b/eazy/__init__.py @@ -23,15 +23,15 @@ 'git+https://github.com/gbrammer/dust_extinction.git') -# Hot fix for importing prospector without SPS_HOME variable set -try: - from prospect.utils.smoothing import smoothspec -except (FileNotFoundError, TypeError): - if 'SPS_HOME' not in os.environ: - sps_home = 'xxxxdummyxxxx' #os.path.dirname(__file__) - print(f'Warning: setting environment variable SPS_HOME={sps_home} ' - 'to be able to import prospect.') - os.environ['SPS_HOME'] = sps_home +## Hot fix for importing prospector without SPS_HOME variable set +# try: +# from prospect.utils.smoothing import smoothspec +# except (FileNotFoundError, TypeError): +# if 'SPS_HOME' not in os.environ: +# sps_home = 'xxxxdummyxxxx' #os.path.dirname(__file__) +# print(f'Warning: setting environment variable SPS_HOME={sps_home} ' +# 'to be able to import prospect.') +# os.environ['SPS_HOME'] = sps_home def fetch_eazy_photoz(): diff --git a/eazy/photoz.py b/eazy/photoz.py index 14c65b4..65b4939 100644 --- a/eazy/photoz.py +++ b/eazy/photoz.py @@ -4539,7 +4539,9 @@ def sps_parameters(self, UBVJ=DEFAULT_UBVJ_FILTERS, extra_rf_filters=DEFAULT_RF_ pad_width=rf_pad_width, max_err=rf_max_err, percentiles=[16,50,84], - verbose=False, simple=simple) + verbose=False, + n_proc=n_proc, + simple=simple) extra_tempfilt, extra_lc, extra_rest = _ex diff --git a/eazy/sps.py b/eazy/sps.py index 343301c..12f03ad 100644 --- a/eazy/sps.py +++ b/eazy/sps.py @@ -553,6 +553,10 @@ class ExtendedFsps(StellarPopulation): line_av_func = None + # Smoothing parameters for fitting + lsf_func = None + FFT_SMOOTH = False + #_meta_bands = ['v'] @property @@ -1389,7 +1393,9 @@ def fit_spec(self, wave_obs, flux_obs, err_obs, mask=None, plist=['tage', 'Av', kwargs[p] = theta0[i] # Model test - margs = (self, plist, wave_obs, flux_obs, err_obs, mask, bspl, kwargs, 'model') + margs = (self, plist, wave_obs, flux_obs, err_obs, + mask, bspl, kwargs, 'model') + flux_model, Anorm, chi2_init = self.objfun_fitspec(theta0, *margs) if show: @@ -1398,7 +1404,14 @@ def fit_spec(self, wave_obs, flux_obs, err_obs, mask=None, plist=['tage', 'Av', plt.close('all') fig = plt.figure(figsize=(12, 6)) - plt.errorbar(wave_obs[mask], flux_obs[mask], err_obs[mask], color='k', alpha=0.5, linestyle='None', marker='.') + plt.errorbar(wave_obs[mask], + flux_obs[mask], + err_obs[mask], + color='k', + alpha=0.5, + linestyle='None', + marker='.') + plt.plot(wave_obs, flux_model, color='pink', linewidth=2, alpha=0.8) else: fig = None @@ -1407,8 +1420,14 @@ def fit_spec(self, wave_obs, flux_obs, err_obs, mask=None, plist=['tage', 'Av', #lsq_kwargs['diff_step'] = np.array(steps)/2. #lsq_kwargs['diff_step'] = 0.05 lsq_kwargs['diff_step'] = steps - lmargs = (self, plist, wave_obs, flux_obs, err_obs, mask, bspl, kwargs, 'least_squares verbose') - _res = least_squares(self.objfun_fitspec, theta0, bounds=bounds, args=lmargs, **lsq_kwargs) + lmargs = (self, plist, wave_obs, flux_obs, err_obs, + mask, bspl, kwargs, 'least_squares verbose') + + _res = least_squares(self.objfun_fitspec, + theta0,bounds=bounds, + args=lmargs, + **lsq_kwargs, + ) fit_model, Anorm, chi2_fit = self.objfun_fitspec(_res.x, *margs) @@ -1441,11 +1460,18 @@ def objfun_fitspec(theta, self, plist, wave_obs, flux_obs, err_obs, mask, bspl, kwargs[p] = theta[i] templ = self.get_full_spectrum(**kwargs) - flux_model = templ.resample(wave_obs, z=self.params['zred'], - in_place=False, - return_array=True, interp_func=interp_func) - - flux_model = flux_model.flatten() + + flux_model = templ.to_observed_frame(z=self.params['zred'], + lsf_func=self.lsf_func, + clip_wavelengths=None, + wavelengths=wave_obs, + smoothspec_kwargs={'fftsmooth':self.FFT_SMOOTH}, + ) + # flux_model = templ.resample(wave_obs, z=self.params['zred'], + # in_place=False, + # return_array=True, interp_func=interp_func) + + flux_model = np.squeeze(flux_model.flux_flam()) if mask is None: mask = np.isfinite(flux_model+flux_obs+err_obs) & (err_obs > 0) diff --git a/eazy/templates.py b/eazy/templates.py index b2135d4..8ef5a43 100644 --- a/eazy/templates.py +++ b/eazy/templates.py @@ -688,9 +688,9 @@ def set_fnu(self): pass - def smooth_velocity(self, velocity_smooth, in_place=True, raise_error=False): + def smooth_velocity(self, velocity_smooth, in_place=True, raise_error=False, smoothspec_kwargs={'fftsmooth':True}): """ - Smooth template in velocity using ``astro-prospector`` + Smooth template in velocity using ``sedpy`` (formerly in ``astro-prospector``) Parameters ---------- @@ -702,15 +702,15 @@ def smooth_velocity(self, velocity_smooth, in_place=True, raise_error=False): return a new `~eazy.templates.Template` object. raise_error : bool - If ``from prospect.utils.smoothing import smooth_vel`` fails, + If ``from sedpy.smoothing import smooth_vel`` fails, raise an exception or die quietly. """ try: - from prospect.utils.smoothing import smooth_vel + from sedpy.smoothing import smoothspec except: if raise_error: - raise ImportError("Couldn't import `prospect.utils.smoothing") + raise ImportError("Couldn't import `sedpy.smoothing") else: return None @@ -720,8 +720,13 @@ def smooth_velocity(self, velocity_smooth, in_place=True, raise_error=False): else: return self - sm_flux = np.array([smooth_vel(self.wave, self.flux[i,:], self.wave, - velocity_smooth) + sm_flux = np.array([smoothspec(self.wave, + self.flux[i,:], + outwave=self.wave, + resolution=velocity_smooth, + smoothtype='vel', + **smoothspec_kwargs + ) for i in range(self.NZ)]) sm_flux[~np.isfinite(sm_flux)] = 0. @@ -737,13 +742,13 @@ def smooth_velocity(self, velocity_smooth, in_place=True, raise_error=False): redshifts=self.redshifts) - def to_observed_frame(self, z=0, scalar=1., extra_sigma=0, lsf_func=None, to_air=False, wavelengths=None, smoothspec_kwargs={'fftsmooth':False}, include_igm=True, clip_wavelengths=[4500,9400]): + def to_observed_frame(self, z=0, scalar=1., extra_sigma=0, lsf_func=None, to_air=False, wavelengths=None, smoothspec_kwargs={'fftsmooth':False}, include_igm=True, clip_wavelengths=[4500,9400], as_template=True): """ Smooth and resample to observed-frame wavelengths, including an optional Line Spread Function (LSF) Note that the smoothing is performed with - `prospect.utils.smoothing.smoothspec `_, + `sedpy.smoothing.smoothspec `_, which doesn't integrate precisely over "pixels" for spectral resolutions that are similar to or less than the target smoothing factor. @@ -780,8 +785,8 @@ def to_observed_frame(self, z=0, scalar=1., extra_sigma=0, lsf_func=None, to_air (e.g., MUSE) spectrum smoothspec_kwargs : dict - Extra keyword arguments to pass to the Prospector smoothing - function `prospect.utils.smoothing.smoothspec `_. + Extra keyword arguments to pass to the `sedpy`/`prospector` smoothing + function `sedpy.smoothing.smoothspec `_. When testing with very high resolution templates around a specific wavelength, ``smoothspec_kwargs = {'fftsmooth':True}`` did not always work as expected, so be careful with this option (which is @@ -793,15 +798,16 @@ def to_observed_frame(self, z=0, scalar=1., extra_sigma=0, lsf_func=None, to_air clip_wavelengths : [float, float] Trim the full observed-frame wavelength array before convolving. The defaults bracket the nominal MUSE range. - + Returns ------- - tobs : `~eazy.template.Template` - Smoothed and resampled `~eazy.template.Template` object + tobs : `~eazy.template.Template` or array-like + Smoothed and resampled `~eazy.template.Template` object or flux array, + depending on ``is_array`` """ from astropy.stats import gaussian_sigma_to_fwhm - from prospect.utils.smoothing import smoothspec + from sedpy.smoothing import smoothspec wobs = self.wave*(1+z) @@ -829,7 +835,10 @@ def to_observed_frame(self, z=0, scalar=1., extra_sigma=0, lsf_func=None, to_air 'Angstroms') else: clip = wobs > 0 - + if wavelengths is not None: + clip &= wobs > 0.95*wavelengths.min() + clip &= wobs < 1.05*wavelengths.max() + if lsf_func in ['Bacon']: # UDF-10 LSF from Bacon et al. 2017 bacon_lsf_fwhm = lambda w: 5.866e-8 * w**2 - 9.187e-4*w + 6.04 @@ -841,7 +850,7 @@ def to_observed_frame(self, z=0, scalar=1., extra_sigma=0, lsf_func=None, to_air elif hasattr(lsf_func, '__call__'): lsf_sigma = lsf_func(wobs[clip])/wobs[clip]*3.e5 lsf_func_name = 'user' - + else: lsf_sigma = 0. lsf_func_name = None @@ -852,22 +861,30 @@ def to_observed_frame(self, z=0, scalar=1., extra_sigma=0, lsf_func=None, to_air # In Angstroms smooth_lambda = vel_sigma / 3.e5 * wobs[clip] - # Do the smoothing - flux_smooth = smoothspec(wobs[clip], - (self.flux_flam(z=z)*igmz*scalar)[clip], - resolution=smooth_lambda, - smoothtype='lsf', **smoothspec_kwargs) + _templ_flux = self.flux_flam(z=z) + if np.allclose(smooth_lambda, 0.): + flux_smooth = (_templ_flux*igmz*scalar)[clip] + else: + flux_smooth = smoothspec(wobs[clip], + (_templ_flux*igmz*scalar)[clip], + resolution=smooth_lambda, + smoothtype='lsf', + **smoothspec_kwargs, + ) + newname = self.name + f' z={z:.3f}' if lsf_func_name is not None: newname += ' + ' + lsf_func_name if extra_sigma > 0: newname += ' + {extra_sigma:.1f} km/s' - + tobs = Template(arrays=(wobs[clip], flux_smooth), - name=newname, resample_wave=wavelengths, - redshifts=[z]) + name=newname, + resample_wave=wavelengths, + redshifts=[z], + ) return tobs diff --git a/eazy/tests/test_templates.py b/eazy/tests/test_templates.py index ee13ed3..4e92a3e 100644 --- a/eazy/tests/test_templates.py +++ b/eazy/tests/test_templates.py @@ -307,9 +307,7 @@ def test_template_smoothing(): """ from astropy.stats import gaussian_sigma_to_fwhm - - from prospect.utils.smoothing import smooth_vel - + #### Template with delta function line xtest = np.linspace(6550, 6576, 1024) ytest = xtest*0; ytest[len(xtest)//2] = 1 diff --git a/requirements.txt b/requirements.txt index 6401895..64d0b75 100644 --- a/requirements.txt +++ b/requirements.txt @@ -17,7 +17,7 @@ git+https://github.com/gbrammer/dust_extinction.git git+https://github.com/gbrammer/dust_attenuation.git # Needed for template smoothing, doesn't bring whole fsps distro -# astro-sedpy +astro-sedpy>=0.3 # astro-prospector diff --git a/setup.cfg b/setup.cfg index 704386c..69cb0cc 100644 --- a/setup.cfg +++ b/setup.cfg @@ -27,8 +27,7 @@ install_requires = peakutils tqdm h5py - astro-sedpy - astro-prospector + astro-sedpy>=0.3 packages = find: include_package_data = True