v1.4.1: engine additions for probabilistic forecasting

CHANGELOG.md

@@ -15,11 +15,27 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 > **[0.0.6]** corresponds to the legacy `weakl` package available on PyPI.
 
 ---
-## [Unreleased]
+## [1.4.1] - 2026-10-08
+
+Patch release for probabilistic forecasting: new methods and options only, plus one fix. Existing calls give the same forecasts; the only change of behaviour is that a distributional or mixture `StaticTAM` with the default log target now raises on a training target with values `<= 0` instead of returning a meaningless model (see "Fixed").
 
 - The `DCO` sign-off is required on every commit of a pull request (`DCO`, `CONTRIBUTING.md`, `docs/contributing/dco_howto.md`); a `REUSE` workflow runs `reuse lint` on each push and pull request.
 - The `pid()` page of `architecture/core/06_the_spectrum_api.md` says when `d_pen` acts and that the derivative is computed along the frame given to `predict()`; three tests pin `pid` behaviour.
 
+### Added
+- **Limiting the CPU threads is documented and tested** (`architecture/core/04_hardware_memory.md`, README): set `OMP_NUM_THREADS` and `MKL_NUM_THREADS` before starting Python. The page also describes a hang of `torch.set_num_threads(n)` seen on one PyTorch / oneMKL build, which reproduces without `tam`; the results of `tam` are unchanged.
+- **`loss="negative_binomial"` and `loss="tweedie"`** on `StaticTAM` (log link): overdispersed counts, and positive targets with exact zeros (sparse retail, rain, night-time solar). The negative binomial dispersion is estimated by maximum likelihood after the IRLS fit (or fixed with `loss_kwargs={"dispersion": alpha}`); the Tweedie power is `loss_kwargs={"power": p}` with `1 < p < 2` (1.5 by default) and its dispersion is the Pearson estimate. The estimate is `model.dispersion_`. A negative target raises a `ValueError`.
+- **`StaticTAM.predict_quantiles(data, taus)` for `poisson`, `negative_binomial` and `tweedie`**: discrete quantiles for the counts and the quantiles of the compound Poisson-gamma law for Tweedie, which are `0` where the probability of zero exceeds the level; non-negative and non-decreasing in the level by construction. The Poisson loss is otherwise unchanged.
+- **`KalmanTAM.predict_quantiles(data, taus)`**: Gaussian quantiles of the online forecast from the predictive variance of the filter, `x' P x + R` in the target scale (for `horizon_steps > 1`, the covariance of the delayed state grown by the process noise of the delay). Returns a DataFrame with `q<tau>` columns, like the distributional `StaticTAM`; like `predict_online`, `data` holds the target. The constructor, `predict_online`, `predict` and `fit` are unchanged. The intervals are as calibrated as `observation_noise_var` and `process_noise_var` are.
+- **`AdaptiveTAM.predict_quantiles(data, taus)`**: online quantile forecasts when the `base_model` is a distributional `StaticTAM` (a dict formula). The adaptive correction moves the location; the scale of the base is rescaled at each update by the root mean square of the standardized residuals of the last `training_window_periods` rows that are at least `horizon_steps` old (past residuals only), and the quantiles go back through the target transform of the base. Returns a DataFrame with `q<tau>` columns, like the distributional `StaticTAM`; like `predict_online`, `data` holds the target. A plain base raises a `ValueError`. `predict_online`, `predict` and `fit` give the same results as before for a plain base, and now also accept a distributional one.
+- **`dist_kwargs={"target_transform": "asinh" | "logit"}`** on a distributional `StaticTAM`: the location-scale law sits on `arsinh((y - c) / s)` (centre and robust scale estimated on the training target; a skewed, heavy-tailed target with negative values, such as electricity prices) or on the logit of `(y - a) / (b - a)` (a target bounded in `bounds=(a, b)`, estimated from the training range when not given). `predict_median`, `predict_quantile(s)`, `cdf`, `crps`, `anomaly_score` and the conformal intervals work through the transform, and quantiles stay inside the bounds for the logit. `target_transform="log"` is the default and the same as `log_target=True`; `"none"` the same as `log_target=False`. A mixture accepts `"none"` and `"log"` only.
+- **`f(x, ..., period=...)`**: the range of one cycle of a Fourier input, on which the feature is normalised in `StaticTAM` and in each rolling window of `AdaptiveTAM`. `period=None` (default): the minimum and maximum of the training rows, as before. `period='auto'`: first value to last value plus one step, the step read from the training rows (hours 0 to 23 give (0, 24), months 1 to 12 give (1, 13)). `period=(low, high)`: given, for a history covering only part of the cycle (`(0, 1)` for a time of year trained on half a year). Without it the training range is still taken as one period, which glues the last value of an integer-coded position onto the first (`f(hour, cyclic=True)` with hours 0 to 23 forces 23:00 = 0:00; Sunday = Monday for 0 to 6; December = January for 1 to 12) and reads a history covering half a year (or a 60-day adaptive window) as a whole year: use `'auto'` for the first case and a given range for the second. Results without `period` are unchanged.
+- **`OperaTAM(..., loss_type='pinball', tau=0.9)`**: online aggregation of quantile forecasts with the pinball loss of level `tau` (MLpol and EWA). `tam.model.opera.aggregate_quantiles(df, target_col, expert_cols_by_level)` aggregates every level and sorts them on each row. The square and absolute losses are unchanged.
+- **`SafetyTAM.predict_quantiles(y_pred, levels, scale=None, signed=False)`**: conformal quantiles at several levels in one call (`q<level>` columns, sorted on each row). By default the symmetric bands of `predict_intervals`; `signed=True` calibrates each tail on the signed residuals kept by `calibrate()` (one-sided split conformal), for skewed errors. `calibrate(..., verbose=False)` silences the sample-count line, printed by default as before.
+
+### Fixed
+- **A distributional or mixture `StaticTAM` with the default log target rejects a training target with values `<= 0`** (`ValueError` naming the count and `dist_kwargs={"log_target": False}`). They were clipped to a tiny constant before the log, i.e. silently turned into huge negative outliers, and the fitted location, scale and quantiles were meaningless. Positive targets are unchanged.
+
 ## [1.4.0] - 2026-10-06
 
 First release on PyPI and Zenodo since 1.3.1; it also carries everything of 1.3.2 to 1.3.4, which were tags only. **`AdaptiveTAM`, `KalmanTAM` and `NeuralTAM` [EXP] forecasts change** (see 1.3.3, 1.3.4 and the sections below); `StaticTAM` changes only for `te()` with GCV (-2%), a mixture mean (-0.6%) and a prediction without the group column. To get the 1.3.1 behaviour: `pip install tam-ml==1.3.1`.
