# (c) 2026 DTU Wind Energy
"""
Measure-Correlate-Predict MCP Module.
"""
import numpy as np
from ..wind import wd_to_sector
from .._errors import WindkitValidationError
from ._validation import assert_timesteps_equal, assert_valid_values, ltc_validate
from .regression import LinearRegression, VarianceRatioLinearRegression
class MCPBase:
"""Sectorwise Regressor base-class used as a blueprint
for sectorwise MCP methods based on regression models.
Parameters
----------
n_sectors : int
Number of sectors to fit independant regression models to.
Default is 12 sectors.
quantiles : bool
Whether to use a equal probability sectors or constant
sector widths. Default is to use constant sector widths.
Equal probability sectors are computed on the reference data
passed to ``fit``, pooled over all points, and sector 0 starts
at north rather than being centred on it.
Attributes
----------
sector_edges : numpy.ndarray or None
Interior sector cut points cached by :py:meth:`wd_to_sector` and
reused at predict time: ``n_sectors`` ceilings for constant width
sectors, ``n_sectors - 1`` interior edges for equal probability
sectors. ``None`` until the first call.
"""
def __init__(self, n_sectors=12, quantiles=False):
self.n_sectors = n_sectors
self.quantiles = quantiles
self.sector_edges = None
def wd_to_sector(self, wd, recalc_sectors=False):
"""
Convert wind direction to sector indices using the object's
n_sectors and quantiles settings.
The sector cut points are computed on ``wd`` when ``recalc_sectors``
is True or no cut points are cached yet, and cached on the object.
Later calls reuse the cached cut points, so predictions are binned
with the sectors the model was fitted with.
Parameters
----------
wd : array, float, shape (n_samples,)
Wind directions.
recalc_sectors : bool
Recompute and cache the sector cut points from ``wd``.
Default is False.
Returns
-------
array, float, shape (n_samples, )
Wind direction sector indices, NaN where ``wd`` is not finite.
"""
if recalc_sectors or self.sector_edges is None:
_, sector_coords = wd_to_sector(
wd,
sectors=self.n_sectors,
output_type="indices",
quantiles=self.quantiles,
)
if self.quantiles:
# Quantile sector 0 starts at north, so the cuts are the
# floors of sectors 1..n_sectors-1 and no wrap-around is needed.
self.sector_edges = sector_coords.sector_floor.values[1:]
else:
# Constant width sector 0 straddles north, so every ceiling is a
# cut and the last bin wraps back onto sector 0 via the modulo.
self.sector_edges = sector_coords.sector_ceil.values
# Fit and predict both bin against the cached cuts, so they agree by
# construction. Wind directions are in [0, 360], so np.digitize returns
# values in [0, len(self.sector_edges)]; the modulo only acts for constant
# width sectors, where wd >= last edge maps back to sector 0 (North).
wd = np.asarray(wd, dtype=np.float64)
sector = (np.digitize(wd, self.sector_edges) % self.n_sectors).astype(
np.float64
)
sector[~np.isfinite(wd)] = np.nan
return sector
class MCPRegressor(MCPBase):
"""
Sectorwise Regressor base-class used as a blueprint
for sectorwise MCP methods based on regression models.
Parameters
----------
ws_cutoff : float
model_kws : dict
kwargs to input when the model is instanciated.
Attributes
----------
model: what model to use in each sector
models : list, shape (n_sectors,)
one instanciated model per sector
"""
model = None
def __init__(self, ws_cutoff=0.0, model_kws=None, **kwargs):
super().__init__(**kwargs)
if model_kws is None:
model_kws = {}
self.ws_cutoff = ws_cutoff
self.model_kws = model_kws
self.models_ = [self.model(**model_kws) for _ in range(self.n_sectors)]
@ltc_validate(n_args=2)
def fit(self, ds_ref, ds_target, **kwargs):
"""
Fit the reference data to the target data by looping
through the sectors and fitting the model to the
data for that sector.
Parameters
----------
ds_ref: xarray.Dataset
Windkit time series wind climate dataset with wind speed
and wind direction for the reference data.
ds_target: xarray.Dataset
Windkit time series wind climate dataset with wind speed
and wind direction for the reference data.
Returns
-------
self : returns an instance of self.
"""
# check valid values
assert_valid_values(ds_ref)
assert_valid_values(ds_target)
assert_timesteps_equal(ds_ref, ds_target)
sector_ref = self.wd_to_sector(
ds_ref.wind_direction.values.flatten(), recalc_sectors=True
)
for i_sec in range(self.n_sectors):
x = ds_ref.wind_speed.values.flatten()
y = ds_target.wind_speed.values.flatten()
mask_sec = sector_ref == i_sec
x = x[mask_sec]
y = y[mask_sec]
mask_na = np.isnan(x) | np.isnan(y)
x = x[~mask_na]
y = y[~mask_na]
mask_ws = x < self.ws_cutoff
x = x[~mask_ws]
y = y[~mask_ws]
X = x[:, np.newaxis]
self.models_[i_sec].fit(X, y)
return self
@ltc_validate(n_args=1)
def predict(self, ds_ref, return_std=False, **kwargs):
"""Predict wind speed from the reference data.
Loops through sectors and predicts with the sector-specific model.
Parameters
----------
ds_ref : xarray.Dataset
Time series wind climate dataset.
return_std : bool
If True, return a tuple ``(ds_pred, std)`` where ``std`` is a
flat array of residual standard deviations (one per element in
the flattened wind speed array) mapped from the sector models.
Default is False.
Returns
-------
xarray.Dataset or tuple of (xarray.Dataset, numpy.ndarray)
Predicted wind speeds and directions. If *return_std* is True,
a tuple ``(ds_pred, std)`` is returned instead.
"""
assert_valid_values(ds_ref)
ws_pred = np.full_like(
ds_ref.wind_speed.values.flatten(), np.nan, dtype=np.float64
)
ds_pred = ds_ref.copy()
sector_ref = self.wd_to_sector(
ds_ref.wind_direction.values.flatten(), recalc_sectors=False
)
if return_std:
if not getattr(self.models_[0], "supports_return_std", False):
raise NotImplementedError(
f"return_std is not supported for {type(self).__name__} "
f"because {type(self.models_[0]).__name__} does not "
f"compute residual_std during fit()."
)
std_pred = np.full_like(ws_pred, np.nan)
for i_sec in range(self.n_sectors):
x = ds_ref.wind_speed.values.flatten()
mask_sec = sector_ref == i_sec
mask_na = np.isnan(x)
mask = mask_sec & ~mask_na
X = x[mask, np.newaxis]
ws_pred[mask] = self.models_[i_sec].predict(X)
if return_std:
std_pred[mask] = self.models_[i_sec].residual_std
ds_pred["wind_speed"].values = ws_pred.reshape(ds_pred.wind_speed.shape)
if return_std:
return ds_pred, std_pred
return ds_pred
[docs]
class LinRegMCP(MCPRegressor):
"""Sectorwise MCP using Scipy's linear regression.
Parameters
----------
n_sectors : int
Number of sectors to fit independant regression models to.
Default is 12 sectors.
quantiles : bool
Use equal probability sectors instead of constant widths. The
sector edges are the quantiles of the reference wind directions
passed to ``fit``, pooled over all points, and are reused at
``predict``. Sector 0 starts at north rather than being centred
on it. Default is False.
ws_cutoff : float
model_kws : dict
kwargs to input when the model is instanciated.
Attributes
----------
model: what model to use in each sector
models : list, shape (n_sectors,)
one instanciated model per sector
"""
model = LinearRegression
[docs]
def __init__(self, ws_cutoff=0.0, **kwargs):
super().__init__(ws_cutoff=ws_cutoff, **kwargs)
[docs]
def predict_with_noise(self, ds_ref, seed=0):
"""Predict wind speed with added Gaussian noise.
Calls :meth:`predict` to obtain the deterministic prediction and the
per-timestep residual standard deviation, then samples noise from
``N(0, residual_std)`` and adds it to the predicted wind speed.
Negative wind speeds are clipped to zero.
Parameters
----------
ds_ref : xarray.Dataset
Time series wind climate dataset.
seed : int
Random seed for reproducibility. Default is 0.
Returns
-------
xarray.Dataset
Predicted wind speeds with added noise.
Raises
------
WindkitValidationError
If any sector has an undefined residual standard deviation
(too few training samples).
"""
unfitted = [i for i, m in enumerate(self.models_) if m.residual_std is None]
if unfitted:
raise WindkitValidationError(
"Cannot add noise: model has not been fitted yet."
)
nan_sectors = [
i for i, m in enumerate(self.models_) if np.isnan(m.residual_std)
]
if nan_sectors:
raise WindkitValidationError(
f"Cannot add noise: sectors {nan_sectors} have undefined "
f"residual_std (too few training samples)."
)
ds_pred, std = self.predict(ds_ref, return_std=True)
rng = np.random.default_rng(seed)
noise = rng.normal(0, std)
ws_noisy = ds_pred.wind_speed.values.flatten() + noise
ws_noisy = np.clip(ws_noisy, 0, None)
ds_pred["wind_speed"].values = ws_noisy.reshape(ds_pred.wind_speed.shape)
return ds_pred
[docs]
class VarRatMCP(MCPRegressor):
"""
Sectorwise MCP using variance ratio linear regression.
Parameters
----------
n_sectors : int
Number of sectors to fit independant regression models to.
Default is 12 sectors.
quantiles : bool
Use equal probability sectors instead of constant widths. The
sector edges are the quantiles of the reference wind directions
passed to ``fit``, pooled over all points, and are reused at
``predict``. Sector 0 starts at north rather than being centred
on it. Default is False.
ws_cutoff : float
model_kws : dict
kwargs to input when the model is instanciated.
Attributes
----------
model: what model to use in each sector
models : list, shape (n_sectors,)
one instanciated model per sector
"""
model = VarianceRatioLinearRegression
[docs]
def __init__(self, fit_intercept=True, ws_cutoff=0.0, **kwargs):
model_kws = {"fit_intercept": fit_intercept}
super().__init__(model_kws=model_kws, ws_cutoff=ws_cutoff, **kwargs)
