Skip to content

Commit

Permalink
Merge pull request #27 from gbrammer/sedpy-smoothing
Browse files Browse the repository at this point in the history 
Use smoothing functions from `sedpy`
  • Loading branch information
@gbrammer
gbrammer committed Jan 25, 2023
2 parents dce4fc4 + d1adef8 commit 1ea4978
Show file tree
Hide file tree
Showing 7 changed files with 92 additions and 50 deletions.
18 changes: 9 additions & 9 deletions eazy/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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():
Expand Down
4 changes: 3 additions & 1 deletion eazy/photoz.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
44 changes: 35 additions & 9 deletions eazy/sps.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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:
Expand All @@ -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
Expand All @@ -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)

Expand Down Expand Up @@ -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)
Expand Down
67 changes: 42 additions & 25 deletions eazy/templates.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
----------
Expand All @@ -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

Expand All @@ -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.
Expand All @@ -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 <https://prospect.readthedocs.io/en/latest/api/utils_api.html>`_,
`sedpy.smoothing.smoothspec <https://prospect.readthedocs.io/en/latest/api/utils_api.html>`_,
which doesn't integrate precisely over "pixels" for spectral
resolutions that are similar to or less than the target smoothing
factor.
Expand Down Expand Up @@ -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 <https://prospect.readthedocs.io/en/latest/api/utils_api.html>`_.
Extra keyword arguments to pass to the `sedpy`/`prospector` smoothing
function `sedpy.smoothing.smoothspec <https://prospect.readthedocs.io/en/latest/api/utils_api.html>`_.
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
Expand All @@ -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)

Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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

Expand Down
4 changes: 1 addition & 3 deletions eazy/tests/test_templates.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion requirements.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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


3 changes: 1 addition & 2 deletions setup.cfg
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,8 +27,7 @@ install_requires =
peakutils
tqdm
h5py
astro-sedpy
astro-prospector
astro-sedpy>=0.3
packages = find:
include_package_data = True

Expand Down

0 comments on commit 1ea4978

Please sign in to comment.