@@ -308,7 +324,8 @@ This version introduced the Formula API and the first object-oriented refactorin
 ### Added
 * Initial project setup based on the original `weakl` v0.0.6 package.
 
-[Unreleased]: https://github.com/EDF-Lab/tam/compare/v1.4.0...HEAD
+[Unreleased]: https://github.com/EDF-Lab/tam/compare/v1.4.1...HEAD
+[1.4.1]: https://github.com/EDF-Lab/tam/releases/tag/v1.4.1
 [1.4.0]: https://github.com/EDF-Lab/tam/releases/tag/v1.4.0
 [1.3.4]: https://github.com/EDF-Lab/tam/releases/tag/v1.3.4
 [1.3.3]: https://github.com/EDF-Lab/tam/releases/tag/v1.3.3

CITATION.cff

@@ -21,8 +21,8 @@ identifiers:
 
 
 repository-code: "https://github.com/EDF-Lab/tam"
-version: 1.4.0
-date-released: 2026-10-06
+version: 1.4.1
+date-released: 2026-10-08
 
 keywords:
   - "Python"

README.md

@@ -2,7 +2,7 @@
 
 [🧪 `EXAMPLES`](EXAMPLES.md) | [📚 `THEORY`](THEORY.md) | [📄 `PAPER`](paper.md)
 
