SeasonalSimilarity¶

class SeasonalSimilarity(BaseSimilarity):
    r"""Temporal similarity using Fourier features for weighting observations.

    Computes observation weights by measuring the distance between
    cyclic temporal features extracted from prediction timestamps.
    Observations at similar seasonal positions (e.g. same day of week,
    same month of year) receive higher weights.

    Timestamps are converted to step indices relative to the first
    observed timestamp, then encoded as sin/cos pairs at the specified
    seasonal periods. Distances in this feature space are converted to
    weights using the same softmax formula as ``DistanceSimilarity``.

    Parameters
    ----------
    seasonality : list of float or None, default=None
        Seasonal periods in time steps (e.g. ``[7.0, 365.25]`` for
        weekly and yearly cycles on daily data). ``None`` is accepted at
        construction but invalid; passing ``None`` raises ``ValueError``
        at ``fit`` time, so a list must be provided before fitting.
    harmonics : dict mapping float to list of int, or None, default=None
        Harmonics to include per seasonality period. Keys must match
        entries in ``seasonality``. Each value is a list of positive
        integers specifying which harmonics to use. Defaults to
        ``{s: [1]}`` for each ``s`` in ``seasonality``.
    metric : str, default="euclidean"
        Distance metric for ``scipy.spatial.distance.cdist``.
    metric_params : dict or None, default=None
        Additional keyword arguments forwarded to the distance metric.

    Attributes
    ----------
    first_time_ : datetime
        Reference timestamp from the first calibration prediction.
    interval_td_ : timedelta
        Time interval between consecutive timestamps, auto-detected
        from calibration data. When ``fit`` receives a single timestamp the
        interval cannot be inferred and is set to ``timedelta(0)``; in that
        case ``_extract_features`` leaves the time axis in raw seconds, so
        ``seasonality`` must be expressed in seconds to remain meaningful.

    Notes
    -----
    Sin/cos encoding ensures that cyclic distances are correctly
    captured (e.g. December 31 is close to January 1). Multiple
    seasonalities combine naturally by concatenating feature vectors.

    The weight normalisation matches ``DistanceSimilarity`` exactly, via the
    shared [`BaseSimilarity._to_weights`][yohou.interval.base.BaseSimilarity]:
    a numerically-stable softmax of negative distances with uniform mass
    reserved for the test point, ``w_{ji} = raw_{ji} / (1 + \sum_k raw_{jk})``
    where ``raw_{ji} = \exp(-(d_{ji} - \max_k d_{jk}))``. Each row is
    non-negative and sums below 1, following the non-exchangeable conformal
    construction (Barber et al., 2023).

    References
    ----------
    [1] Barber, R.F., Candes, E.J., Ramdas, A., & Tibshirani, R.J.
        (2023). "Conformal prediction beyond exchangeability." Annals of
        Statistics, 51(2), 816-845.
        https://doi.org/10.1214/23-AOS2276

    See Also
    --------
    - [`DistanceSimilarity`][yohou.interval.similarity.DistanceSimilarity] : Value-based distance similarity.
    - [`BaseSimilarity`][yohou.interval.base.BaseSimilarity] : Abstract similarity base class.

    Examples
    --------
    >>> from datetime import datetime, timedelta
    >>> import polars as pl
    >>> import numpy as np
    >>> from yohou.interval.similarity import SeasonalSimilarity
    >>>
    >>> # Daily data with 3 weeks of calibration
    >>> dates = [datetime(2021, 1, 1) + timedelta(days=i) for i in range(21)]
    >>> y = pl.DataFrame({"time": dates, "value": np.random.randn(21)})
    >>> y_pred = pl.DataFrame({"time": dates, "value": np.random.randn(21)})
    >>>
    >>> # Fit with weekly seasonality
    >>> sim = SeasonalSimilarity(seasonality=[7.0])
    >>> _ = sim.fit(y, y_pred)
    >>>
    >>> # Predict weights for a new Monday
    >>> new_date = [datetime(2021, 1, 22)]
    >>> y_pred_new = pl.DataFrame({"time": new_date, "value": [0.5]})
    >>> weights = sim.predict(y_pred_new)
    >>> weights.shape
    (1, 21)

    """

    _parameter_constraints: dict = {
        "seasonality": [list],
        "harmonics": [dict, None],
        "metric": [str],
        "metric_params": [dict, None],
    }

    def __init__(
        self,
        seasonality: list[float] | None = None,
        harmonics: dict[float, list[int]] | None = None,
        metric: str = "euclidean",
        metric_params: dict[str, object] | None = None,
    ) -> None:
        self.seasonality = seasonality
        self.harmonics = harmonics
        self.metric = metric
        self.metric_params = metric_params

    def _resolve_harmonics(self) -> dict[float, list[int]]:
        """Resolve harmonics, defaulting to first harmonic per seasonality.

        Returns
        -------
        dict[float, list[int]]
            Mapping from seasonality period to list of harmonic indices.

        """
        if self.harmonics is not None:
            return self.harmonics
        return {s: [1] for s in self.seasonality}  # ty: ignore[not-iterable]

    def _extract_features(self, times: pl.Series) -> np.ndarray:
        """Extract Fourier features from a datetime series.

        Parameters
        ----------
        times : pl.Series
            Datetime series from which to compute features.

        Returns
        -------
        np.ndarray
            Feature matrix of shape ``(len(times), n_features)``.

        """
        time_diff = times - self.first_time_
        t = time_diff.dt.total_seconds().to_numpy().astype(np.float64)
        interval_seconds = self.interval_td_.total_seconds()
        if interval_seconds != 0:
            t = t / interval_seconds

        harmonics = self._resolve_harmonics()
        features = []
        for s in self.seasonality:  # ty: ignore[not-iterable]
            for k in harmonics.get(s, [1]):
                angle = 2.0 * math.pi * k * t / s
                features.append(np.sin(angle))
                features.append(np.cos(angle))

        return np.column_stack(features)

    @_fit_context(prefer_skip_nested_validation=True)
    def fit(
        self,
        y: pl.DataFrame,
        y_pred: pl.DataFrame,
        X_actual: pl.DataFrame | None = None,
    ) -> "SeasonalSimilarity":
        """Fit the temporal similarity from calibration predictions.

        Auto-detects the time interval from consecutive timestamps in
        ``y_pred`` and stores a reference timestamp and Fourier feature
        matrix for later distance computation.

        Parameters
        ----------
        y : pl.DataFrame
            Target time series (unused, accepted for API consistency).
        y_pred : pl.DataFrame
            Point forecast time series with a ``"time"`` column.
        X_actual : pl.DataFrame or None, default=None
            Exogenous features (unused, accepted for API consistency).

        Returns
        -------
        self

        Notes
        -----
        When ``y_pred`` contains only a single timestamp, the time interval
        cannot be inferred and ``interval_td_`` is set to ``timedelta(0)``;
        the feature extractor then works in raw seconds and ``seasonality``
        must also be expressed in seconds to remain meaningful.

        Raises
        ------
        ValueError
            If ``seasonality`` is ``None`` or empty.

        """
        if self.seasonality is None or len(self.seasonality) == 0:
            raise ValueError("seasonality must be a non-empty list of floats")

        times = y_pred["time"]
        self.first_time_ = times[0]

        if len(times) > 1:
            self.interval_td_ = times[1] - times[0]
        else:
            self.interval_td_ = timedelta(0)

        self._features_observed = self._extract_features(times)

        return self

    def observe(
        self,
        y: pl.DataFrame,
        y_pred: pl.DataFrame,
        X_actual: pl.DataFrame | None = None,
    ) -> "SeasonalSimilarity":
        """Observe new data and extend the reference feature matrix.

        Parameters
        ----------
        y : pl.DataFrame
            New target observations (unused, accepted for API
            consistency).
        y_pred : pl.DataFrame
            New predictions with a ``"time"`` column.
        X_actual : pl.DataFrame or None, default=None
            Exogenous features (unused, accepted for API consistency).

        Returns
        -------
        self

        """
        check_is_fitted(self, "first_time_")
        new_features = self._extract_features(y_pred["time"])
        self._features_observed = np.vstack([self._features_observed, new_features])
        return self

    def rewind(
        self,
        y: pl.DataFrame,
        y_pred: pl.DataFrame,
        X_actual: pl.DataFrame | None = None,
    ) -> "SeasonalSimilarity":
        """Rewind the most recently observed data.

        Removes the last ``len(y)`` rows from the internal feature
        matrix, reversing the effect of the corresponding ``observe()``
        call.

        Parameters
        ----------
        y : pl.DataFrame
            Target observations to rewind (used only for row count).
        y_pred : pl.DataFrame
            Predictions to rewind (used only for row count).
        X_actual : pl.DataFrame or None, default=None
            Exogenous features to rewind (unused).

        Returns
        -------
        self

        """
        check_is_fitted(self, "first_time_")
        n_rewind = len(y)
        self._features_observed = self._features_observed[: len(self._features_observed) - n_rewind]
        return self

    def predict(
        self,
        y_pred: pl.DataFrame,
        X_actual: pl.DataFrame | None = None,
    ) -> np.ndarray[tuple[int, int], np.dtype[np.floating[Any]]]:
        """Compute temporal similarity weights for new predictions.

        Parameters
        ----------
        y_pred : pl.DataFrame
            New predictions with a ``"time"`` column.
        X_actual : pl.DataFrame or None, default=None
            Exogenous features (unused).

        Returns
        -------
        np.ndarray
            Weight matrix of shape ``(n_predictions, n_calibration)``.

        """
        check_is_fitted(self, "first_time_")
        new_features = self._extract_features(y_pred["time"])

        distances: np.ndarray = cdist(
            new_features,
            self._features_observed,
            metric=self.metric,
            **(self.metric_params or {}),
        )  # ty: ignore[no-matching-overload]
        return self._to_weights(distances)