Merge pull request #189 from GeminiDRSoftware/more-apfixing

More apfixing

geminidr/core/primitives_spect.py

@@ -2928,8 +2928,8 @@ def averaging_func(data, mask=None, variance=None):
 
                     # Recalculate aperture limits after rectification
                     apcoords = fit1d.evaluate(np.arange(ext.shape[dispaxis]))
-                    this_aptable["aper_lower"] = aperture["aper_lower"] + (location - apcoords.min())
-                    this_aptable["aper_upper"] = aperture["aper_upper"] - (apcoords.max() - location)
+                    this_aptable["aper_lower"] = aperture["aper_lower"]  #+ (location - apcoords.min())
+                    this_aptable["aper_upper"] = aperture["aper_upper"]  #- (apcoords.max() - location)
                     all_tables.append(this_aptable)
 
                 new_aptable = vstack(all_tables, metadata_conflicts="silent")

geminidr/gmos/tests/spect/test_trace_apertures.py

@@ -65,7 +65,8 @@ def test_regression_trace_apertures(ad, change_working_dir, ref_ad_factory):
         assert input_table['aper_lower'][0] <= 0
         assert input_table['aper_upper'][0] >= 0
 
-        keys = ext.APERTURE.colnames
+        # aper_lower and aper_upper aren't part of traceApertures
+        keys = [k for k in ext.APERTURE.colnames if not k.startswith("aper_")]
         actual = np.array([input_table[k] for k in keys])
         desired = np.array([reference_table[k] for k in keys])
         np.testing.assert_allclose(desired, actual, atol=0.05)

gempy/library/astrotools.py

@@ -253,6 +253,43 @@ def get_corners(shape):
 
     return corners
 
