fix miscellaneous bugs

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)
 ------------------

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)

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)

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

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:

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()

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: