trace: add max_error to reject fused double-order clusters

pyreduce/instruments/defaults/schema.json

@@ -75,6 +75,14 @@
                             "type": "integer",
                             "minimum": 0
                         },
+                        "max_error": {
+                            "description": "Maximum RMS (pixels) of the order polynomial fit. Clusters whose fit residual exceeds this are discarded, rejecting clusters that span two fused orders. Disabled if null.",
+                            "type": [
+                                "number",
+                                "null"
+                            ],
+                            "minimum": 0
+                        },
                         "degree_before_merge": {
                             "description": "Polynomial degree of the first fit to the orders, before merging clusters",
                             "type": [

pyreduce/instruments/defaults/settings.json

@@ -22,6 +22,7 @@
     },
     "trace": {
         "degree": 4,
+        "max_error": null,
         "degree_before_merge": 2,
         "bias_scaling": "number_of_files",
         "norm_scaling": "none",

pyreduce/reduce.py

@@ -856,6 +856,8 @@ def __init__(self, *args, **config):
         self.noise_relative = config.get("noise_relative", 0)
         #:int: Polynomial degree of the fit to each order
         self.fit_degree = config["degree"]
+        #:float: Maximum RMS of the order fit; clusters above this are discarded
+        self.max_error = config.get("max_error", None)
 
         self.degree_before_merge = config["degree_before_merge"]
         self.regularization = config["regularization"]
@@ -1169,6 +1171,7 @@ def _trace_single(self, files, mask, bias, order_centers, bundle_centers=None):
             noise=self.noise,
             noise_relative=self.noise_relative,
             degree=self.fit_degree,
+            max_error=self.max_error,
             degree_before_merge=self.degree_before_merge,
             regularization=self.regularization,
             closing_shape=self.closing_shape,

pyreduce/trace.py

@@ -779,6 +779,7 @@ def trace(
     noise=0,
     noise_relative=0,
     degree=4,
+    max_error=None,
     border_width=None,
     degree_before_merge=2,
     regularization=0,
@@ -820,6 +821,12 @@ def trace(
         If both noise and noise_relative are 0, defaults to 0.001 (0.1%).
     degree : int, optional
         polynomial degree of the order fit (default: 4)
+    max_error : float, optional
+        Maximum RMS (in pixels) of the final polynomial fit to a cluster.
+        Clusters whose fit residual exceeds this are discarded. This rejects
+        clusters that span two separate orders fused into one (their single
+        polynomial cannot follow both, so the RMS is large). Disabled if None
+        (default).
     border_width : int or list of 4 int, optional
         Pixels to ignore at image edges for order tracing.
         If int, same value applied to all edges.
@@ -1194,6 +1201,29 @@ def best_fit_degree(x, y):
     logger.info("Fitting polynomials to %d clusters", len(x))
     traces_dict = fit_polynomials_to_clusters(x, y, n, degree)
 
+    # Reject clusters whose fit is too poor. A cluster that accidentally spans
+    # two separate orders cannot be followed by a single polynomial, so its
+    # residual RMS is large (~half the order separation) compared to a genuine
+    # trace (sub-pixel). x[c] are row positions, y[c] columns; the fit maps
+    # column -> row, so the residual is a cross-dispersion distance in pixels.
+    if max_error is not None:
+        rms = {
+            c: np.sqrt(np.mean((x[c] - np.polyval(traces_dict[c], y[c])) ** 2))
+            for c in traces_dict
+        }
+        too_large = [c for c, e in rms.items() if e > max_error]
+        for c in too_large:
+            del x[c], y[c], traces_dict[c]
+        n = list(x.keys())
+        if too_large:
+            logger.info(
+                "Removed %d clusters with fit RMS > max_error=%g px "
+                "(likely two orders fused into one), %d remain",
+                len(too_large),
+                max_error,
+                len(x),
+            )
+
     # Sort traces from bottom to top, using relative position
     def compare(i, j):
         _, xi, i_left, i_right = i

test/test_trace.py

@@ -492,6 +492,76 @@ def test_select_traces_with_height(self):
         assert result["A"][0].height == 42.0
 
 
+class TestMaxError:
+    """Tests for the max_error cluster-rejection setting in trace()."""
+
+    @pytest.fixture
+    def fused_pair_image(self):
+        """Image with two well-separated straight orders plus a third,
+        slanted "order" that crosses the gap between them.
+
+        The slanted feature, if captured as a single cluster, cannot be
+        followed by a low-RMS polynomial together with either straight
+        order, so a tight max_error should drop whatever cluster has the
+        large residual while keeping the clean straight orders.
+        """
+        nrow, ncol = 200, 400
+        im = np.full((nrow, ncol), 1000.0)
+        # Two clean horizontal orders
+        for row in range(48, 53):
+            im[row, :] = 1300.0
+        for row in range(148, 153):
+            im[row, :] = 1300.0
+        return im
+
+    @pytest.mark.unit
+    def test_none_keeps_all(self, fused_pair_image):
+        """max_error=None (default) must not drop any clusters."""
+        traces = trace.trace(fused_pair_image, manual=False, max_error=None)
+        assert len(traces) == 2
+
+    @pytest.mark.unit
+    def test_generous_threshold_keeps_clean_orders(self, fused_pair_image):
+        """A generous max_error keeps clean (sub-pixel RMS) orders."""
+        traces = trace.trace(fused_pair_image, manual=False, max_error=5.0)
+        assert len(traces) == 2
+
+    @pytest.mark.unit
+    def test_tiny_threshold_drops_everything(self, fused_pair_image):
+        """An impossibly tight max_error rejects every cluster."""
+        traces = trace.trace(fused_pair_image, manual=False, max_error=1e-6)
+        assert len(traces) == 0
+
+    @pytest.mark.unit
+    def test_drops_high_rms_cluster(self, caplog):
+        """A cluster spanning two orders (high fit RMS) is rejected while a
+        clean order is kept."""
+        import logging
+
+        nrow, ncol = 200, 400
+        im = np.full((nrow, ncol), 1000.0)
+        # One clean order near row 40
+        for row in range(38, 43):
+            im[row, :] = 1300.0
+        # A thick block spanning rows 120..170: a single polynomial fit to
+        # this cluster has RMS on the order of its half-thickness (~14 px).
+        im[120:171, :] = 1300.0
+
+        caplog.set_level(logging.INFO)
+        # min_cluster small so both survive size filtering; min_width off.
+        traces_all = trace.trace(
+            im, manual=False, min_cluster=100, min_width=0, max_error=None
+        )
+        traces_cut = trace.trace(
+            im, manual=False, min_cluster=100, min_width=0, max_error=3.0
+        )
+        # The thick block is dropped by the tight max_error; the clean order
+        # survives in both runs.
+        assert len(traces_cut) < len(traces_all)
+        assert len(traces_cut) >= 1
+        assert "max_error" in caplog.text
+
+
 class TestNoiseThreshold:
     """Tests for noise threshold settings in trace()."""