-**Version:** 1.4.0  
+**Version:** 1.4.1  
 **License:** LGPL-3.0-or-later.  
 **Copyright (c):** 2023-2026 EDF (Electricité De France).  
 **Copyright (c):** 2023-2025 Sorbonne Université.  
@@ -79,6 +79,8 @@ Raw data flows through dynamic normalizers before being processed by the `Static
 pip install tam-ml
 ```
 
+> **Threads:** to limit the CPU threads, set `OMP_NUM_THREADS` and `MKL_NUM_THREADS` in the environment before starting Python (see `architecture/core/04_hardware_memory.md`).
+
 > ⚠️ **Hardware Note:** TAM relies heavily on PyTorch for VRAM management and tensor acceleration. Ensure you have installed the correct PyTorch distribution (CUDA for NVIDIA, or MPS for Apple Silicon) for your specific hardware to fully leverage the framework's GPU capabilities.
 
 ---
@@ -167,10 +169,10 @@ If you use these packages in your research, please cite them using their permane
 **TAM**
 ```bibtex
 @misc{tam2026package,
-  title={TAM: Time series Additive Model (v1.4.0)},
+  title={TAM: Time series Additive Model (v1.4.1)},
   author={Allioux, Yann and Goude, Yannig},
   year={2026},
-  doi={10.5281/zenodo.23171729},
+  doi={10.5281/zenodo.23248249},
   note={With foundational WeaKL research by Nathan Doumèche, and tensorisation and adaptive model contributions by Éloi Bedek}
 }
 ```

architecture/core/04_hardware_memory.md

@@ -24,6 +24,18 @@ Crucially, it manages disaster recovery via the `handle_oom` method. When an ope
 
 ---
 
+## Limiting the CPU threads
+
+To share a machine between several fits (a benchmark, an AutoTAM search run in parallel), limit the threads **through the environment, before Python starts**:
+
+```bash
+OMP_NUM_THREADS=2 MKL_NUM_THREADS=2 python my_script.py
+```
+
+On Windows PowerShell: `$env:OMP_NUM_THREADS=2; $env:MKL_NUM_THREADS=2; python my_script.py`. `torch.get_num_threads()` then returns 2 (`tests/test_thread_limits.py` runs a grouped fit this way in a fresh process).
+
+**Known issue on one build.** With the PyTorch 2.14.0 CPU wheel on Windows (Intel oneMKL 2026.1), calling `torch.set_num_threads(n)` with `n >= 2` made every later batched LU factorisation of matrices of a few hundred rows (`torch.linalg.solve`, `lu_factor`, `inv`, which the grouped solves use) print `Intel oneMKL ERROR: Parameter 6 was incorrect on entry to DLASWP` and hang. It reproduces with PyTorch alone, without `tam`, so it is not specific to this library (Cholesky factorisations, unbatched solves and `n = 1` were not affected). Other platforms and builds were not tested. If a fit prints that message on your setup, remove the `torch.set_num_threads` call and use the environment variables above.
+
 ## The Memory Oracle and Safe Chunking
 
 To proactively avoid invoking the OOM handler, the framework utilizes `_memory.py` as an advanced predictive oracle. 

architecture/core/06_the_spectrum_api.md

@@ -65,6 +65,19 @@ The linear effect applies direct spatial scaling, and its penalty is a standard
 :end-before: "#: </parse_fourier>"
 ```
 
+**Declared period.** `f(x, ..., period=...)` makes the feature normalised on the range of one cycle instead of the minimum and
+maximum of the training rows: `None` (default) keeps the training range, `'auto'` uses (first value, last value + one step), the
+step being the smallest gap between the distinct training values (resolved once on all training rows by `_resolve_fixed_ranges`,
+so every group and adaptive window shares it), and `(low, high)` is used as given (a history covering part of the cycle). This applies in `StaticTAM` and in every rolling window
+of `AdaptiveTAM` (`KalmanTAM` fits through a `StaticTAM`). Without it the training range is the period, as before: an integer-coded position
+has its last value glued onto its first (hours 0 to 23 give 23:00 = 0:00; use `period='auto'` or `(0, 24)`), and a history covering half a cycle is read as a whole cycle. `fixed_ranges` collects the declared periods (tensor products included) and
+`_fit_normalization_params` / `_transform_data_adaptive` apply them; two different periods for one feature raise. The default
+`extrapolate='continue'` lets an input beyond the period wrap around the cycle (`'constant'` would clamp it).
+```{literalinclude} ../../../../src/tam/model/spectrum/_factory.py
+:language: python
+:pyobject: fixed_ranges
+```
+
 **2. Feature Map**
 Evaluates batched trigonometric functions, dynamically scaling angular frequencies based on the `cyclic` boundary flag to prevent endpoint distortions.
 ```{literalinclude} ../../../../src/tam/model/spectrum/_fourier.py

