apply single-line solution with warning

changes to code in wavecal plot to allow for this

Author: chris-simpson

geminidr/core/primitives_spect.py

@@ -2909,6 +2909,7 @@ def determineWavelengthSolution(self, adinputs=None, **params):
                     wavecal.update_wcs_with_solution(ext, fit1d, other, config)
             else:
                 bad_solutions = []
+                single_line_solutions = []
                 for ext in ad:
                     input_data, fit1d = wavecal.get_automated_fit(
                         ext, uiparams, p=self, linelist=linelist, bad_bits=DQ.not_signal,
@@ -2921,6 +2922,8 @@ def determineWavelengthSolution(self, adinputs=None, **params):
                             title=f"{ad.filename}:{ext.id}",
                             absorption=absorption))
                         any_good = True
+                        if fit1d.image.size == 1:
+                            single_line_solutions.append(ext.id)
                     else:  # no line matches so not an acceptable solution
                         bad_solutions.append(ext.id)
 
@@ -2935,6 +2938,16 @@ def determineWavelengthSolution(self, adinputs=None, **params):
                         raise RuntimeError(msg)
                     log.warning(msg)
 
+                if single_line_solutions:
+                    msg = f"{ad.filename}: a single-line solution has been applied"
+                    if len(ad) > 1:
+                        if len(ad) > len(single_line_solutions):
+                            msg += "in extensions " + ", ".join(str(i) for i in single_line_solutions)
+                        else:
+                            msg += " in all extensions"
+                    log.warning(msg)
+
+
             ad.update_filename(suffix=sfx, strip=True)
             if figures:
                 plot_filename = ad.filename.replace('.fits', '.pdf')

gempy/library/wavecal.py

@@ -814,7 +814,7 @@ def find_solution(init_models, config, peaks=None, peak_weights=None,
 
             # We allow this continue if we only matched 1 line but only 1 line is
             # expected
-            if len(matched) > 1: # or (len(matched) * nlines_expected == 1 and len(peaks) <= 2):
+            if len(matched) > 1 or (len(matched) * nlines_expected == 1 and len(peaks) <= 2):
                 m_init = models.Chebyshev1D(
                     degree=min(config["order"], (len(matched) + 1) // 2),
                     domain=domain)
@@ -829,59 +829,69 @@ def find_solution(init_models, config, peaks=None, peak_weights=None,
                 matched_peaks = peaks[matched]
                 matched_arc_lines = arc_lines[matches[matched]]
                 m_final = fit_it(m_init, matched_peaks, matched_arc_lines)
-                print(repr(m_final))
+                dw = abs(np.diff(m_final(m_final.domain))[0] / np.diff(m_final.domain)[0])
                 #for p, l in zip(matched_peaks, matched_arc_lines):
                 #    print(f"{p:.2f} => {l:.2f}")
 
-                # We're close to the correct solution, perform a KDFit
-                m_init = models.Chebyshev1D(m_final.degree, domain=domain)
-                for p, v in zip(m_final.param_names, m_final.parameters):
-                    setattr(m_init, p, v)
-                dw = abs(np.diff(m_final(m_final.domain))[0] / np.diff(m_final.domain)[0])
-                fit_it = matching.KDTreeFitter(sigma=2 * abs(dw), maxsig=5,
-                                               k=k, method='Nelder-Mead')
-                m_final = fit_it(m_init, peaks, arc_lines, in_weights=peak_weights,
-                                 ref_weights=arc_weights)
-                logit(f'{repr(m_final)} {fit_it.statistic}')
-
-                # And then recalculate the matches
-                match_radius = 4 * fwidth * abs(m_final.c1) / len_data  # 2*fwidth pixels
-                try:
-                    matched = matching.match_sources(m_final(peaks), arc_lines,
-                                                     radius=match_radius)
-                    incoords, outcoords = zip(*[(peaks[i], arc_lines[m])
-                                                for i, m in enumerate(matched) if m > -1])
-                    # Probably incoords and outcoords as defined here should go to
-                    # the interactive fitter, but cull to derive the "best" model
-                    fit1d = fit_1D(outcoords, points=incoords, function="chebyshev",
-                                   order=min(m_final.degree, (len(incoords) + 1) // 2),
-                                   domain=m_final.domain,
-                                   niter=config["niter"], sigma_lower=config["lsigma"],
-                                   sigma_upper=config["hsigma"])
-                    fit1d.image = np.asarray(outcoords)
-                except ValueError:
-                    log.warning("Line-matching failed")
-                    continue
-                nmatched = np.sum(~fit1d.mask)
-                logit(f"{filename} {repr(fit1d.model)} {nmatched} {fit1d.rms}")
-
-                # A posteriori check that the fit is within the bounds
-                # (LinearLSQFitter does not allow bounds so we can't force
-                # the model to be bounded).
-                try:
-                    for p, v in zip(fit1d.model.param_names, fit1d.model.parameters):
-                        bounds = model_bounds.get(p, (-np.inf, np.inf))
-                        new_bounds = (np.asarray(bounds) +
-                                      np.array([-0.5, 0.5]) * np.diff(bounds)[0])
-                        assert new_bounds[0] <= v <= new_bounds[1]
-                except AssertionError:
-                    for p, v in zip(fit1d.model.param_names, fit1d.model.parameters):
-                        bounds = model_bounds.get(p, (-np.inf, np.inf))
-                        new_bounds = (np.asarray(bounds) +
-                                      np.array([-0.5, 0.5]) * np.diff(bounds)[0])
-                        print(p, v, new_bounds)
-                    print("*** REJECTED ***")
-                    continue
+                if len(matched) > 1:
+                    # We're close to the correct solution, perform a KDFit
+                    m_init = models.Chebyshev1D(m_final.degree, domain=domain)
+                    for p, v in zip(m_final.param_names, m_final.parameters):
+                        setattr(m_init, p, v)
+                    fit_it = matching.KDTreeFitter(sigma=2 * abs(dw), maxsig=5,
+                                                   k=k, method='Nelder-Mead')
+                    m_final = fit_it(m_init, peaks, arc_lines, in_weights=peak_weights,
+                                     ref_weights=arc_weights)
+                    logit(f'{repr(m_final)} {fit_it.statistic}')
+
+                    # And then recalculate the matches
+                    match_radius = 4 * fwidth * abs(m_final.c1) / len_data  # 2*fwidth pixels
+                    try:
+                        matched = matching.match_sources(m_final(peaks), arc_lines,
+                                                         radius=match_radius)
+                        incoords, outcoords = zip(*[(peaks[i], arc_lines[m])
+                                                    for i, m in enumerate(matched) if m > -1])
+                        # Probably incoords and outcoords as defined here should go to
+                        # the interactive fitter, but cull to derive the "best" model
+                        fit1d = fit_1D(outcoords, points=incoords, function="chebyshev",
+                                       order=min(m_final.degree, (len(incoords) + 1) // 2),
+                                       domain=m_final.domain,
+                                       niter=config["niter"], sigma_lower=config["lsigma"],
+                                       sigma_upper=config["hsigma"])
+                        fit1d.image = np.asarray(outcoords)
+                    except ValueError:
+                        log.warning("Line-matching failed")
+                        continue
+                    nmatched = np.sum(~fit1d.mask)
+                    logit(f"{filename} {repr(fit1d.model)} {nmatched} {fit1d.rms}")
+
+                    # A posteriori check that the fit is within the bounds
+                    # (LinearLSQFitter does not allow bounds so we can't force
+                    # the model to be bounded). We allow the tolerances to be twice
+                    # as large as they were during the piecewise fit, as this seems
+                    # to work well empirically.
+                    try:
+                        for p, v in zip(fit1d.model.param_names, fit1d.model.parameters):
+                            bounds = model_bounds.get(p, (-np.inf, np.inf))
+                            new_bounds = (np.asarray(bounds) +
+                                          np.array([-0.5, 0.5]) * np.diff(bounds)[0])
+                            assert new_bounds[0] <= v <= new_bounds[1]
+                    except AssertionError:
+                        for p, v in zip(fit1d.model.param_names, fit1d.model.parameters):
+                            bounds = model_bounds.get(p, (-np.inf, np.inf))
+                            new_bounds = (np.asarray(bounds) +
+                                          np.array([-0.5, 0.5]) * np.diff(bounds)[0])
+                            print(p, v, new_bounds)
+                        continue
+                else:
+                    # Hack a fit1D object
+                    fit1d = fit_1D(np.arange(5), function="chebyshev", order=1,
+                                        niter=0)
+                    fit1d._models = m_final
+                    fit1d.image = np.array(matched_arc_lines)
+                    fit1d.points = np.array(matched_peaks)
+                    fit1d.mask = np.array([False], dtype=bool)
+                    nmatched = 1
 
                 # Wavelength solution models need to be monotonic. Make that check.
                 waves = fit1d.evaluate(np.arange(len_data))
@@ -1023,7 +1033,7 @@ def perform_piecewise_fit(model, peaks, arc_lines, pixel_start, kdsigma,
                     # automatically removes old (bad) match
                     matches[i] = m
                     found_new_matches = True
-                    print(f"Pixel {p} => was {m_init(p):.4f} now {m_this(p):.4f} {arc_lines[m]}")
+                    #print(f"Pixel {p} => was {m_init(p):.4f} now {m_this(p):.4f} {arc_lines[m]}")
         try:
             p_lo = peaks[matches > -1].min()
         except ValueError:
@@ -1304,14 +1314,15 @@ def create_pdf_plot(input_data, peaks, arc_lines, title="",
         data = input_data["spectrum"]
         spacing = 0.01
         vert_align = "bottom"
-    xmin, xmax = input_data["init_models"][0].domain
+    init_model = input_data["init_models"][0]
+    xmin, xmax = init_model.domain
     pixels = np.arange(xmin, xmax + 1)
     data_max = data.max()
     spacing *= data_max
     fig, ax = plt.subplots()
     ax.plot(pixels, data, 'b-')
     ax.set_ylim(0, data_max * 1.1)
-    if len(arc_lines) and np.diff(arc_lines)[0] / np.diff(peaks)[0] < 0:
+    if len(arc_lines) and init_model(xmin) > init_model(xmax):
         ax.set_xlim(xmax + 1, xmin - 1)
     else:
         ax.set_xlim(xmin - 1, xmax + 1)