+
+def get_spline3_extrema(spline):
+    """
+    Find the locations of the minima and maxima of a cubic spline.
+
+    Parameters
+    ----------
+    spline: a callable spline object
+
+    Returns
+    -------
+    minima, maxima: 1D arrays
+    """
+    derivative = spline.derivative()
+    try:
+        knots = derivative.get_knots()
+    except AttributeError:  # for BSplines
+        knots = derivative.k
+
+    minima, maxima = [], []
+    # We take each pair of knots and map to interval [-1,1]
+    for xm, xp in zip(knots[:-1], knots[1:]):
+        ym, y0, yp = derivative([xm, 0.5*(xm+xp), xp])
+        # a*x^2 + b*x + c
+        a = 0.5 * (ym+yp) - y0
+        b = 0.5 * (yp-ym)
+        c = y0
+        for root in np.roots([a, b, c]):
+            if np.isreal(root) and abs(root) <= 1:
+                x = 0.5 * (root * (xp-xm) + (xp+xm))  # unmapped from [-1, 1]
+                if 2*a*root + b > 0:
+                    minima.append(x)
+                else:
+                    maxima.append(x)
+    return np.array(minima), np.array(maxima)
+
+
 def rotate_2d(degs):
     """
     Little helper function to return a basic 2-D rotation matrix.

gempy/library/tests/test_tracing.py

@@ -67,8 +67,8 @@ def test_get_limits():
 
     limits = tracing.get_limits(y, None, peaks=[CENT], threshold=0.01, method='peak')[0]
     for l in limits:
-        assert abs(l - CENT) - (SIG * np.sqrt(2 * np.log(100))) < 0.05 * SIG
+        assert abs(l - CENT) - (SIG * np.sqrt(2 * np.log(100))) < 0.1 * SIG
 
     limits = tracing.get_limits(y, None, peaks=[CENT], threshold=0.01, method='integral')[0]
     for l in limits:
-        assert abs(l - CENT) - (SIG * 2.576) < 0.05 * SIG
+        assert abs(l - CENT) - (SIG * 2.576) < 0.1 * SIG

gempy/library/tracing.py

@@ -155,7 +155,7 @@ def standard_extraction(self, data, mask, var, aper_lower, aper_upper):
             self.var[:] = [result.variance for result in results]
 
     def optimal_extraction(self, data, mask, var, aper_lower, aper_upper,
-                           cr_rej=5, max_iters=None):
+                           cr_rej=5, max_iters=None, degree=3):
         """Optimal extraction following Horne (1986, PASP 98, 609)"""
         BAD_BITS = DQ.bad_pixel | DQ.cosmic_ray | DQ.no_data | DQ.unilluminated
 
@@ -167,8 +167,9 @@ def optimal_extraction(self, data, mask, var, aper_lower, aper_upper,
         ix1 = max(int(min(all_x1) + 0.5), 0)
         ix2 = min(int(max(all_x2) + 1.5), slitlength)
 
-        fit_it = fitting.FittingWithOutlierRemoval(fitting.LinearLSQFitter(),
-                                                   outlier_func=sigma_clip, sigma_upper=3, sigma_lower=None)
+        fit_it = fitting.FittingWithOutlierRemoval(
+            fitting.LinearLSQFitter(), outlier_func=sigma_clip, sigma_upper=3,
+            sigma_lower=None)
 
         # If we don't have a VAR plane, assume uniform variance based
         # on the pixel-to-pixel variations in the data
@@ -178,9 +179,12 @@ def optimal_extraction(self, data, mask, var, aper_lower, aper_upper,
             var_mask = np.zeros_like(var, dtype=bool)
         else:
             mvar_init = models.Polynomial1D(degree=1)
-            var_model, var_mask = fit_it(mvar_init, np.ma.masked_where(mask.ravel(), abs(data).ravel()), var.ravel())
+            var_model, var_mask = fit_it(
+                mvar_init, np.ma.masked_where(mask.ravel(),
+                                              abs(data).ravel()), var.ravel())
             var_mask = var_mask.reshape(var.shape)[ix1:ix2]
             var = np.where(var_mask, var[ix1:ix2], var_model(data[ix1:ix2]))
+        var[var < 0] = 0
 
         if mask is None:
             mask = np.zeros((ix2 - ix1, npix), dtype=DQ.datatype)
@@ -191,26 +195,29 @@ def optimal_extraction(self, data, mask, var, aper_lower, aper_upper,
         # Step 4; first calculation of spectrum. We don't do any masking
         # here since we need all the flux
         spectrum = data.sum(axis=0)
-        weights = np.where(var > 0, var, 0)
+        inv_var = at.divide0(1., var)
         unmask = np.ones_like(data, dtype=bool)
 
         iter = 0
         while True:
             # Step 5: construct spatial profile for each wavelength pixel
             profile = np.divide(data, spectrum,
-                                out=np.zeros_like(data, dtype=np.float32), where=spectrum > 0)
+                                out=np.zeros_like(data, dtype=np.float32),
+                                where=spectrum > 0)
             profile_models = []
-            for row, wt_row in zip(profile, weights):
-                m_init = models.Chebyshev1D(degree=3, domain=[0, npix - 1])
-                m_final, _ = fit_it(m_init, pixels, row, weights=wt_row)
+            for row, ivar_row in zip(profile, inv_var):
+                m_init = models.Chebyshev1D(degree=degree, domain=[0, npix - 1])
+                m_final, _ = fit_it(m_init, pixels, row,
+                                    weights=np.sqrt(ivar_row) * np.where(spectrum > 0, spectrum, 0))
                 profile_models.append(m_final(pixels))
             profile_model_spectrum = np.array([np.where(pm < 0, 0, pm) for pm in profile_models])
             sums = profile_model_spectrum.sum(axis=0)
             model_profile = divide0(profile_model_spectrum, sums)
 
             # Step 6: revise variance estimates
-            var = np.where(var_mask | mask & BAD_BITS, var, var_model(abs(model_profile * spectrum)))
-            weights = divide0(1.0, var)
+            var = np.where(var_mask | mask & BAD_BITS, var,
+                           var_model(abs(model_profile * spectrum)))
+            inv_var = divide0(1.0, var)
 
             # Step 7: identify cosmic ray hits: we're (probably) OK
             # to flag more than 1 per wavelength
@@ -223,8 +230,8 @@ def optimal_extraction(self, data, mask, var, aper_lower, aper_upper,
 
             last_unmask = unmask
             unmask = (mask & BAD_BITS) == 0
-            spec_numerator = np.sum(unmask * model_profile * data * weights, axis=0)
-            spec_denominator = np.sum(unmask * model_profile ** 2 * weights, axis=0)
+            spec_numerator = np.sum(unmask * model_profile * data * inv_var, axis=0)
+            spec_denominator = np.sum(unmask * model_profile ** 2 * inv_var, axis=0)
             self.data = divide0(spec_numerator, spec_denominator)
             self.var = divide0(np.sum(unmask * model_profile, axis=0), spec_denominator)
             self.mask = np.bitwise_and.reduce(mask, axis=0)
@@ -693,21 +700,14 @@ def get_limits(data, mask, variance=None, peaks=[], threshold=0, method=None):
     else:
         w = divide0(1.0, np.sqrt(variance))
 
-    # We need to fit a quartic spline since we want to know its
-    # minima (roots of its derivative), and can only find the
-    # roots of a cubic spline
-    # TODO: Quartic splines look bad with outlier removal
-    #spline = astromodels.UnivariateSplineWithOutlierRemoval(x, y, w=w, k=4)
-    spline = interpolate.UnivariateSpline(x, y, w=w, k=4)
-
-    derivative = spline.derivative(n=1)
-    extrema = derivative.roots()
-    second_derivatives = spline.derivative(n=2)(extrema)
-    minima = [ex for ex, second in zip(extrema, second_derivatives) if second > 0]
+    # TODO: Consider outlier removal
+    #spline = am.UnivariateSplineWithOutlierRemoval(x, y, w=w, k=3)
+    spline = interpolate.UnivariateSpline(x, y, w=w, k=3)
+    minima, maxima = at.get_spline3_extrema(spline)
 
     all_limits = []
     for peak in peaks:
-        tweaked_peak = extrema[np.argmin(abs(extrema - peak))]
+        tweaked_peak = maxima[np.argmin(abs(maxima - peak))]
 
         # Now find the nearest minima above and below the peak
         for upper in minima: