"""This module contains the noise models available in Pastas.
A Noise model may be used to transform the residual series into a noise series that
better represents white noise.
Examples
--------
By default, a noise model is added to a Pastas model. It is possible to replace the
default model with different models as follows:
>>> n = ps.ArmaModel()
>>> ml.add_noisemodel(n)
or, to delete the noise model from the model:
>>> ml.del_noisemodel()
See Also
--------
pastas.model.Model.add_noisemodel
"""
# Type Hinting
from typing import Optional
import numpy as np
from pandas import DataFrame, Series, Timedelta
from pastas.typing import ArrayLike
from .decorators import njit, set_parameter
__all__ = ["NoiseModel", "ArmaModel"]
[docs]class NoiseModelBase:
_name = "NoiseModelBase"
[docs] def __init__(self) -> None:
self.nparam = 1
self.name = "noise"
self.parameters = DataFrame(columns=["initial", "pmin", "pmax", "vary", "name"])
[docs] def set_init_parameters(self, oseries: Optional[Series] = None) -> None:
if oseries is not None:
pinit = np.diff(oseries.index.to_numpy()) / Timedelta("1D")
pinit = np.median(pinit)
else:
pinit = 14.0
self.parameters.loc["noise_alpha"] = (pinit, 1e-5, 5000.0, True, "noise")
@set_parameter
def _set_initial(self, name: str, value: float) -> None:
"""Internal method to set the initial parameter value.
Notes
-----
The preferred method for parameter setting is through the model.
"""
self.parameters.loc[name, "initial"] = value
@set_parameter
def _set_pmin(self, name: str, value: float) -> None:
"""Internal method to set the minimum value of the noisemodel.
Notes
-----
The preferred method for parameter setting is through the model.
"""
self.parameters.loc[name, "pmin"] = value
@set_parameter
def _set_pmax(self, name: str, value: float) -> None:
"""Internal method to set the maximum parameter values.
Notes
-----
The preferred method for parameter setting is through the model.
"""
self.parameters.loc[name, "pmax"] = value
@set_parameter
def _set_vary(self, name: str, value: float) -> None:
"""Internal method to set if the parameter is varied.
Notes
-----
The preferred method for parameter setting is through the model.
"""
self.parameters.loc[name, "vary"] = value
[docs] def to_dict(self) -> dict:
"""Method to return a dict to store the noise model"""
data = {"class": self._name, "norm": self.norm}
return data
[docs] @staticmethod
def weights(res, p) -> int:
return 1
[docs]class NoiseModel(NoiseModelBase):
"""Noise model with exponential decay of the residuals and weighting.
Parameters
----------
norm: boolean, optional
Boolean to indicate whether weights are normalized according to the Von
Asmuth and Bierkens (2005) paper. Default is True.
Notes
-----
Calculates the noise :cite:t:`von_asmuth_modeling_2005` according to:
.. math::
v(t_1) = r(t_1) - r(t_0) * \\exp(- \\Delta t / \\alpha)
Calculates the weights as
.. math::
w = 1 / \\sqrt{(1 - \\exp(-2 \\Delta t / \\alpha))}
The units of the alpha parameter is always in days. The first value of the noise
is the residual ($v(t=0=r(t=0)$). First weight is 1 / sig_residuals (i.e.,
delt = infty). Normalization of weights as in :cite:t:`von_asmuth_modeling_2005`,
optional.
"""
_name = "NoiseModel"
[docs] def __init__(self, norm: bool = True) -> None:
NoiseModelBase.__init__(self)
self.norm = norm
self.nparam = 1
self.set_init_parameters()
[docs] @staticmethod
def simulate(res: Series, p: ArrayLike) -> Series:
"""Simulate noise from the residuals.
Parameters
----------
res: pandas.Series
The residual series.
p: array_like
array_like object with the values as floats representing the model
parameters. Here, Alpha parameter used by the noisemodel.
Returns
-------
noise: pandas.Series
Series of the noise.
"""
alpha = p[0]
odelt = np.diff(res.index.to_numpy()) / Timedelta("1D")
v = np.append(
res.values[0], res.values[1:] - np.exp(-odelt / alpha) * res.values[:-1]
)
return Series(data=v, index=res.index, name="Noise")
[docs] def weights(self, res: Series, p: ArrayLike) -> Series:
"""Method to calculate the weights for the noise.
Parameters
----------
res: pandas.Series
Pandas Series with the residuals to compute the weights for. The Series
index must be a DatetimeIndex.
p: array_like
NumPy array with the parameters used in the noise model.
Returns
-------
w: pandas.Series
Series of the weights.
Notes
-----
Weights are
.. math:: w = 1 / sqrt((1 - exp(-2 \\Delta t / \\alpha)))
which are then normalized so that sum(w) = len(res).
"""
alpha = p[0]
# large for first measurement
odelt = np.append(1e12, np.diff(res.index.to_numpy()) / Timedelta("1D"))
exp = np.exp(-2.0 / alpha * odelt) # Twice as fast as 2*odelt/alpha
w = 1 / np.sqrt(1.0 - exp) # weights of noise, not noise^2
if self.norm:
w *= np.exp(1.0 / (2.0 * odelt.size) * np.sum(np.log(1.0 - exp)))
return Series(data=w, index=res.index, name="noise_weights")
[docs] def to_dict(self) -> dict:
"""Method to return a dict to store the noise model"""
data = {"class": self._name, "norm": self.norm}
return data
[docs]class ArmaModel(NoiseModelBase):
"""ARMA(1,1) Noise model to simulate the noise as defined in
:cite:t:`collenteur_estimation_2021`.
Notes
-----
Calculates the noise according to:
.. math::
\\upsilon_t = r_t - r_{t-1} e^{-\\Delta t/\\alpha} - \\upsilon_{t-1}
e^{-\\Delta t/\\beta}
The units of the alpha and beta parameters are always in days.
Warnings
--------
This model has only been tested on regular time steps and should not be used for
irregular time steps yet.
"""
_name = "ArmaModel"
[docs] def __init__(self) -> None:
NoiseModelBase.__init__(self)
self.nparam = 2
self.set_init_parameters()
[docs] def set_init_parameters(self, oseries: Series = None) -> None:
if oseries is not None:
pinit = np.diff(oseries.index.to_numpy()) / Timedelta("1D")
pinit = np.median(pinit)
else:
pinit = 14.0
self.parameters.loc["noise_alpha"] = (pinit, 1e-9, 5000.0, True, "noise")
self.parameters.loc["noise_beta"] = (1.0, -np.inf, np.inf, True, "noise")
[docs] def simulate(self, res: Series, p: ArrayLike) -> Series:
alpha = p[0]
beta = p[1]
# Calculate the time steps
odelt = np.diff(res.index.to_numpy()) / Timedelta("1D")
a = self.calculate_noise(res.values, odelt, alpha, beta)
return Series(index=res.index, data=a, name="Noise")
[docs] @staticmethod
@njit
def calculate_noise(
res: ArrayLike, odelt: ArrayLike, alpha: float, beta: float
) -> ArrayLike:
# Create an array to store the noise
a = np.zeros_like(res)
a[0] = res[0]
if beta == 0.0: # Prevent division by zero errors
beta = 1e-24
pm = beta / np.abs(beta)
# We have to loop through each value
for i in range(1, res.size):
a[i] = (
res[i]
- res[i - 1] * np.exp(-odelt[i - 1] / alpha)
- a[i - 1] * pm * np.exp(-odelt[i - 1] / np.abs(beta))
)
return a