architecture/meta/01_adaptive_code.md

@@ -96,6 +96,28 @@ The `grid_search_fit()` method implements the "Multi-Start" hyperparameter solve
 :caption: src/tam/model/adaptative.py (Coordinate Descent Algorithm for Hyperparameters)
 ```
 
+## Online Quantiles from a Distributional Base
+
+`AdaptiveTAM(base_model=...)` also accepts a fitted distributional `StaticTAM` (a dict formula). The correction is then fitted on the model scale of its location sub-model (the target `__mu__`, the target after its transform): the constructor redirects the base to that sub-model, maps `Residual<target>` in `adaptive_formula` to `Residual__mu__`, and `_with_model_target` adds the model-scale column to any frame that holds the target. The forecasts of a plain base are untouched.
+
+`predict_quantiles(data, taus)` is the probabilistic counterpart of `predict_online`, so `data` holds the target. It runs the online simulation, then rescales the scale of the base with `_online_scale`, a rolling root mean square of the standardized residuals computed per group from cumulative sums (vectorised), and goes back to the response scale with the transform of the base:
+
+```{literalinclude} ../../../../src/tam/model/adaptative.py
+:language: python
+:start-after: "#: <online_scale>"
+:end-before: "#: </online_scale>"
+:caption: src/tam/model/adaptative.py (rolling scale of the standardized residuals)
+```
+
+```{literalinclude} ../../../../src/tam/model/adaptative.py
+:language: python
+:start-after: "#: <adaptive_quantiles>"
+:end-before: "#: </adaptive_quantiles>"
+:caption: src/tam/model/adaptative.py (online quantiles)
+```
+
+The first window forecasts the row `training_window_periods + horizon_steps - 1` of a group, and the scale is renewed every `update_interval_periods` rows; each renewal reads the residuals of the `training_window_periods` rows that are at least `horizon_steps` old, so no residual of a later row enters an earlier quantile. A renewal with fewer than five usable residuals keeps the scale of the base. A plain base raises a `ValueError`.
+
 ## Separation of Concerns: Simulation vs. Inference
 
 To guarantee speed and safety in operational production pipelines, `AdaptiveTAM` strictly separates the continuous historical simulation from out-of-sample inference.

architecture/meta/02_kalman_torchscript.md

@@ -18,3 +18,16 @@
 **Process noise per term (verified by `tests/test_kalman_per_term_noise.py`):** the TorchScript kernel takes the diagonal of $Q$ as a tensor (`process_noise_diag`, one variance per design column) instead of a scalar. `KalmanTAM._term_spans` maps each formula term (the names of `decomposition_names`, plus `"offset"` for the two constant columns: the filter's systematic offset and the formula's intercept) to its design columns, and `_noise_diagonal` fills the diagonal from `process_noise_var` (a float, or a dict `{term: q, "offset": q, "default": q}`). The boost on the offset's noise (`offset_q_boost`) is `offset_boost` unless the offset has a variance of its own; the initial covariance of the offset keeps `offset_boost`. A float call builds the same diagonal as before, so its forecasts are unchanged.
 
 **Base-model components (verified by `tests/test_base_effect_columns.py`):** with `add_base_effects=True`, one `l(effect_<name>)` term per base-model component is appended to the Kalman formula, with the column names `decompose_prediction` emits (`decomposition_names()`).
+
+## Predictive quantiles
+
+The compiled loop can also return the predictive variance of every forecast row (`with_variance`, off by default: the filter and its outputs are then identical to before). Per block it adds one batched product `X P Xᵀ` on the a priori covariance, and for `horizon_steps > 1` the covariance of the delayed state grown by `(h - 1) Q`, written at the row it forecasts. `_run_filter` keeps it in `_predictive_variance_`; `_online` is the former body of `predict_online`, which now calls it, and `predict_quantiles` converts the standard deviation to the target scale with the reference scale and adds `z_tau` standard deviations to the online forecast:
+
+```{literalinclude} ../../../../src/tam/model/kalman.py
+:language: python
+:start-after: "#: <kalman_quantiles>"
+:end-before: "#: </kalman_quantiles>"
+:caption: src/tam/model/kalman.py (Gaussian quantiles from the filter variance)
+```
+
+The constructor and the outputs of `predict_online`, `predict` and `fit` are unchanged. One limit of the filter carries over: a row without a target leaves the state unchanged but still reduces the covariance, so the variance after a long gap of missing targets is understated.

architecture/meta/05_opera_gpu.md

@@ -60,6 +60,16 @@ A `torch.where` mask safely intercepts any groups with perfectly zero targets to
 :caption: src/tam/model/opera.py (TorchScript Compiled EWA Loop)
 ```
 
+`loss_type='pinball'` with `tau` aggregates quantile forecasts: the MLpol loop uses the subgradient `1{y_hat > y} - tau` in place
+of the square or absolute one, the EWA loop the pinball loss of each expert. `aggregate_quantiles` runs one `OperaTAM` per level (each
+with its own expert columns) and sorts the levels on each row:
+
+```{literalinclude} ../../../../src/tam/model/opera.py
+:language: python
+:start-after: "#: <aggregate_quantiles>"
+:end-before: "#: </aggregate_quantiles>"
+```
+
 ---
 
 ## Causal Boundary Enforcement (Horizon Shifting)

architecture/statistics/01_reweighted_estimation_code.md

@@ -16,6 +16,28 @@ The `estimation/` package mirrors `spectrum/`: a base contract, one file per los
 :end-before: "#: </glm_family>"
 ```
 
+## Negative binomial and Tweedie
+
+`negative_binomial` and `tweedie` are `GLMFamily` subclasses in `_count_families.py`, registered in `build_strategy` (`StaticTAM(..., loss="negative_binomial", loss_kwargs={"dispersion": 0.5})` fixes the dispersion, `loss_kwargs={"power": 1.3}` sets the Tweedie power). The negative binomial variance reads the current `dispersion`, so the same IRLS loop serves any value of it:
+
+```{literalinclude} ../../../../src/tam/model/statistics/estimation/_count_families.py
+:language: python
+:start-after: "#: <negative_binomial>"
+:end-before: "#: </negative_binomial>"
+```
+
+After the first fit, `BaseTAM.fit` asks a family that has `estimate_dispersion` for the dispersion at the fitted means and, if it moved by more than 0.1%, refits the means with the new variance (at most ten rounds). The Tweedie dispersion does not enter the weights, so its estimate returns without a refit. `validate_target` rejects a negative target before the fit. The estimate is exposed as `StaticTAM.dispersion_` (None for the other losses).
+
+The Tweedie family computes the compound Poisson-gamma distribution function with the Poisson jump probabilities up to a tail of 1e-10, and its quantiles by bisection on that function (cost proportional to rows times jumps times 45 steps; a quantile below the mass at zero is exactly 0):
+
+```{literalinclude} ../../../../src/tam/model/statistics/estimation/_count_families.py
+:language: python
+:start-after: "#: <tweedie>"
+:end-before: "#: </tweedie>"
+```
+
+`StaticTAM.predict_quantiles(data, taus)` on a plain model whose loss has a `quantiles` method (`poisson`, `negative_binomial`, `tweedie`) returns the `q<tau>` columns of `predict_count_quantiles`; on any other plain loss it raises as before.
+
 ## Robust M-estimators and expectiles
 
 ```{literalinclude} ../../../../src/tam/model/statistics/estimation/_robust.py