Repo-wide error audit: fix statistical and robustness bugs (#500)

* chore: start repo-wide error audit and fix series

Ruthless audit of statistical correctness and robustness across the
package. Fixes land as micro-commits on this branch; the PR description
carries the ranked findings.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EM4uAK1eM5YqBtpoLsSseM

* fix(multivariate): correct statistical errors found in the repo audit

- TSR PCA: the N <= K SVD branch returned an already-oriented loading
  matrix which the shared transpose then mangled: an IndexError for
  N < K and a silently wrong imputation regression for N == K. Fit
  scores were also centred while transform() projects uncentred, so
  scores_, SPE and R2 disagreed with transform on the same data; the
  sign convention is now applied like the SVD/NIPALS paths, and the
  EM loop is skipped entirely for complete data.
- Score-plot T2 ellipse: use the bivariate limit (2 degrees of
  freedom) instead of the full model's A; the old ellipse was ~42
  percent too wide per axis at A=5, N=50, hiding genuine outliers.
- spe_plot / t2_plot: compute the confidence limit at the plotted
  component count instead of always the last component, and restore
  the y-axis title (it previously showed the limit legend text).
- PCA.select_n_components: the 1-SE band compared total PRESS against
  a per-fold standard error (~n_folds too narrow, degenerating the
  1se rule to min); the Q2 null model now uses the centred sum of
  squares instead of sum(x^2).
- PLS.select_n_components: same n_folds rescaling for the Q2 SE band.
- PLS.cross_validate: K-fold beta confidence intervals now use the
  delete-a-block jackknife standard error (the plain sample SD was
  (K-1)/sqrt(K) times too small); Q2 uses nanmean for the Y centre.
- PLS/PCA NIPALS: the max-iterations warning could never fire (itern
  is capped AT the maximum); PCA previously had no warning at all.
- Target projection / selectivity ratio: the projection direction now
  uses the scaled-space regression vector and maps X through the
  model's own scaler; the raw-units beta_coefficients_ vector is not
  a direction in the internal space when scale=True (the default).
- TPLS.diagnose: zero out missing cells after building the presence
  maps, as fit() does; NaN * 0 is NaN, so one missing F/Z cell
  previously poisoned that observation's scores, T2, SPE and
  predictions.
- Hotelling's T2 in fit(): skip components with ~zero score variance
  instead of dividing by ~zero (rank-deficient fits produced inf/NaN
  T2 for every observation); validate n_components >= 1.
- MCUVScaler: a column with fewer than two observed values has NaN
  nanstd which the ==0 guard missed, emitting an all-NaN column on
  transform; non-finite centres/scales are now treated as constant.
- center()/scale(): axis=1 broadcast the row statistic across
  columns (ValueError for rectangular input, silently wrong for
  square); the statistic is now reshaped to a column vector.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EM4uAK1eM5YqBtpoLsSseM

* fix(multivariate): satisfy mypy and CodeQL on the audit fixes

Use a dedicated DataFrame variable in _target_projection_arrays (mypy
union-attr), and hoist the ekf PRESS scale multiplier out of the branch
so CodeQL cannot see an uninitialized local.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EM4uAK1eM5YqBtpoLsSseM

* fix(univariate, monitoring): correct statistical errors found in the repo audit

- Generalized ESD: the outlier count is the LARGEST i with R_i >
  lambda_i (NIST/Rosner); the code took the first crossing, which
  under-reports exactly in the masking scenarios the test exists for.
- Generalized ESD: robust_variant now defaults to False. The MAD-scaled
  statistic has no upper bound while the critical values are derived
  for the mean/std statistic (bounded by (N-1)/sqrt(N)), so the robust
  variant declares outliers in clean data; it stays available as an
  explicitly documented screening heuristic.
- Robust confidence interval (metrics and the agent tool): the interval
  is for the median, whose asymptotic standard error is
  sigma*sqrt(pi/2)/sqrt(n); the missing factor gave ~87% coverage for
  a nominal 95% interval.
- variance_decomposition: between_stddev now reports the between-group
  variance COMPONENT sqrt((MS_between - MS_within)/n0) instead of
  sqrt(MS_between), which mixed the within-group noise into the
  between number (the docstring example itself showed the wrong value).
- biweight_midvariance: use the midvariance tuning constant c = 9;
  c = 6 is the biweight location constant and biased the scale low.
- Holt-Winters chart: the biweight rho conflated the consistency
  constant with the cutoff k = 2.52, so every scale estimate was 12%
  too small and the +/-3S limits were really +/-2.63 sigma (~3x the
  nominal false-alarm rate). Warm-up residuals now subtract the trend
  beta_0*t rather than the constant beta_0. The lambda grid search is
  NaN-aware (row 0 has no error value, so for 10 <= N < 20 every grid
  cell was NaN and (0.1, 0.1) always won silently). An explicit
  ld_1=0.0 is respected instead of being treated as unset. Unknown
  chart variants are rejected at construction with a clear message.
- The agent-facing control_chart tool no longer advertises a CUSUM
  chart type that always failed with a misleading error; the package
  docstring's CUSUM/EWMA claim is corrected too.
- calculate_cpk: rsd is now the relative standard deviation of the
  data (spread over the data centre), not spread over the
  distance-to-spec centre, which changed value when the spec moved;
  the docstring documents that the overall-sigma statistic is
  Ppk-style. The capability tool reports an undefined Cpk as 'could
  not be computed' instead of 'Poor capability'.
- Residual diagnostics: p-values that underflow to exactly 0.0 are the
  most significant result possible; use 'is not None' instead of
  truthiness so they are no longer rendered as unavailable.
- Test suite: pins that encoded the pre-fix constants are updated with
  derivations in comments; the registry test no longer depends on
  sibling tests having run on the same xdist worker.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EM4uAK1eM5YqBtpoLsSseM

* style: satisfy ruff on the univariate and monitoring audit fixes

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EM4uAK1eM5YqBtpoLsSseM

* docs: fix RST list structure in the variance_decomposition docstring

The Note block was never a real RST list (no blank line after the
heading), so the new multi-line bullet's continuation line failed the
strict Sphinx build. Promote it to a proper NumPy-style Notes section.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EM4uAK1eM5YqBtpoLsSseM

* fix(experiments): correct DOE errors found in the repo audit

- Clear effects now follow Wu & Hamada: an effect is clear only when
  every alias has order >= 3. The previous 'higher order than the
  effect' rule declared every main effect of a resolution-III design
  clear (A = BC has order 2 > 1), the exact case the concept exists to
  flag.
- Explicit fractional-factorial generators: pyDOE3 returns columns as
  (bases..., derived...) while the caller assigns positionally to the
  factor list, so a generator on a non-last factor (B=AC) silently
  swapped factor columns; raw factor names were also lower-cased into
  pyDOE3's single-letter notation, so multi-character names were
  misread as products of letters. Generators are now parsed against
  the real factor names (same convention as evaluate._parse_word, but
  raising on unparseable content), translated to canonical letters,
  and the columns re-ordered back to the caller's factor order.
  Negative generators (D=-ABC) are supported and inconsistent
  generator sets are rejected with clear messages.
- Column.to_coded / to_realworld: an explicit center=0 (falsy) was
  silently replaced by the stored pi_center; missing or zero-width
  ranges now raise a clear ValueError instead of TypeError / silent
  inf.
- gather(): positional arguments were accepted by the signature and
  silently discarded; they are now folded in via their own column
  names, and a nameless positional argument raises.
- D-optimal point exchange: the scorer de-duplicated the design before
  computing |X'X| (replicated runs carry real information); an
  improving swap onto the row with index label 0 was discarded by a
  truthiness test; the shuffle now takes a random_state for
  reproducibility.
- Lack-of-fit test: replicate groups are found on the MODEL's factor
  columns with rounded numeric values. Grouping on the whole frame
  meant the unique-per-row RunOrder column made every group a
  singleton, so no generated design ever had detectable replicates
  and the test always reported 'No replicated points'.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EM4uAK1eM5YqBtpoLsSseM

* fix(regression, batch): correct numerical errors found in the repo audit

- Robust regression: the degenerate-x guard branch itself divided the
  ~zero x deviations by the ~zero x sum-of-squares, poisoning leverage
  (and influence) with NaN/inf; with no variation in x the leverage is
  exactly 1/N.
- DTW alignment: the reported 'distance' summed the CUMULATIVE cost
  matrix entries along the warping path (a sum of prefix sums that
  grows super-linearly with path length); the DTW distance is the
  accumulated cost D[-1, -1]. The normalized distance and the
  per-batch alignment-quality numbers inherit the fix.
- Kassidas batch alignment weights: a variable whose trajectories
  align near-perfectly (SSQ ~ 0) must receive a LARGE weight (weights
  are inversely proportional to the SSQ); the previous guard
  substituted the scale-dependent magic value 10000 for a near-zero
  SSQ, giving the best-aligned variables a weight of ~1e-4, the exact
  opposite. The SSQ is now floored relative to the largest observed
  SSQ.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EM4uAK1eM5YqBtpoLsSseM

* fix(infra): config caching, JSON serialisation, tool discovery, pool thread-safety, CI hygiene

- Settings: the setdefault pattern evaluated the env read on every
  access, so nothing was actually cached and a knob served successfully
  could raise later if the env var went bad; a _get helper now caches
  genuinely on first access. Numeric knobs must be positive, and
  boolean env vars reject unrecognized values instead of silently
  reading them as false (a typo in PROCESS_IMPROVE_MCP_SAFE_MODE
  previously disabled the security-relevant safe mode with no error).
- clean(): handle np.bool_ (a subclass of neither np.integer nor
  bool), numpy scalar dict keys (pandas groupby labels), sets, and the
  remaining numpy scalar types via np.generic.item(); each previously
  surfaced as a generic internal error at the MCP boundary.
- discover_tools: only tolerate a ModuleNotFoundError whose missing
  module is third-party; a typo'd or renamed first-party module now
  propagates instead of silently dropping a whole tool category with a
  'missing dependency' warning.
- tool_safety: the module-level worker pool was created and torn down
  with no lock while the MCP server calls tools from executor threads;
  one thread's teardown could SIGKILL the worker running another
  thread's task (mis-reported as a memory-limit kill) or leak an
  orphaned worker. The default path now runs each call in a private
  per-call pool (same cost as the old per-call recycling, no shared
  state), and the remaining module-pool helpers are lock-guarded.
- tests/fuzz: the boundary fuzzer now excludes test-only tools
  (leading underscore). Depending only on import order, the registry
  snapshot could include test_tool_safety's deliberate infinite-loop
  tool and the fuzzer then ran it in-process with no timeout,
  hanging the run; this reproduced locally during this audit.
- raincloud: without the plotting extra, raise the documented
  'install the extra' ImportError at the call site instead of an
  AttributeError from the module stub.
- CI: drop the no-op create trigger (the create event ignores
  branch/tag filters, so the full matrix ran on every branch
  creation); grant run-tests contents:read only; move the Pages
  deploy scopes off the docs build job, which executes PR code.
- Version 1.67.0: many of the audit fixes change numerical results
  (and one default), so this is a MINOR bump. The changelog also gains
  the missing 1.66.2 section: pyproject and CITATION already claimed
  1.66.2 while its entries still sat under Unreleased, which would
  have made the tag-gated release notes extraction silently fall back
  to auto-generated notes.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EM4uAK1eM5YqBtpoLsSseM

* refactor(tool_safety): fold the pool globals into one state tuple

Addresses the CodeQL alerts on the previous commit: the split
_pool/_pool_memory_mb globals read as an unused variable to the
scanner, and the test-local module import mixed import styles. A
single _pool_state tuple is also harder to update inconsistently.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EM4uAK1eM5YqBtpoLsSseM

---------

Co-authored-by: Claude <noreply@anthropic.com>

Author: kgdunn

.github/workflows/docs.yml

@@ -7,10 +7,12 @@ on:
     branches: ['**']
   workflow_dispatch:  # allow manual trigger
 
+# Only the read scope at workflow level. The Pages-deploy scopes are
+# granted to the `deploy` job alone: the `build` job executes arbitrary
+# dependency and notebook code on every same-repo PR and must not carry
+# deploy credentials while doing so.
 permissions:
   contents: read
-  pages: write
-  id-token: write
 
 concurrency:
   group: ${{ github.workflow }}-${{ github.event.pull_request.number || github.ref }}
@@ -19,6 +21,9 @@ concurrency:
 jobs:
   build:
     runs-on: ubuntu-latest
+    permissions:
+      contents: read
+      pages: read  # configure-pages (main-branch pushes only) queries the Pages config
     steps:
       - uses: actions/checkout@v7
       - name: Install uv
@@ -44,6 +49,9 @@ jobs:
     if: github.event_name == 'push' && github.ref == 'refs/heads/main'
     needs: build
     runs-on: ubuntu-latest
+    permissions:
+      pages: write
+      id-token: write
     environment:
       name: github-pages
       url: ${{ steps.deployment.outputs.page_url }}

.github/workflows/run-tests.yml

@@ -3,17 +3,24 @@
 
 name: Unit tests
 
+# NOTE: no `create` trigger. The create event does not support
+# branches/tags filters (they are silently ignored), so the old
+#   create: {branches: [main], tags: ['**']}
+# ran the full matrix on EVERY branch creation, duplicating the
+# push / pull_request runs.
 on:
   push:
     branches: [main]
   pull_request:
     branches: ['**']
-  create:
-    branches: [main]
-    tags: ['**']
   # schedule:
   #   - cron: "0 4 * * *"
 
+# The workflow only reads the repository; the default token needs no
+# write scopes (codecov uploads use their own OIDC/token path).
+permissions:
+  contents: read
+
 concurrency:
   group: ${{ github.workflow }}-${{ github.event.pull_request.number || github.ref }}
   # Only cancel redundant in-progress runs for pull requests. Pushes to main

CHANGELOG.md

@@ -11,6 +11,145 @@ those changes.
 
 ## [Unreleased]
 
+## [1.68.0] - 2026-08-22
+
+A repo-wide correctness audit. Most entries below change numerical results
+because the previous values were statistically or numerically wrong; the
+regression tests in `tests/test_audit_regressions_*.py` pin each fix.
+
+### Fixed
+
+- PCA with `algorithm="tsr"`: the missing-data imputation used a wrongly
+  transposed singular-vector matrix whenever the data had at least as many
+  columns as rows (an `IndexError` for wide matrices, silently wrong imputation
+  for square ones), and the fitted scores were centred while `transform()`
+  projects uncentred, so `scores_`, SPE and R2 disagreed with `transform` on
+  the same data. TSR loadings now also follow the same sign convention as the
+  SVD and NIPALS fits.
+- The score-plot T2 confidence ellipse (`ellipse_coordinates`) now uses the
+  bivariate limit (2 degrees of freedom) instead of the full model's component
+  count. The old ellipse was ~42 percent too wide per axis for a 5-component
+  model on 50 observations, hiding genuine outliers.
+- `spe_plot(with_a=...)` and `t2_plot(with_a=...)` now draw the confidence
+  limit computed at the plotted component count instead of always the full
+  model's; the y-axis title no longer shows the limit legend text.
+- `PCA.select_n_components`: the 1-SE selection band compared the total PRESS
+  against a per-fold standard error (about n_folds times too narrow, silently
+  degenerating the `"1se"` rule to `"min"`), and the Q2 null model used the
+  uncentred sum of squares. `PLS.select_n_components` had the same too-narrow
+  Q2 standard-error band.
+- `PLS.cross_validate` with K-fold resampling: beta confidence intervals now
+  use the delete-a-block jackknife standard error; the plain sample standard
+  deviation previously used is `(K-1)/sqrt(K)` times too small (1.79x for
+  K=5), over-declaring significance.
+- PLS and PCA NIPALS now actually warn when the iteration cap is reached (the
+  PLS warning condition could never fire; PCA had no warning at all).
+- `target_projection` and `selectivity_ratio` used the raw-units regression
+  vector as a direction in the model's internal (scaled) space; with
+  `scale=True` (the default) the direction was wrong whenever the X columns
+  have different raw scales. The projection now uses the scaled-space beta and
+  maps X through the model's own scaler.
+- `TPLS.diagnose` poisoned an observation's scores, T2, SPE and predictions
+  with NaN if a single F or Z cell was missing (fit() handled the same case
+  correctly).
+- Hotelling's T2 in `fit()` skips components with ~zero score variance instead
+  of dividing by ~zero (a rank-deficient fit produced inf/NaN T2 for every
+  observation); `n_components < 1` is now rejected.
+- `MCUVScaler` treats a column with fewer than two observed values (NaN
+  standard deviation) as constant instead of emitting an all-NaN column;
+  `center()`/`scale()` with `axis=1` no longer broadcast the row statistic
+  across columns (a `ValueError` for rectangular input, silently wrong for
+  square input).
+- Generalized ESD outlier test (`detect_outliers_esd`): the number of outliers
+  is the LARGEST i with `R_i > lambda_i` (NIST/Rosner); the first crossing was
+  used, under-reporting in exactly the masking scenarios the test exists for.
+- Robust confidence interval for the median (metrics and the agent tool): the
+  missing `sqrt(pi/2)` factor on the median's standard error gave ~87 percent
+  coverage for a nominal 95 percent interval.
+- `variance_decomposition`: `between_stddev` now reports the between-group
+  variance component `sqrt((MS_between - MS_within)/n0)` instead of
+  `sqrt(MS_between)`, which mixed the within-group noise into the "between"
+  number (the docstring example itself showed the wrong value).
+- `biweight_midvariance` now uses the midvariance tuning constant c = 9; the
+  previous c = 6 is the biweight location constant and biased the scale low.
+- Holt-Winters control chart: the biweight rho function conflated the
+  consistency constant with the cutoff k = 2.52, making every derived scale
+  estimate 12 percent too small (nominal +/-3S limits were really +/-2.63
+  sigma, about 3x the false-alarm rate). Warm-up residuals now subtract the
+  fitted trend `beta_0 * t` rather than the constant `beta_0`. The lambda grid
+  search is NaN-aware (for 10 <= N < 20 every grid cell was NaN and (0.1, 0.1)
+  always won silently). An explicit `ld_1=0.0` is respected. Unknown chart
+  variants are rejected at construction with a clear message.
+- `calculate_cpk`: `rsd` is now the relative standard deviation of the data
+  itself; it previously divided by the distance-to-spec centre and changed
+  value when the specification moved. The capability tool reports an undefined
+  Cpk as "could not be computed" instead of "Poor capability".
+- Clear effects (`evaluate_design(metric="clear_effects")`) now follow Wu and
+  Hamada: an effect is clear only when every alias has order >= 3. The old
+  rule declared every main effect of a resolution-III design clear.
+- Explicit fractional-factorial generators: a generator on a non-last factor
+  (for example `"B=AC"`) silently swapped factor columns, and multi-character
+  factor names were misread as products of single letters. Generators are now
+  parsed against the real factor names and the columns mapped back to the
+  requested factor order; negative generators are supported and inconsistent
+  generator sets are rejected.
+- `Column.to_coded`/`to_realworld` no longer ignore an explicit `center=0`
+  (falsy); missing or zero-width ranges raise a clear error. `gather()` no
+  longer silently discards positional arguments.
+- D-optimal point exchange: the scorer no longer de-duplicates the design
+  before computing `|X'X|` (replicated runs carry information), and an
+  improving swap onto the row with index label 0 is no longer discarded;
+  `point_exchange` accepts `random_state`.
+- The lack-of-fit test can now find replicates on generated designs: it groups
+  on the model's factor columns (with rounded numeric values) instead of the
+  whole frame, whose unique-per-row `RunOrder` column made every group a
+  singleton.
+- Robust regression: the degenerate-x guard branch itself divided by the
+  ~zero x sum-of-squares; the leverage there is exactly 1/N.
+- Batch DTW: the reported alignment distance summed the cumulative cost matrix
+  entries along the warping path (not a distance); it is now the accumulated
+  cost `D[-1, -1]`. Kassidas alignment weights now up-weight (rather than
+  effectively zero out) variables whose trajectories align near-perfectly.
+- Residual diagnostics (`analyze_experiment`): p-values that underflow to
+  exactly 0.0 (the most significant possible result) are no longer rendered as
+  "not available".
+- `Settings` (config) now genuinely caches on first access (the `setdefault`
+  pattern re-read the environment on every access and re-raised later if the
+  env var went bad after a successful read); numeric knobs must be positive.
+- `clean()` now serialises `numpy.bool_`, numpy-keyed dicts, sets and the
+  remaining numpy scalar types, closing "internal error" failures at the MCP
+  boundary.
+- `discover_tools` no longer swallows a missing FIRST-PARTY module as an
+  optional dependency: a typo'd or renamed `process_improve` module now
+  propagates instead of silently dropping a whole tool category.
+- `safe_execute_tool_call` runs each default-path call in a private worker
+  pool. The shared module pool was not thread-safe: with concurrent calls
+  (the MCP server runs tools on executor threads) one thread's teardown could
+  kill the worker running another thread's task, mis-diagnosed as a memory
+  limit kill, or leak an orphaned worker.
+- `raincloud` without the `plotting` extra now raises the documented
+  "install the extra" `ImportError` at the call site instead of an
+  `AttributeError` from the module stub.
+
+### Changed
+
+- `detect_outliers_esd` and the `detect_outliers` agent tool now default to
+  the classical (mean/std) statistic, `robust_variant=False`. The MAD-scaled
+  variant is tested against critical values derived for the classical
+  statistic and declares outliers in clean data; it stays available as an
+  explicitly documented screening heuristic. This changes results for callers
+  who relied on the old default.
+- `PROCESS_IMPROVE_*` boolean environment variables now reject unrecognized
+  values with a `ValueError` instead of silently reading them as false; a typo
+  in `PROCESS_IMPROVE_MCP_SAFE_MODE` previously disabled safe mode with no
+  error.
+- The agent-facing `control_chart` tool no longer advertises a `cusum` chart
+  type: it never existed and always failed with a misleading error. The
+  package docstring's "CUSUM, EWMA" claim is corrected likewise.
+- CI: the no-op `create` trigger is removed from the test workflow, and both
+  workflows now grant least-privilege token scopes (the docs build job no
+  longer carries Pages-deploy credentials while executing PR code).
+
 ## [1.67.1] - 2026-08-22
 
 ### Fixed
@@ -3150,7 +3289,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.67.1...HEAD
+[Unreleased]: https://github.com/kgdunn/process-improve/compare/v1.68.0...HEAD
+[1.68.0]: https://github.com/kgdunn/process-improve/compare/v1.67.1...v1.68.0
 [1.67.1]: https://github.com/kgdunn/process-improve/compare/v1.67.0...v1.67.1
 [1.67.0]: https://github.com/kgdunn/process-improve/compare/v1.66.1...v1.67.0
 [1.66.1]: https://github.com/kgdunn/process-improve/compare/v1.66.0...v1.66.1

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.67.1
+version: 1.68.0
 date-released: "2026-08-22"
 keywords:
   - chemometrics

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "process-improve"
-version = "1.67.1"
+version = "1.68.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/__init__.py

@@ -13,7 +13,8 @@
 experiments
     Factorial and response surface experiment designs.
 monitoring
-    Control charts (Shewhart, CUSUM, EWMA).
+    Control charts (Shewhart, and a robust Holt-Winters chart that blends
+    Shewhart and CUSUM-like behaviour) and process capability metrics.
 batch
     Batch process data analysis.
 regression

src/process_improve/batch/alignment_helpers.py

@@ -58,38 +58,40 @@ def distance_matrix(test: np.ndarray, ref: np.ndarray, weight_matrix: np.ndarray
 
 @jit(nopython=True)
 def backtrack_optimal_path(D: np.ndarray) -> tuple[np.ndarray, float]:
-    """Backtrack through the distance matrix to find the optimal warping path."""
+    """Backtrack through the distance matrix to find the optimal warping path.
+
+    Returns the path and the DTW distance, which is the ACCUMULATED cost at
+    the end of the alignment, ``D[-1, -1]``. An earlier version summed the
+    cumulative ``D`` entries along the path - a sum of prefix sums that grows
+    super-linearly with path length and is not a distance; the per-batch
+    alignment-quality numbers built from it were meaningless.
+    """
     nr, nt = D.shape
     nr -= 1
     nt -= 1
-    path_sum = 0.0
+    distance = float(D[nr, nt])
     path = [
         [nr, nt],
     ]
     while (nt + nr) != 0:
         if nt == 0:
             nr -= 1
-            path_sum = path_sum + float(D[nr, nt])
         elif nr == 0:
             nt -= 1
-            path_sum = path_sum + float(D[nr, nt])
         else:
             # Commented-code here is to read, but for Numba JIT, the other code is able to be
             # compiled. They give the same results in regular Python.
             # number = np.argmin([D[nr - 1, nt - 1], D[nr, nt - 1], D[nr - 1, nt]])
             a, b, c = D[nr - 1, nt - 1], D[nr, nt - 1], D[nr - 1, nt]
             if (a <= b) & (a <= c):
                 # assert number == 0
-                path_sum = path_sum + D[nr - 1, nt - 1]
                 nt -= 1
                 nr -= 1
             elif (b <= a) & (b <= c):
                 # assert number == 1
-                path_sum = path_sum + D[nr, nt - 1]
                 nt -= 1
             elif (c <= a) & (c <= b):
                 # assert number == 2
-                path_sum = path_sum + D[nr - 1, nt]
                 nr -= 1
             else:
                 raise AssertionError
@@ -98,4 +100,4 @@ def backtrack_optimal_path(D: np.ndarray) -> tuple[np.ndarray, float]:
 
     # All done:
     path.reverse()
-    return np.array(path), path_sum
+    return np.array(path), distance

src/process_improve/batch/preprocessing.py

@@ -402,7 +402,16 @@ def batch_dtw(  # noqa: C901, PLR0915
             # TODO: make this a configurable setting
             # problematic_threshold = dist_df["Distance"].quantile(0.95)
 
-        next_weights = 1.0 / np.where(next_weights > epsqrt, next_weights, 10000)
+        # Kassidas: each variable's weight is inversely proportional to its
+        # summed squared deviation from the average trajectory, so a variable
+        # that aligns consistently (SSQ ~ 0) must get a LARGE weight. The
+        # previous guard substituted the magic value 10000 for a near-zero
+        # SSQ, handing the best-aligned variables a weight of ~1e-4 - the
+        # exact opposite - and the constant was scale-dependent. Floor the
+        # SSQ (relative to the largest observed SSQ) instead, so the weight
+        # stays large but finite.
+        ssq_floor = max(epsqrt, 1e-6 * float(np.max(next_weights)))
+        next_weights = 1.0 / np.maximum(next_weights, ssq_floor)
         weight_vector = (next_weights / np.sum(next_weights) * len(columns_to_align)).ravel()
         # If change in delta_weight is small, we terminate early; no need to fine-tune excessively.
         delta_weight = np.diag(weight_matrix) - weight_vector  # old - new