feat: surface tsam's cross-column concurrency metrics

docs/data-model.md

@@ -36,6 +36,7 @@ graph LR
     D --- F5[".segment_durations<br/><i>cluster, timestep, *slice_dims | None</i>"]
 
     Meta --- A[".accuracy<br/><b>→ AccuracyMetrics</b>"]
+    Meta --- CC[".concurrency<br/><b>→ ConcurrencyMetrics | None</b>"]
     Meta --- C[".clustering<br/><b>→ ClusteringResult</b>"]
 ```
 
@@ -121,6 +122,22 @@ Per-column metrics as DataArrays, plus weighted scalars.
 | `weighted_mae` | float | Scalar MAE weighted by column weights |
 | `weighted_rmse_duration` | float | Scalar duration RMSE weighted by column weights |
 
+## ConcurrencyMetrics
+
+Scalars measuring how well the joint structure *across* the clustered columns
+— which values co-occur in time — survives aggregation, where
+`AccuracyMetrics` measures each column on its own. Lower is better; both are
+`NaN` when a single column is clustered.
+
+`result.concurrency` is `None` on tsam < 4, which does not compute these.
+
+| Field | Type | Description |
+|-------|------|-------------|
+| `correlation_error` | float | Frobenius norm of the difference between the Pearson correlation matrices of the original and the reconstructed columns |
+| `rank_correlation_error` | float | The same for Spearman rank correlation, a copula proxy invariant to monotone changes in the marginals |
+
+With slice dims, both are DataArrays over `(*slice_dims)` rather than scalars.
+
 ## Glossary
 
 | Term | Meaning |

src/tsam_xarray/__init__.py

@@ -3,7 +3,11 @@
 from tsam_xarray._clustering import ClusteringInfo, ClusteringResult
 from tsam_xarray._core import aggregate
 from tsam_xarray._dim_names import DimNames
-from tsam_xarray._result import AccuracyMetrics, AggregationResult
+from tsam_xarray._result import (
+    AccuracyMetrics,
+    AggregationResult,
+    ConcurrencyMetrics,
+)
 from tsam_xarray._tuning import (
     TuningResult,
     find_best_combination,
@@ -19,6 +23,7 @@
     "AggregationResult",
     "ClusteringInfo",
     "ClusteringResult",
+    "ConcurrencyMetrics",
     "DimNames",
     "TuningResult",
     "aggregate",

src/tsam_xarray/_clustering.py

@@ -641,6 +641,7 @@ def _apply_single(
 
     from tsam_xarray._core import (
         _cluster_counts,
+        _concurrency_metrics,
         _metric_to_da,
         _reconstructed_to_da,
         _representatives_to_da,
@@ -707,6 +708,7 @@ def _make_accuracy() -> AccuracyMetrics:
         segment_durations=seg_durations,
         _accuracy_factory=_make_accuracy,
         _reconstructed_factory=_make_reconstructed,
+        _concurrency_factory=lambda: _concurrency_metrics(tsam_result),
         original=da,
         clustering=clustering_info,
         is_transferred=True,

src/tsam_xarray/_core.py

@@ -13,12 +13,26 @@
 import xarray as xr
 
 from tsam_xarray._dim_names import DimNames
-from tsam_xarray._result import AccuracyMetrics, AggregationResult
+from tsam_xarray._result import AccuracyMetrics, AggregationResult, ConcurrencyMetrics
 
 Weights = dict[str, float] | dict[str, dict[str, float]] | None
 ClusterOn = str | Sequence[str] | dict[str, Sequence[str]] | None
 
 
+def _concurrency_metrics(tsam_result: Any) -> ConcurrencyMetrics | None:
+    """Cross-column concurrency metrics, or None on a tsam that lacks them.
+
+    tsam v4 added ``AggregationResult.concurrency``; v3 has no equivalent.
+    """
+    concurrency = getattr(tsam_result, "concurrency", None)
+    if concurrency is None:
+        return None
+    return ConcurrencyMetrics(
+        correlation_error=xr.DataArray(concurrency.correlation_error),
+        rank_correlation_error=xr.DataArray(concurrency.rank_correlation_error),
+    )
+
+
 def _cluster_counts(tsam_result: Any) -> dict[int, float]:
     """Per-cluster occurrence counts from a tsam result, across tsam versions.
 
@@ -778,6 +792,7 @@ def _make_accuracy() -> AccuracyMetrics:
         segment_durations=seg_durations,
         _accuracy_factory=_make_accuracy,
         _reconstructed_factory=_make_reconstructed,
+        _concurrency_factory=lambda: _concurrency_metrics(tsam_result),
         original=da,
         clustering=clustering_info,
     )
@@ -819,6 +834,25 @@ def _recursive_concat(node: Any, dims: list[str]) -> xr.DataArray:
     return _recursive_concat(nested, slice_dims)
 
 
+def _concat_concurrency(
+    results: list[AggregationResult],
+    slice_dims: list[str],
+    slice_coords: dict[str, Any],
+) -> ConcurrencyMetrics | None:
+    """Concatenate per-slice concurrency metrics, or None if any slice lacks them."""
+    metrics = [m for m in (r.concurrency for r in results) if m is not None]
+    if len(metrics) != len(results):
+        return None
+    return ConcurrencyMetrics(
+        correlation_error=_concat_along_dims(
+            [m.correlation_error for m in metrics], slice_dims, slice_coords
+        ),
+        rank_correlation_error=_concat_along_dims(
+            [m.rank_correlation_error for m in metrics], slice_dims, slice_coords
+        ),
+    )
+
+
 def _concat_results(
     results: list[AggregationResult],
     slice_dims: list[str],
@@ -884,6 +918,9 @@ def _acc_field(field_name: str) -> xr.DataArray:
             ),
         ),
         _reconstructed_factory=lambda: _field("reconstructed"),
+        _concurrency_factory=lambda: _concat_concurrency(
+            results, slice_dims, slice_coords
+        ),
         original=_field("original"),
         clustering=merged_clustering,
         is_transferred=first.is_transferred,

src/tsam_xarray/_result.py

@@ -17,6 +17,13 @@
     from tsam_xarray._dim_names import DimNames
 
 
+def _fmt_metric(da: xr.DataArray) -> str:
+    mean = float(da.mean())
+    if da.size <= 1:
+        return f"{mean:.4f}"
+    return f"{mean:.4f} [{float(da.min()):.4f}-{float(da.max()):.4f}]"
+
+
 @dataclass(frozen=True, repr=False)
 class AccuracyMetrics:
     """Accuracy metrics from time series aggregation.
@@ -45,18 +52,44 @@ class AccuracyMetrics:
     weighted_rmse_duration: xr.DataArray
 
     def __repr__(self) -> str:
-        def _fmt(da: xr.DataArray) -> str:
-            mean = float(da.mean())
-            if da.size <= 1:
-                return f"{mean:.4f}"
-            return f"{mean:.4f} [{float(da.min()):.4f}-{float(da.max()):.4f}]"
-
         return (
             f"AccuracyMetrics("
-            f"weighted_rmse={_fmt(self.weighted_rmse)}, "
-            f"weighted_mae={_fmt(self.weighted_mae)}, "
+            f"weighted_rmse={_fmt_metric(self.weighted_rmse)}, "
+            f"weighted_mae={_fmt_metric(self.weighted_mae)}, "
             f"weighted_rmse_duration="
-            f"{_fmt(self.weighted_rmse_duration)})"
+            f"{_fmt_metric(self.weighted_rmse_duration)})"
+        )
+
+
+@dataclass(frozen=True, repr=False)
+class ConcurrencyMetrics:
+    """Cross-column concurrency metrics from time series aggregation.
+
+    Measures how well the joint structure across the clustered columns --
+    which values co-occur in time -- survives aggregation, complementing the
+    per-column error in `AccuracyMetrics`. Lower is better; both values are
+    ``NaN`` for a single clustered column.
+
+    Requires tsam >= 4. See `AggregationResult.concurrency`.
+
+    Attributes:
+        correlation_error: Frobenius norm of the difference between the
+            Pearson correlation matrices of the original and the
+            reconstructed columns. Dims: ``(*slice_dims)`` or scalar.
+        rank_correlation_error: The same for the Spearman rank-correlation
+            matrices, a copula proxy invariant to monotone changes in the
+            marginals. Dims: ``(*slice_dims)`` or scalar.
+    """
+
+    correlation_error: xr.DataArray
+    rank_correlation_error: xr.DataArray
+
+    def __repr__(self) -> str:
+        return (
+            f"ConcurrencyMetrics("
+            f"correlation_error={_fmt_metric(self.correlation_error)}, "
+            f"rank_correlation_error="
+            f"{_fmt_metric(self.rank_correlation_error)})"
         )
 
 
@@ -81,6 +114,9 @@ class AggregationResult:
             Computed on first access; on a tsam that defers
             metric computation (v4), never reading it skips
             the computation entirely.
+        concurrency: Cross-column concurrency metrics, or
+            ``None`` on tsam < 4, which does not compute them.
+            Computed on first access, like ``accuracy``.
         reconstructed: Reconstructed time series
             (same shape and dim order as ``original``).
             Computed on first access, like ``accuracy``.
@@ -104,6 +140,9 @@ class AggregationResult:
     _reconstructed_factory: Callable[[], xr.DataArray] = field(
         kw_only=True, repr=False, compare=False
     )
+    _concurrency_factory: Callable[[], ConcurrencyMetrics | None] = field(
+        kw_only=True, repr=False, compare=False
+    )
 
     @cached_property
     def accuracy(self) -> AccuracyMetrics:
@@ -115,6 +154,15 @@ def reconstructed(self) -> xr.DataArray:
         """Reconstructed series on the original time axis, computed on first access."""
         return self._reconstructed_factory()
 
+    @cached_property
+    def concurrency(self) -> ConcurrencyMetrics | None:
+        """Cross-column concurrency metrics, computed on first access.
+
+        ``None`` on tsam < 4, which does not compute them. See
+        `ConcurrencyMetrics`.
+        """
+        return self._concurrency_factory()
+
     def __repr__(self) -> str:
         c = self.clustering
         slices = f", slice_dims={c.slice_dims}" if c.slice_dims else ""

test/test_concurrency.py

@@ -0,0 +1,123 @@
+"""AggregationResult.concurrency, tsam's cross-column concurrency metrics."""
+
+from __future__ import annotations
+
+import numpy as np
+import pandas as pd
+import pytest
+import tsam
+import xarray as xr
+
+import tsam_xarray
+from tsam_xarray import ConcurrencyMetrics, aggregate
+
+HAS_CONCURRENCY = hasattr(tsam.AggregationResult, "concurrency")
+
+requires_concurrency = pytest.mark.skipif(
+    not HAS_CONCURRENCY,
+    reason="tsam < 4 does not compute concurrency metrics",
+)
+
+
+def _data(n_slices: int = 1, variables: list[str] | None = None) -> xr.DataArray:
+    if variables is None:
+        variables = ["a", "b", "c"]
+    rng = np.random.default_rng(0)
+    n_t = 14 * 24
+    dims = ["variable", "time"]
+    shape: tuple[int, ...] = (len(variables), n_t)
+    coords: dict[str, object] = {
+        "variable": variables,
+        "time": pd.date_range("2023-01-01", periods=n_t, freq="h"),
+    }
+    if n_slices > 1:
+        dims = ["scenario", *dims]
+        shape = (n_slices, *shape)
+        coords["scenario"] = [f"s{i}" for i in range(n_slices)]
+    return xr.DataArray(rng.random(shape), dims=dims, coords=coords, name="load")
+
+
+def _aggregate(da: xr.DataArray, n_clusters: int = 4):
+    return aggregate(da, time_dim="time", cluster_dim="variable", n_clusters=n_clusters)
+
+
+def test_concurrency_metrics_is_exported():
+    assert tsam_xarray.ConcurrencyMetrics is ConcurrencyMetrics
+
+
+@pytest.mark.skipif(HAS_CONCURRENCY, reason="tsam >= 4 computes concurrency metrics")
+def test_none_without_tsam_support():
+    assert _aggregate(_data()).concurrency is None
+
+
+@requires_concurrency
+def test_scalar_without_slice_dims():
+    concurrency = _aggregate(_data()).concurrency
+
+    assert isinstance(concurrency, ConcurrencyMetrics)
+    for metric in (concurrency.correlation_error, concurrency.rank_correlation_error):
+        assert metric.dims == ()
+        assert float(metric) >= 0
+
+
+@requires_concurrency
+def test_dims_follow_slice_dims():
+    da = _data(n_slices=3)
+    concurrency = _aggregate(da).concurrency
+
+    for metric in (concurrency.correlation_error, concurrency.rank_correlation_error):
+        assert metric.dims == ("scenario",)
+        xr.testing.assert_identical(metric.coords["scenario"], da.coords["scenario"])
+
+
+@requires_concurrency
+def test_exact_reconstruction_has_no_concurrency_error():
+    da = _data()
+    n_periods = da.sizes["time"] // 24
+    concurrency = _aggregate(da, n_clusters=n_periods).concurrency
+
+    assert float(concurrency.correlation_error) == pytest.approx(0, abs=1e-9)
+    assert float(concurrency.rank_correlation_error) == pytest.approx(0, abs=1e-9)
+
+
+@requires_concurrency
+def test_nan_for_a_single_clustered_column():
+    concurrency = _aggregate(_data(variables=["a"])).concurrency
+
+    assert np.isnan(float(concurrency.correlation_error))
+    assert np.isnan(float(concurrency.rank_correlation_error))
+
+
+@requires_concurrency
+@pytest.mark.parametrize("n_slices", [1, 3])
+def test_deferred_until_accessed(n_slices):
+    result = _aggregate(_data(n_slices))
+
+    assert "concurrency" not in result.__dict__
+
+    concurrency = result.concurrency
+    assert "concurrency" in result.__dict__
+    assert result.concurrency is concurrency
+
+
+@requires_concurrency
+@pytest.mark.parametrize("n_slices", [1, 3])
+def test_apply_reports_concurrency(n_slices):
+    da = _data(n_slices)
+    result = _aggregate(da)
+    transferred = result.clustering.apply(da)
+
+    assert transferred.is_transferred
+    xr.testing.assert_allclose(
+        transferred.concurrency.correlation_error,
+        result.concurrency.correlation_error,
+    )
+
+
+@requires_concurrency
+def test_repr_reports_both_metrics():
+    text = repr(_aggregate(_data()).concurrency)
+
+    assert text.startswith("ConcurrencyMetrics(")
+    assert "correlation_error=" in text
+    assert "rank_correlation_error=" in text