Align only the columns that carry a trajectory, and compute the DTW cost per row

CHANGELOG.md

@@ -11,6 +11,42 @@ those changes.
 
 ## [Unreleased]
 
+## [1.93.0] - 2026-09-12
+
+### Changed
+
+- **`columns_to_align=None` now resolves to the columns that carry a trajectory**, not
+  every column of the first batch. Two kinds are skipped: columns that are not numeric,
+  and numeric columns that are flat through every batch. The second is what an identifier
+  column looks like, and a dtype test alone does not catch it, because the identifier
+  `melted_to_dict` leaves in place is usually an integer.
+
+  This was not harmless. On the dryer data `batch_id` was scaled as though it were a tag.
+  It took almost no weight itself (0.000051, since its forced range of 1.0 leaves the raw
+  identifiers deviating wildly from the average), but it joined the distance the alignment
+  minimises and the normalisation that follows, moving the real weights a long way:
+  measured over ten batches, `JacketTemperatureSP` went from 0.132 to 0.472.
+
+  A non-numeric column named explicitly in `columns_to_align` now raises, naming the
+  column and its dtype, rather than being silently dropped: the caller asked for it by
+  name, and whether to carry, encode or discard it is theirs to decide. A constant column
+  named explicitly is left alone, since constant over one set of batches does not mean
+  constant in general. (#199)
+
+- **The DTW cost matrix is computed per row instead of through a square product.**
+  `np.diag(A @ W @ A.T)` built an h-by-h product to keep its h diagonal entries, quadratic
+  in the reference length per test sample and cubic over the whole matrix. Summing
+  `(A @ W) * A` along the tags computes only the diagonal, in `h * J^2` work rather than
+  `h^2 * J`, with no large temporary. Measured unconstrained on random series: 2.08 ms to
+  1.49 ms at 100 samples, 24.98 to 11.69 at 300, 272.65 to 82.14 at 700, and 885.63 to
+  240.06 at 1200.
+
+  Results are unchanged in practice. Individual cost values differ by about 5e-16
+  relative, one unit in the last place, because the summation order changes; over nine
+  combinations of length and tag count the warping path is identical, and the dryer
+  alignment reproduces its final weights exactly.
+
+
 ## [1.92.0] - 2026-09-12
 
 ### Added
@@ -4987,7 +5023,8 @@ this entry records them together.
 - Reworked the README with a sharper value proposition and a
   "Why not scikit-learn?" comparison table.
 
-[Unreleased]: https://github.com/kgdunn/process-improve/compare/v1.92.0...HEAD
+[Unreleased]: https://github.com/kgdunn/process-improve/compare/v1.93.0...HEAD
+[1.93.0]: https://github.com/kgdunn/process-improve/compare/v1.92.0...v1.93.0
 [1.92.0]: https://github.com/kgdunn/process-improve/compare/v1.89.0...v1.92.0
 [1.89.0]: https://github.com/kgdunn/process-improve/compare/v1.88.0...v1.89.0
 [1.88.0]: https://github.com/kgdunn/process-improve/compare/v1.87.1...v1.88.0

CITATION.cff

@@ -12,7 +12,7 @@ authors:
 repository-code: "https://github.com/kgdunn/process-improve"
 url: "https://kgdunn.github.io/process-improve/"
 license: MIT
-version: 1.92.0
+version: 1.93.0
 date-released: "2026-09-12"
 keywords:
   - chemometrics

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "process-improve"
-version = "1.92.0"
+version = "1.93.0"
 description = 'Designed Experiments; Latent Variables (PCA, PLS, multivariate methods with missing data); Process Monitoring; Batch data analysis.'
 readme = "README.md"
 license = "MIT"

src/process_improve/batch/alignment_helpers.py

@@ -427,23 +427,25 @@ def _banded_distance_matrix(
 
     # Mahalanobis distance, computed only where the band admits it. Filling every cell
     # first and then constraining only the accumulation would leave the whole function
-    # quadratic in the two batch lengths however narrow the band: on random 700-sample
-    # series a 10% band was then 7x faster, against 359x once the cost is restricted too.
+    # quadratic in the two batch lengths however narrow the band.
     #
-    # The two cases are separate loops rather than one loop over runtime bounds. The
-    # bounds are monotone, so these two entries decide whether the band spans everything;
-    # when it does, the original whole-array expression is used. A slice of the same
-    # extent taken with runtime bounds is not free: numba cannot prove it contiguous and
-    # compiles a slower matmul, which cost the unconstrained default 4x at 700 samples.
-    if band[nt - 1, 0] == 0 and band[0, 1] == nr:
-        for idx in np.arange(nt):
-            deviation = test[idx] - ref
-            dist[:, idx] = np.diag(deviation @ weight_matrix @ deviation.T)
-    else:
-        for idx in np.arange(nt):
-            lower, upper = band[idx, 0], band[idx, 1]
-            deviation = test[idx] - ref[lower:upper]
-            dist[lower:upper, idx] = np.diag(deviation @ weight_matrix @ deviation.T)
+    # The per-row form, not `np.diag(A @ W @ A.T)`: that builds an h-by-h product to keep
+    # its h diagonal entries, so it is quadratic in the number of reference rows and cubic
+    # over the whole matrix. `(A @ W) * A` summed along the tags computes only the
+    # diagonal, in `h * J^2` work rather than `h^2 * J`, and needs no large temporary. It
+    # also removes the contiguity problem the matmul had, where a slice taken with runtime
+    # bounds compiled to a slower kernel than the whole array, so one loop now serves both
+    # the banded and the unconstrained case.
+    #
+    # The summation order differs from the matmul's, so individual cost values move in the
+    # last bits: measured at about 5e-16 relative, one unit in the last place. That is
+    # below the resolution of anything downstream. Over nine combinations of length and
+    # tag count the warping path came back identical, and the dryer alignment reproduces
+    # its weights exactly, so no pinned value moved.
+    for idx in np.arange(nt):
+        lower, upper = band[idx, 0], band[idx, 1]
+        deviation = test[idx] - ref[lower:upper]
+        dist[lower:upper, idx] = np.sum((deviation @ weight_matrix) * deviation, axis=1)
 
     D = np.zeros((nr, nt)) * np.nan
     D[0, 0] = dist[0, 0]

src/process_improve/batch/preprocessing.py

@@ -28,6 +28,25 @@
 epsqrt = np.sqrt(np.finfo(float).eps)
 
 
+def _varies_within_any_batch(batches: dict[str, pd.DataFrame], column: object) -> bool:
+    """
+    Report whether ``column`` takes more than one value inside at least one batch.
+
+    A column that is flat through every batch carries no trajectory for the alignment to
+    work with. Every batch is checked rather than only the first, so a tag that happens to
+    be flat in the first batch but moves in a later one is kept.
+    """
+    return any(_has_spread(batch[column].to_numpy()) for batch in batches.values() if column in batch)
+
+
+def _has_spread(values: np.ndarray) -> bool:
+    """Report whether a numeric array holds more than one distinct usable value."""
+    usable = values[np.isfinite(values)]
+    # min against max rather than counting distinct values: one pass, no hashing, and the
+    # question is only whether the column moves at all.
+    return bool(usable.size and usable.min() != usable.max())
+
+
 def _resolve_columns_to_align(
     batches: dict[str, pd.DataFrame],
     columns_to_align: list | pd.Index | None,
@@ -48,7 +67,8 @@ def _resolve_columns_to_align(
     batches : dict[str, pd.DataFrame]
         Batch data, in the standard format (keyed by batch identifier).
     columns_to_align : list, pd.Index, or None
-        Passed through when given; resolved from the first batch when ``None``.
+        Checked when given; resolved from the numeric columns of the first batch when
+        ``None``.
     caller : str
         Name of the calling function, used in the error messages.
 
@@ -62,7 +82,31 @@ def _resolve_columns_to_align(
     TypeError
         If ``batches`` is a DataFrame, or is not a dict.
     ValueError
-        If ``batches`` is an empty dict, so there is no batch to take columns from.
+        If ``batches`` is an empty dict, so there is no batch to take columns from; if a
+        batch holds no numeric column at all; or if an explicit ``columns_to_align``
+        names a column that is not numeric.
+
+    Notes
+    -----
+    Only numeric columns are aligned. Resolving from every column swept in whatever else
+    a batch frame carried, most often the batch identifier that
+    :func:`~process_improve.batch.data_input.melted_to_dict` leaves in place: constant
+    within a batch, so it contributed a zero range, and meaningless to scale. A column
+    that is not numeric is not a trajectory, and what to do with it (carry it, drop it,
+    encode it) is the caller's decision, not something to guess at here.
+
+    A numeric column that holds one value throughout every batch is dropped for the same
+    reason. That is what an identifier column looks like, and the identifier
+    :func:`~process_improve.batch.data_input.melted_to_dict` leaves in place is usually a
+    number, so a dtype test alone does not catch it. Including it is not harmless: on the
+    dryer data it takes almost no weight itself (0.000051) but moves the others
+    substantially, ``JacketTemperatureSP`` from 0.132 to 0.472, because it joins the
+    distance the alignment minimises and the normalisation that follows.
+
+    An explicit ``columns_to_align`` naming a non-numeric column raises rather than
+    quietly dropping it, since the caller asked for it by name. A constant column named
+    explicitly is left alone: constant over this particular set of batches does not mean
+    constant in general, and the caller may know better.
     """
     if isinstance(batches, pd.DataFrame):
         raise TypeError(
@@ -78,16 +122,38 @@ def _resolve_columns_to_align(
             f"identifier; got {type(batches).__name__}."
         )
 
-    if columns_to_align is not None:
-        return columns_to_align
-
     if not batches:
         raise ValueError(
             f"{caller} cannot resolve `columns_to_align` from an empty `batches` dict; "
             "pass `columns_to_align` explicitly, or supply at least one batch."
         )
 
-    return batches[next(iter(batches))].columns
+    first_batch = batches[next(iter(batches))]
+    if columns_to_align is not None:
+        non_numeric = {
+            str(column): str(first_batch[column].dtype)
+            for column in columns_to_align
+            if column in first_batch and not pd.api.types.is_numeric_dtype(first_batch[column])
+        }
+        if non_numeric:
+            listed = ", ".join(f"{name!r} ({dtype})" for name, dtype in sorted(non_numeric.items()))
+            raise ValueError(
+                f"{caller} can only align numeric columns, but `columns_to_align` names: {listed}. "
+                "Drop them from `columns_to_align`; a non-numeric column is not a trajectory, and "
+                "carrying, encoding or discarding it is yours to decide."
+            )
+        return columns_to_align
+
+    numeric = first_batch.select_dtypes(include="number").columns
+    varying = [column for column in numeric if _varies_within_any_batch(batches, column)]
+    if not varying:
+        raise ValueError(
+            f"{caller} found no column to align in the first batch, whose columns are "
+            f"{list(first_batch.columns)}. A column must be numeric and must vary within at least "
+            "one batch. Pass `columns_to_align` explicitly, or supply batches holding at least one "
+            "numeric trajectory."
+        )
+    return varying
 
 
 #: Quantile pair behind each ``settings["robust_range"]`` choice in :func:`determine_scaling`.

tests/batch/test_preprocessing.py

@@ -877,14 +877,24 @@ class TestDetermineScalingZeroRangeWarning:
     value for over half a batch has no interquartile spread at all.
     """
 
-    def test_a_constant_tag_is_reported_once_for_the_whole_call(self, dryer_data: dict) -> None:
+    def test_a_tag_flat_in_some_batches_is_reported_once_for_the_whole_call(self, dryer_data: dict) -> None:
         with pytest.warns(UserWarning, match="zero range and substituted 1.0") as caught:
-            determine_scaling(dryer_data)  # every column, including the constant batch_id
+            determine_scaling(dryer_data)
 
         assert len(caught) == 1, "one aggregated warning, not one per batch"
-        message = str(caught[0].message)
-        assert "'batch_id' (71 of 71 batches)" in message
-        assert "'DifferentialPressure' (16 of 71 batches)" in message
+        assert "'DifferentialPressure' (16 of 71 batches)" in str(caught[0].message)
+
+    def test_the_identifier_column_is_not_reported_because_it_is_not_aligned(self, dryer_data: dict) -> None:
+        """It used to head this warning on every default call, before it was excluded.
+
+        `batch_id` is flat in every batch, so column resolution now drops it and it never
+        reaches the zero-range substitution. `DifferentialPressure` is the real case: flat
+        in 16 of the 71 batches and varying in the rest, so it is aligned and reported.
+        """
+        with pytest.warns(UserWarning, match="zero range and substituted 1.0") as caught:
+            determine_scaling(dryer_data)
+
+        assert "batch_id" not in str(caught[0].message)
 
     def test_tags_that_vary_produce_no_warning(self, dryer_data: dict) -> None:
         with warnings.catch_warnings():
@@ -1216,3 +1226,140 @@ def test_the_settings_reach_determine_scaling(self, dryer_data: dict) -> None:
 
         assert (iqr.scale_df_["Range"] < wide.scale_df_["Range"]).all()
         assert not np.array_equal(raw.scale_df_["Range"].to_numpy(), wide.scale_df_["Range"].to_numpy())
+
+
+class TestCostMatrixDiagonal:
+    """The cost matrix computes only the diagonal it needs (#199 follow-up)."""
+
+    @staticmethod
+    def _textbook(test: np.ndarray, ref: np.ndarray, weights: np.ndarray) -> np.ndarray:
+        """Build the cost matrix from the literal Mahalanobis expression, as a reference."""
+        columns = [np.diag((row - ref) @ weights @ (row - ref).T) for row in test]
+        return np.column_stack(columns)
+
+    def test_it_agrees_with_the_literal_mahalanobis_expression(self) -> None:
+        rng = np.random.default_rng(0)
+        test, ref = rng.normal(size=(40, 5)), rng.normal(size=(47, 5))
+        weights = np.diag(rng.uniform(0.5, 2.0, 5))
+
+        accumulated = distance_matrix(test, ref, weights)
+        # `distance_matrix` returns the accumulated cost, so compare the entry that has
+        # not yet been accumulated into: the top-left corner is the raw cost there.
+        expected = self._textbook(test, ref, weights)
+        assert accumulated[0, 0] == pytest.approx(expected[0, 0], rel=1e-12)
+        # And the first row accumulates left to right over the raw costs.
+        assert accumulated[0, :].tolist() == pytest.approx(np.cumsum(expected[0, :]).tolist(), rel=1e-12)
+
+    @pytest.mark.parametrize(("n_test", "n_ref", "n_tags"), [(50, 57, 3), (200, 193, 5), (120, 120, 11)])
+    def test_the_warping_path_is_unchanged_by_the_row_wise_form(self, n_test: int, n_ref: int, n_tags: int) -> None:
+        """The per-row diagonal sums in a different order, so this pins that it does not matter."""
+        rng = np.random.default_rng(n_test)
+        test, ref = rng.normal(size=(n_test, n_tags)), rng.normal(size=(n_ref, n_tags))
+        weights = np.diag(rng.uniform(0.5, 2.0, n_tags))
+
+        expected_cost = self._textbook(test, ref, weights)
+        accumulated = np.full((n_ref, n_test), np.nan)
+        accumulated[0, 0] = expected_cost[0, 0]
+        accumulated[0, 1:] = np.cumsum(expected_cost[0, :])[1:]
+        accumulated[1:, 0] = np.cumsum(expected_cost[:, 0])[1:]
+        for column in range(1, n_test):
+            for row in range(1, n_ref):
+                accumulated[row, column] = expected_cost[row, column] + np.nanmin(
+                    [accumulated[row, column - 1], accumulated[row - 1, column - 1], accumulated[row - 1, column]]
+                )
+
+        by_kernel, _ = backtrack_optimal_path(distance_matrix(test, ref, weights))
+        by_textbook, _ = backtrack_optimal_path(accumulated)
+        assert np.array_equal(by_kernel, by_textbook)
+
+    def test_it_does_not_build_the_square_product(self) -> None:
+        """A tall reference against few tags would allocate gigabytes under the old form."""
+        rng = np.random.default_rng(1)
+        test, ref = rng.normal(size=(30, 2)), rng.normal(size=(4000, 2))
+        weights = np.eye(2)
+
+        accumulated = distance_matrix(test, ref, weights)
+
+        assert accumulated.shape == (4000, 30)
+        assert np.isfinite(accumulated[-1, -1])
+
+
+class TestColumnResolution:
+    """Only columns that carry a trajectory are aligned (#199)."""
+
+    @staticmethod
+    def _batches(**columns: object) -> dict:
+        """Two batches, each holding the given columns; scalars become constant columns."""
+        return {
+            batch_id: pd.DataFrame(
+                {name: (value if isinstance(value, list) else [value] * 4) for name, value in columns.items()}
+            )
+            for batch_id in (1, 2)
+        }
+
+    def test_a_non_numeric_column_is_skipped_by_default(self) -> None:
+        batches = self._batches(temp=[10.0, 11.0, 12.0, 13.0], operator="alice")
+
+        scale_df = determine_scaling(batches)
+
+        assert list(scale_df.index) == ["temp"]
+
+    def test_a_non_numeric_column_named_explicitly_raises(self) -> None:
+        """The caller asked for it by name, so dropping it silently would hide the mistake."""
+        batches = self._batches(temp=[10.0, 11.0, 12.0, 13.0], operator="alice")
+
+        with pytest.raises(ValueError, match=r"can only align numeric columns.*'operator' \(object\)"):
+            determine_scaling(batches, columns_to_align=["temp", "operator"])
+
+    def test_a_column_flat_in_every_batch_is_skipped(self) -> None:
+        """An identifier looks exactly like this, and it is numeric, so dtype alone misses it."""
+        batches = self._batches(temp=[10.0, 11.0, 12.0, 13.0], batch_id=7)
+
+        scale_df = determine_scaling(batches)
+
+        assert list(scale_df.index) == ["temp"]
+
+    def test_a_column_flat_in_only_the_first_batch_is_kept(self) -> None:
+        """Every batch is checked, so a tag that starts flat but moves later still counts."""
+        batches = {
+            1: pd.DataFrame({"temp": [10.0, 11.0, 12.0], "pressure": [5.0, 5.0, 5.0]}),
+            2: pd.DataFrame({"temp": [10.0, 11.0, 12.0], "pressure": [5.0, 6.0, 7.0]}),
+        }
+
+        scale_df = determine_scaling(batches)
+
+        assert sorted(scale_df.index) == ["pressure", "temp"]
+
+    def test_a_constant_column_named_explicitly_is_left_alone(self) -> None:
+        """Constant over these batches does not mean constant in general; the caller may know."""
+        batches = self._batches(temp=[10.0, 11.0, 12.0, 13.0], setpoint=50.0)
+
+        with pytest.warns(UserWarning, match="zero range"):
+            scale_df = determine_scaling(batches, columns_to_align=["temp", "setpoint"])
+
+        assert sorted(scale_df.index) == ["setpoint", "temp"]
+
+    def test_no_alignable_column_raises_and_lists_what_was_there(self) -> None:
+        batches = self._batches(label="a", identifier=3)
+
+        with pytest.raises(ValueError, match=r"found no column to align.*\['label', 'identifier'\]"):
+            determine_scaling(batches)
+
+    def test_an_all_nan_column_counts_as_flat(self) -> None:
+        """No usable value means no spread, so it carries nothing to align."""
+        batches = self._batches(temp=[10.0, 11.0, 12.0, 13.0], broken=[np.nan] * 4)
+
+        scale_df = determine_scaling(batches)
+
+        assert list(scale_df.index) == ["temp"]
+
+    def test_a_column_absent_from_some_batches_is_judged_on_the_rest(self) -> None:
+        batches = {
+            1: pd.DataFrame({"temp": [10.0, 11.0, 12.0]}),
+            2: pd.DataFrame({"temp": [10.0, 11.0, 12.0], "extra": [1.0, 2.0, 3.0]}),
+        }
+
+        # `extra` is absent from the first batch, so it is not in the resolved set, which
+        # comes from that batch's columns. The point is that checking spread across every
+        # batch does not raise on the batch where the column is missing.
+        assert list(determine_scaling(batches).index) == ["temp"]