add PSF model error; chi2pix doesn't use masked pixels

doc/release_notes.rst

@@ -11,6 +11,7 @@ specter Release Notes
 * Added a documentation page for the specter API.
 * Added full_output option to ex2d to get model image and metrics based upon
   goodness of fit
+* PSFs can specify their model error with PSFERR header keyword; default 0.01
 
 0.4.1 (2016-03-10)
 ------------------

py/specter/extract/ex2d.py

@@ -140,7 +140,7 @@ def ex2d(image, imageivar, psf, specmin, nspec, wavelengths, xyrange=None,
             ivar[iispec, iwave:iwave+wavesize+1] = specivar[:, nlo:-nhi]
 
             if full_output:
-                A = results['A']
+                A = results['A'].copy()
                 xflux = results['xflux']
 
                 #- number of spectra and wavelengths for this sub-extraction
@@ -153,9 +153,23 @@ def ex2d(image, imageivar, psf, specmin, nspec, wavelengths, xyrange=None,
                 #- Fraction of input pixels that are unmasked for each flux bin
                 subpixmask_fraction = 1.0-(A.T.dot(subivar.ravel()>0)).reshape(subnspec, subnwave)
 
-                #- Weighted chi2 of pixels that contribute to each flux bin
-                chi = (subimg - submodel) * np.sqrt(subivar)
-                chi2x = (A.T.dot(chi.ravel()**2) / A.sum(axis=0)).reshape(subnspec, subnwave)
+                #- original weighted chi2 of pixels that contribute to each flux bin
+                # chi = (subimg - submodel) * np.sqrt(subivar)
+                # chi2x = (A.T.dot(chi.ravel()**2) / A.sum(axis=0)).reshape(subnspec, subnwave)
+
+                #- pixel variance including input noise and PSF model errors
+                modelivar = (submodel*psf.psferr + 1e-32)**-2
+                ii = (modelivar > 0) & (subivar > 0)
+                totpix_ivar = np.zeros(submodel.shape)
+                totpix_ivar[ii] = 1.0 / (1.0/modelivar[ii] + 1.0/subivar[ii])
+
+                #- Weighted chi2 of pixels that contribute to each flux bin;
+                #- only use unmasked pixels and avoid dividing by 0
+                chi = (subimg - submodel) * np.sqrt(totpix_ivar)
+                psfweight = A.T.dot(totpix_ivar.ravel()>0)
+                bad = (psfweight == 0.0)
+                chi2x = (A.T.dot(chi.ravel()**2) * ~bad) / (psfweight + bad)
+                chi2x = chi2x.reshape(subnspec, subnwave)
 
                 #- outputs
                 #- TODO: watch out for edge effects on overlapping regions of submodels

py/specter/psf/gausshermite.py

@@ -48,7 +48,13 @@ def __init__(self, filename):
         #- Other necessary keywords
         self.npix_x = hdr['NPIX_X']
         self.npix_y = hdr['NPIX_Y']
-
+
+        #- PSF model error
+        if 'PSFERR' in hdr:
+            self.psferr = hdr['PSFERR']
+        else:
+            self.psferr = 0.01
+
         #- Load the parameters into self.coeff dictionary keyed by PARAM
         #- with values as TraceSets for evaluating the Legendre coefficients
         data = fx[1].data

py/specter/psf/gausshermite2.py

@@ -51,7 +51,13 @@ def __init__(self, filename):
         #- Other necessary keywords
         self.npix_x = hdr['NPIX_X']
         self.npix_y = hdr['NPIX_Y']
-
+
+        #- PSF model error
+        if 'PSFERR' in hdr:
+            self.psferr = hdr['PSFERR']
+        else:
+            self.psferr = 0.01
+
         #- Load the parameters into self.coeff dictionary keyed by PARAM
         #- with values as TraceSets for evaluating the Legendre coefficients
         data = fx[1].data

py/specter/psf/psf.py

@@ -55,6 +55,12 @@ def __init__(self, filename):
         self.npix_y = hdr['NPIX_Y']
         self.nspec  = hdr['NSPEC']
 
+        #- PSF model error
+        if 'PSFERR' in hdr:
+            self.psferr = hdr['PSFERR']
+        else:
+            self.psferr = 0.01
+
         #- Load x, y legendre coefficient tracesets
         xc, hdr = fits.getdata(filename, 'XCOEFF', header=True)
         self._x = TraceSet(xc, domain=(hdr['WAVEMIN'], hdr['WAVEMAX']))

py/specter/util/util.py

@@ -293,9 +293,7 @@ def resample(x, xp, yp, xedges=False, xpedges=False):
 
     edges = x if xedges else pixspline.cen2bound(x)
     return ys.resample(edges)
-
-
-
+