Merge pull request #137 from desihub/snr-narrowlines

require more pixels to estimate the variance in the emission-line amplitude

doc/changes.rst

@@ -11,13 +11,25 @@ Change Log
 * Do not constrain the SPS age by default [`PR #132`_].
 * Bug fix of emission-line subtracted Dn(4000) measurement [`PR #135`_].
 * Update IGM attenuation coefficients [`PR #136`_].
+* Several significant changes [`PR #137`_]:
+
+  * Record the observed-space emission-line amplitude in ``_AMP`` and move the
+    model-space amplitude to ``_MODELAMP``. 
+  * Demand at least 12 pixels to measure the scatter in the pixels under the
+    line (therefore ``_AMP_IVAR`` should be more reliable for narrow lines).
+  * Major bug fix whereby the model emission-line spectra were not being
+    convolved with the resolution matrix.
+  * Redefine ``_CHI2`` for an emission line as the observed not reduced chi2.
+  * Switch from (deprecated) ``pkg_resources`` to ``importlib``.
+  * Updated documentation (data model) and several non-negligible speed-ups.
 
 .. _`PR #115`: https://github.com/desihub/fastspecfit/pull/115
 .. _`PR #116`: https://github.com/desihub/fastspecfit/pull/116
 .. _`PR #120`: https://github.com/desihub/fastspecfit/pull/120
 .. _`PR #132`: https://github.com/desihub/fastspecfit/pull/132
 .. _`PR #135`: https://github.com/desihub/fastspecfit/pull/135
 .. _`PR #136`: https://github.com/desihub/fastspecfit/pull/136
+.. _`PR #137`: https://github.com/desihub/fastspecfit/pull/137
 
 2.1.2 (2023-04-01)
 ------------------

py/fastspecfit/continuum.py

@@ -14,61 +14,6 @@
 from fastspecfit.io import FLUXNORM
 from fastspecfit.util import C_LIGHT
 
-#import numba
-#@numba.jit(nopython=True)
-def nnls(A, b, maxiter=None, eps=1e-7, v=False):
-    # https://en.wikipedia.org/wiki/Non-negative_least_squares#Algorithms
-    # not sure what eps should be set to
-    m, n = A.shape
-    if b.shape[0] != m:
-        raise ValueError()#f"Shape mismatch: {A.shape} {b.shape}. Expected: (m, n) (m) ")
-    if maxiter is None:
-        # this is what scipy does when maxiter is not specified
-        maxiter = 3 * n
-    # init
-    P = np.zeros(n).astype(bool)
-    R = np.ones(n).astype(bool)
-    x = np.zeros(n)
-    s = np.zeros(n)
-    # compute these once ahead of time
-    ATA = A.T.dot(A)
-    ATb = A.T.dot(b)
-    # main loop
-    w = ATb - ATA.dot(x)
-    j = np.argmax(R * w)
-    c = 0 # iteration count
-    # while R != {} and max(w[R]) > eps
-    while np.any(R) and w[j] > eps:
-        #if v: print(f"{c=}", f"{P=}\n {R=}\n {w=}\n {j=}")
-        # add j to P, remove j from R
-        P[j], R[j] = True, False
-        s[P] = np.linalg.inv(ATA[P][:, P]).dot(ATb[P])
-        s[R] = 0
-        d = 0 # inner loop iteration count, for debugging
-        #if v: print(f"{c=}", f"{P=}\n {R=}\n {s=}")
-        # while P != {} and min(s[P]) < eps
-        # make sure P is not empty before checking min s[P]
-        while np.any(P) and np.min(s[P]) < eps:
-            i = P & (s < eps)
-            #if v: print(f" {d=}", f"  {P=}\n  {i=}\n  {s=}\n  {x=}")
-            a = np.min(x[i] / (x[i] - s[i]))
-            x = x + a * (s - x)
-            j = P & (x < eps)
-            R[j], P[j] = True, False
-            s[P] = np.linalg.inv(ATA[P][:, P]).dot(ATb[P])
-            s[R] = 0
-            d += 1
-        x[:] = s
-        # w = A.T.dot(b - A.dot(x))
-        w = ATb - ATA.dot(x)
-        j = np.argmax(R * w)
-        #if v: print(f"{c=}", f"{P=}\n {R=}\n {w=}\n {j=}")
-        c += 1
-        if c >= maxiter:
-            break
-    res = np.linalg.norm(A.dot(x) - b)
-    return x, res
-
 def _smooth_continuum(wave, flux, ivar, redshift, medbin=150, 
                       smooth_window=50, smooth_step=10, maskkms_uv=3000.0, 
                       maskkms_balmer=1000.0, maskkms_narrow=200.0,
@@ -1485,7 +1430,7 @@ def templates2data(self, _templateflux, _templatewave, redshift=0.0, dluminosity
 
     @staticmethod
     def call_nnls(modelflux, flux, ivar, xparam=None, debug=False,
-                  interpolate_coeff=False, xlabel=None, log=None):
+                  interpolate_coeff=False, xlabel=None, png=None, log=None):
         """Wrapper on nnls.
     
         Works with both spectroscopic and photometric input and with both 2D and
@@ -1535,11 +1480,14 @@ def call_nnls(modelflux, flux, ivar, xparam=None, debug=False,
 
         try:
             imin = find_minima(chi2grid)[0]
-            xbest, xerr, chi2min, warn = minfit(xparam[imin-1:imin+2], chi2grid[imin-1:imin+2])
+            if debug:
+                xbest, xerr, chi2min, warn, (a, b, c) = minfit(xparam[imin-1:imin+2], chi2grid[imin-1:imin+2], return_coeff=True)
+            else:
+                xbest, xerr, chi2min, warn = minfit(xparam[imin-1:imin+2], chi2grid[imin-1:imin+2])
         except:
             errmsg = 'A problem was encountered minimizing chi2.'
             log.warning(errmsg)
-            imin, xbest, xerr, chi2min, warn = 0, 0.0, 0.0, 0.0, 1
+            imin, xbest, xerr, chi2min, warn, (a, b, c) = 0, 0.0, 0.0, 0.0, 1, (0., 0., 0.)
 
         if warn == 0:
             xivar = 1.0 / xerr**2
@@ -1565,29 +1513,37 @@ def call_nnls(modelflux, flux, ivar, xparam=None, debug=False,
 
         if debug:
             if xivar > 0:
-                leg = r'${:.3f}\pm{:.3f}\ (\chi^2_{{min}}={:.2f})$'.format(xbest, 1/np.sqrt(xivar), chi2min)
+                leg = r'${:.1f}\pm{:.1f}$'.format(xbest, 1 / np.sqrt(xivar))
+                #leg = r'${:.3f}\pm{:.3f}\ (\chi^2_{{min}}={:.2f})$'.format(xbest, 1/np.sqrt(xivar), chi2min)
             else:
                 leg = r'${:.3f}$'.format(xbest)
 
+            if xlabel:
+                leg = f'{xlabel} = {leg}'
+
             import matplotlib.pyplot as plt
+            import seaborn as sns
+            sns.set(context='talk', style='ticks', font_scale=0.8)
             fig, ax = plt.subplots()
-            ax.scatter(xparam, chi2grid)
-            ax.scatter(xparam[imin-1:imin+2], chi2grid[imin-1:imin+2], color='red')
-            #ax.set_ylim(chi2min*0.95, np.max(chi2grid[imin-1:imin+2])*1.05)
-            #ax.plot(xx, np.polyval([aa, bb, cc], xx), ls='--')
-            ax.axvline(x=xbest, color='k')
-            if xivar > 0:
-                ax.axhline(y=chi2min, color='k')
-            #ax.set_yscale('log')
-            #ax.set_ylim(chi2min, 63.3)
+            ax.scatter(xparam, chi2grid-chi2min, marker='s', s=50, color='gray', edgecolor='k')
+            #ax.scatter(xparam[imin-1:imin+2], chi2grid[imin-1:imin+2]-chi2min, color='red')
+            if not np.all(np.array([a, b, c]) == 0):
+                ax.plot(xparam, np.polyval([a, b, c], xparam)-chi2min, lw=2, ls='--')
+            #ax.axvline(x=xbest, color='k')
+            #if xivar > 0:
+            #    ax.axhline(y=chi2min, color='k')
             if xlabel:
                 ax.set_xlabel(xlabel)
                 #ax.text(0.03, 0.9, '{}={}'.format(xlabel, leg), ha='left',
                 #        va='center', transform=ax.transAxes)
-            ax.text(0.03, 0.9, leg, ha='left', va='center', transform=ax.transAxes)
-            ax.set_ylabel(r'$\chi^2$')
+            ax.text(0.9, 0.9, leg, ha='right', va='center', transform=ax.transAxes)
+            ax.set_ylabel(r'$\Delta\chi^2$')
+            #ax.set_ylabel(r'$\chi^2 - {:.1f}$'.format(chi2min))
             #fig.savefig('qa-chi2min.png')
-            fig.savefig('desi-users/ioannis/tmp/qa-chi2min.png')
+            fig.tight_layout()
+            if png:
+                log.info(f'Writing {png}')
+                fig.savefig(png)
 
         return chi2min, xbest, xivar, bestcoeff
 
@@ -1797,7 +1753,7 @@ def _kcorr_and_absmag(filters_out, band_shift):
         return kcorr, absmag, ivarabsmag, bestmaggies, lums, cfluxes
 
 def continuum_specfit(data, result, templatecache, nophoto=False, constrain_age=False,
-                      fastphot=False, log=None, verbose=False):
+                      fastphot=False, log=None, verbose=False, debug_plots=False):
     """Fit the non-negative stellar continuum of a single spectrum.
 
     Parameters
@@ -1960,9 +1916,10 @@ def _younger_than_universe(age, tuniv, agepad=0.5):
             vdispchi2min, vdispbest, vdispivar, _ = CTools.call_nnls(
                 ztemplateflux_vdisp[Ivdisp, :, :], 
                 specflux[Ivdisp], specivar[Ivdisp],
-                xparam=templatecache['vdisp'], xlabel=r'$\sigma$ (km/s)', debug=False, log=log)
+                xparam=templatecache['vdisp'], xlabel=r'$\sigma$ (km/s)', log=log,
+                debug=debug_plots, png='deltachi2-vdisp.png')
             log.info('Fitting for the velocity dispersion took {:.2f} seconds.'.format(time.time()-t0))
-
+            
             if vdispivar > 0:
                 # Require vdisp to be measured with S/N>1, which protects
                 # against tiny ivar becomming infinite in the output table.