diff --git a/doc/changes.rst b/doc/changes.rst index fceba9dc..fc6c8708 100644 --- a/doc/changes.rst +++ b/doc/changes.rst @@ -5,6 +5,7 @@ Change Log 3.0.0 (not released yet) ------------------------ +* Miscellaneous bug fixes [`PR #185`_]. * Use Numba-based Resolution implementation with faster multiply. Also tweak parameters of continuum least_squares optimizer for speed. [`PR #181`_] @@ -18,6 +19,7 @@ Change Log .. _`PR #179`: https://github.com/desihub/fastspecfit/pull/179 .. _`PR #180`: https://github.com/desihub/fastspecfit/pull/180 .. _`PR #181`: https://github.com/desihub/fastspecfit/pull/181 +.. _`PR #185`: https://github.com/desihub/fastspecfit/pull/185 2.5.2 (2024-04-28) ------------------ diff --git a/py/fastspecfit/continuum.py b/py/fastspecfit/continuum.py index f42742f5..f8de58ad 100644 --- a/py/fastspecfit/continuum.py +++ b/py/fastspecfit/continuum.py @@ -1375,11 +1375,11 @@ def continuum_fastspec(redshift, objflam, objflamivar, CTools, log.critical(errmsg) raise ValueError(errmsg) - ncam = len(data['snr']) + ncam = len(data['cameras']) if ncam == 1: - snrmsg = f"Median S/N_{data['cameras']}={data['snr'][0]:.2f}" + snrmsg = f"Median spectral S/N_{data['cameras']}={data['snr'][0]:.2f}" else: - snrmsg = f"Median S/N_{data['cameras'][0]}={data['snr'][0]:.2f}" + snrmsg = f"Median spectral S/N_{data['cameras'][0]}={data['snr'][0]:.2f}" for icam in np.arange(ncam-1)+1: snrmsg += f" S/N_{data['cameras'][icam]}={data['snr'][icam]:.2f}" log.info(snrmsg) diff --git a/py/fastspecfit/emlines.py b/py/fastspecfit/emlines.py index e5bafe8e..95dcc959 100644 --- a/py/fastspecfit/emlines.py +++ b/py/fastspecfit/emlines.py @@ -909,6 +909,8 @@ def get_boundaries(A, v_lo, v_hi): else: stats = narrow_stats stats.append((linesigma, linez)) + else: + flux, flux_ivar = 0., 0. else: flux, flux_ivar = 0., 0. @@ -937,10 +939,13 @@ def get_boundaries(A, v_lo, v_hi): ew, ewivar = 0., 0. if cmed != 0. and civar != 0.: - if flux > 0. and flux_ivar > 0.: - # add the uncertainties in flux and the continuum in quadrature - ew = flux / cmed / (1. + redshift) # rest frame [A] - ewivar = (1. + redshift)**2 / (1. / (cmed**2 * flux_ivar) + flux**2 / (cmed**4 * civar)) + try: + if flux > 0. and flux_ivar > 0.: + # add the uncertainties in flux and the continuum in quadrature + ew = flux / cmed / (1. + redshift) # rest frame [A] + ewivar = (1. + redshift)**2 / (1. / (cmed**2 * flux_ivar) + flux**2 / (cmed**4 * civar)) + except: + import pdb ; pdb.set_trace() # upper limit on the flux is defined by snrcut*cont_err*sqrt(2*pi)*linesigma fluxlimit = np.sqrt(2. * np.pi) * linesigma_ang / np.sqrt(civar) # * u.erg/(u.second*u.cm**2) diff --git a/py/fastspecfit/io.py b/py/fastspecfit/io.py index ab7b3466..5fc15d06 100644 --- a/py/fastspecfit/io.py +++ b/py/fastspecfit/io.py @@ -11,7 +11,7 @@ import fitsio from astropy.table import Table -from fastspecfit.logger import log +from fastspecfit.logger import log, DEBUG from fastspecfit.singlecopy import sc_data from fastspecfit.photometry import Photometry from fastspecfit.util import FLUXNORM, ZWarningMask @@ -903,7 +903,7 @@ def one_spectrum(iobj, specdata, meta, ebv, fastphot, for icam, camera in enumerate(specdata['cameras']): # Check whether the camera is fully masked. if np.sum(specdata['ivar0'][icam]) == 0: - log.warning(f'Dropping fully masked camera {camera}.') + log.warning(f'Dropping fully masked camera {camera} [specdata["uniqueid"]].') else: ivar = specdata['ivar0'][icam] mask = specdata['mask0'][icam] @@ -917,7 +917,7 @@ def one_spectrum(iobj, specdata, meta, ebv, fastphot, ivar[mask != 0] = 0. if np.all(ivar == 0.): - log.warning(f'Dropping fully masked camera {camera}.') + log.warning(f'Dropping fully masked camera {camera} [specdata["uniqueid"]].') else: # interpolate over pixels where the resolution matrix is masked I = (mask != 0) @@ -1313,7 +1313,7 @@ def one_stacked_spectrum(iobj, specdata, synthphot): for icam, camera in enumerate(specdata['cameras']): # Check whether the camera is fully masked. if np.sum(specdata['ivar0'][icam]) == 0: - log.warning(f'Dropping fully masked camera {camera}.') + log.warning(f'Dropping fully masked camera {camera} [specdata["uniqueid"]].') else: ivar = specdata['ivar0'][icam] mask = specdata['mask0'][icam] @@ -1327,7 +1327,7 @@ def one_stacked_spectrum(iobj, specdata, synthphot): ivar[mask != 0] = 0 if np.all(ivar == 0): - log.warning(f'Dropping fully masked camera {camera}.') + log.warning(f'Dropping fully masked camera {camera} [specdata["uniqueid"]].') else: cameras.append(camera) npixpercamera.append(len(specdata['wave0'][icam])) # number of pixels in this camera diff --git a/py/fastspecfit/linemasker.py b/py/fastspecfit/linemasker.py index d4ebe5d9..c554c4d5 100644 --- a/py/fastspecfit/linemasker.py +++ b/py/fastspecfit/linemasker.py @@ -76,7 +76,7 @@ def _get_contpix(zlinewaves, sigmas, nsigma_factor=2., linemask=None, lya=False) pix = {'linepix': {}, 'contpix': {}} if patchMap is not None: - pix.update({'patch_contpix': {}, 'dropped': [], 'merged': []}) + pix.update({'patch_contpix': {}, 'dropped': [], 'merged': [], 'merged_from': []}) patchids = list(patchMap.keys()) npatch = len(patchids) @@ -177,6 +177,7 @@ def _get_contpix(zlinewaves, sigmas, nsigma_factor=2., linemask=None, lya=False) del pix['patch_contpix'][patchids[ipatch+1]] pix['patch_contpix'][newpatchid] = newcontpix pix['merged'].append(newpatchid) + pix['merged_from'].append([patchid, patchids[ipatch+1]]) if patchMap is not None: if len(pix['dropped']) > 0: @@ -298,8 +299,8 @@ def _fit_patches(continuum_patches, patchMap, linemodel, debug_plots=False, # In the case that patches have been merged, update the # continuum_patches table and patchMap dictionary. if len(pix['merged']) > 0: - for newpatchid in np.atleast_1d(pix['merged']): - oldpatchids = list(newpatchid) + for oldpatchids, newpatchid in zip(np.atleast_1d(pix['merged_from']), + np.atleast_1d(pix['merged'])): O = np.where(np.isin(continuum_patches['patchid'].value, oldpatchids))[0] # update continuum_patches @@ -316,6 +317,7 @@ def _fit_patches(continuum_patches, patchMap, linemodel, debug_plots=False, newI.append(patchMap[oldpatchid][1]) newJ.append(patchMap[oldpatchid][2]) del patchMap[oldpatchid] + patchMap[newpatchid] = (np.hstack(newlines), np.hstack(newI), np.hstack(newJ)) @@ -549,7 +551,7 @@ def _niceline(line): xx.set_ylim(ymin, ymax) xx.legend(loc='upper left', fontsize=8, ncols=nlegcol) xx.set_title(f'Patch {patchid}') - for rem in range(ipatch+1, ncols*nrowspatch): + for rem in range(ipatch+1, ncols*nrows): ax.flat[rem].axis('off') if ax.ndim == 1: diff --git a/py/fastspecfit/logger.py b/py/fastspecfit/logger.py index 62add693..56975b2b 100644 --- a/py/fastspecfit/logger.py +++ b/py/fastspecfit/logger.py @@ -10,6 +10,6 @@ This needs to be in its own file to prevent circular imports with other code. """ -from desiutil.log import get_logger +from desiutil.log import get_logger, DEBUG log = get_logger() diff --git a/py/fastspecfit/photometry.py b/py/fastspecfit/photometry.py index 27ea2edc..75440fa5 100644 --- a/py/fastspecfit/photometry.py +++ b/py/fastspecfit/photometry.py @@ -10,7 +10,7 @@ import fitsio from astropy.table import Table -from fastspecfit.logger import log +from fastspecfit.logger import log, DEBUG from fastspecfit.util import trapz, C_LIGHT, FLUXNORM @@ -213,7 +213,7 @@ def kcorr_and_absmag(self, nanomaggies, ivar_nanomaggies, redshift, dmod, synth_maggies_in = np.zeros(len(maggies)) return kcorr, absmag, ivarabsmag, synth_absmag, synth_maggies_in - maggies = nanomaggies * 1e-9 + maggies = nanomaggies * 1e-9 ivarmaggies = (ivar_nanomaggies / 1e-9**2) * self.bands_to_fit filters_in = self.filters[photsys] @@ -244,16 +244,16 @@ def kcorr_and_absmag(self, nanomaggies, ivar_nanomaggies, redshift, dmod, # m_R = M_Q + DM(z) + K_QR(z) or # M_Q = m_R - DM(z) - K_QR(z) - absmag = -2.5 * np.log10(maggies[oband]) - dmod - kcorr - - C = 0.8483036976765437 # (0.4 * np.log(10.))**2 - ivarabsmag = maggies[oband]**2 * ivarmaggies[oband] * C + absmag = np.copy(synth_absmag) + ivarabsmag = np.zeros_like(absmag) # if we use synthesized photometry then ivarabsmag is zero # (which should never happen?) - I = (maggies[oband] * np.sqrt(ivarmaggies[oband]) <= snrmin) - absmag[I] = synth_absmag[I] - ivarabsmag[I] = 0. + I = (maggies[oband] * np.sqrt(ivarmaggies[oband]) > snrmin) + if np.any(I): + C = 0.8483036976765437 # (0.4 * np.log(10.))**2 + absmag[I] = -2.5 * np.log10(maggies[oband[I]]) - dmod - kcorr[I] + ivarabsmag[I] = maggies[oband[I]]**2 * ivarmaggies[oband[I]] * C return kcorr, absmag, ivarabsmag, synth_maggies_in @@ -334,7 +334,7 @@ def get_ab_maggies(filters, flux, wave): maggies = Photometry.get_ab_maggies_unchecked(filters, flux, wave) except: # pad in case of an object at very high redshift (z > 5.5) - log.warning('Padding model spectrum due to insufficient wavelength coverage to synthesize photometry.') + log.debug('Padding input spectrum due to insufficient wavelength coverage to synthesize photometry.') padflux, padwave = filters.pad_spectrum(flux, wave, axis=0, method='edge') if flux.ndim > 1: