Change simps for simpson.

Extract also average_function_value.

skfda/_utils/ndfunction/_functions.py

@@ -0,0 +1,116 @@
+"""Functions for working with arrays of functions."""
+from __future__ import annotations
+
+import math
+from typing import (
+    TYPE_CHECKING,
+    Any,
+    Literal,
+    Protocol,
+    TypeVar,
+    runtime_checkable,
+)
+
+from ...misc.validation import validate_domain_range
+from .. import nquad_vec
+
+if TYPE_CHECKING:
+    from ...typing._numpy import NDArrayFloat
+    from ...representation import FData
+    from ...typing._base import DomainRangeLike
+
+
+UfuncMethod = Literal[
+    "__call__",
+    "reduce",
+    "reduceat",
+    "accumulate",
+    "outer",
+    "inner",
+]
+
+
+class _SupportsArrayUFunc(Protocol):
+    def __array_ufunc__(
+        self,
+        ufunc: Any,
+        method: UfuncMethod,
+        *inputs: Any,
+        **kwargs: Any,
+    ) -> Any:
+        pass
+
+
+T = TypeVar("T", bound=_SupportsArrayUFunc)
+
+
+class _UnaryUfunc(Protocol):
+
+    def __call__(self, __arg: T) -> T:  # noqa: WPS112
+        pass
+
+
+def _average_function_ufunc(
+    data: FData,
+    ufunc: _UnaryUfunc,
+    *,
+    domain: DomainRangeLike | None = None,
+) -> NDArrayFloat:
+
+    if domain is None:
+        domain = data.domain_range
+    else:
+        domain = validate_domain_range(domain)
+
+    lebesgue_measure = math.prod(
+        (
+            (iterval[1] - iterval[0])
+            for iterval in domain
+        ),
+    )
+
+    try:
+        data_eval = ufunc(data)
+    except TypeError:
+
+        def integrand(*args: NDArrayFloat) -> NDArrayFloat:  # noqa: WPS430
+            f1 = data(args)[:, 0, :]
+            return ufunc(f1)
+
+        return nquad_vec(
+            integrand,
+            domain,
+        ) / lebesgue_measure
+
+    else:
+        return data_eval.integrate(domain=domain) / lebesgue_measure
+
+
+def average_function_value(
+    data: FData,
+    *,
+    domain: DomainRangeLike | None = None,
+) -> NDArrayFloat:
+    r"""
+    Calculate the average function value for each function.
+
+    This is the value that, if integrated over the whole domain of each
+    function, has the same integral as the function itself.
+
+    .. math::
+        \bar{x} = \frac{1}{\text{Vol}(\mathcal{T})}\int_{\mathcal{T}} x(t) dt
+
+    Args:
+        data: Functions where we want to calculate the expected value.
+        domain: Integration domain. By default, the whole domain is used.
+
+    Returns:
+        ndarray of shape (n_dimensions, n_samples) with the values of the
+        expectations.
+
+    See also:
+        `Entry on Wikipedia
+        <https://en.wikipedia.org/wiki/Mean_of_a_function>`_
+
+    """
+    return _average_function_ufunc(data, ufunc=lambda x: x, domain=domain)

skfda/exploratory/depth/_depth.py

@@ -92,7 +92,7 @@ def transform(self, X: FDataGrid) -> NDArrayFloat:  # noqa: D102
         integrand = pointwise_depth
 
         for d, s in zip(X.domain_range, X.grid_points):
-            integrand = scipy.integrate.simps(
+            integrand = scipy.integrate.simpson(
                 integrand,
                 x=s,
                 axis=1,

skfda/exploratory/outliers/_directional_outlyingness.py

@@ -213,7 +213,7 @@ def directional_outlyingness_stats(  # noqa: WPS218
     )
     assert weighted_dir_outlyingness.shape == dir_outlyingness.shape
 
-    mean_dir_outlyingness = scipy.integrate.simps(
+    mean_dir_outlyingness = scipy.integrate.simpson(
         weighted_dir_outlyingness,
         fdatagrid.grid_points[0],
         axis=1,
@@ -230,7 +230,7 @@ def directional_outlyingness_stats(  # noqa: WPS218
         axis=-1,
     ))
     weighted_norm = norm * pointwise_weights
-    variation_dir_outlyingness = scipy.integrate.simps(
+    variation_dir_outlyingness = scipy.integrate.simpson(
         weighted_norm,
         fdatagrid.grid_points[0],
         axis=1,

skfda/exploratory/stats/_fisher_rao.py

@@ -133,7 +133,7 @@ def _fisher_rao_warping_mean(
         # Compute shooting vectors
         for psi_i in psi_data:
 
-            inner = scipy.integrate.simps(mu * psi_i, x=eval_points)
+            inner = scipy.integrate.simpson(mu * psi_i, x=eval_points)
             inner = max(min(inner, 1), -1)
 
             theta = np.arccos(inner)
@@ -143,7 +143,7 @@ def _fisher_rao_warping_mean(
 
         # Mean of shooting vectors
         vmean /= warping.n_samples
-        v_norm = np.sqrt(scipy.integrate.simps(np.square(vmean)))
+        v_norm = np.sqrt(scipy.integrate.simpson(np.square(vmean)))
 
         # Convergence criterion
         if v_norm < tol:
@@ -266,7 +266,7 @@ def fisher_rao_karcher_mean(
     # Initialize with function closest to the L2 mean with the L2 distance
     centered = (srsf.T - srsf.mean(axis=0, keepdims=True).T).T
 
-    distances = scipy.integrate.simps(
+    distances = scipy.integrate.simpson(
         np.square(centered, out=centered),
         eval_points_normalized,
         axis=1,
@@ -304,14 +304,14 @@ def fisher_rao_karcher_mean(
 
         # Convergence criterion
         mu_norm = np.sqrt(
-            scipy.integrate.simps(
+            scipy.integrate.simpson(
                 np.square(mu, out=mu_aux),
                 eval_points_normalized,
             ),
         )
 
         mu_diff = np.sqrt(
-            scipy.integrate.simps(
+            scipy.integrate.simpson(
                 np.square(mu - mu_1, out=mu_aux),
                 eval_points_normalized,
             ),

skfda/exploratory/stats/_stats.py

@@ -124,7 +124,7 @@ def modified_epigraph_index(X: FDataGrid) -> NDArrayFloat:
     integrand = num_functions_above
 
     for d, s in zip(X.domain_range, X.grid_points):
-        integrand = integrate.simps(
+        integrand = integrate.simpson(
             integrand,
             x=s,
             axis=1,

skfda/misc/_math.py

@@ -363,7 +363,7 @@ def _inner_product_fdatagrid(
         # Perform quadrature inside the einsum
         for i, s in enumerate(arg1.grid_points[::-1]):
             identity = np.eye(len(s))
-            weights = scipy.integrate.simps(identity, x=s)
+            weights = scipy.integrate.simpson(identity, x=s)
             index = (slice(None),) + (np.newaxis,) * (i + 1)
             d1 *= weights[index]
 

skfda/misc/metrics/_fisher_rao.py

@@ -234,7 +234,7 @@ def fisher_rao_amplitude_distance(
         penalty = np.sqrt(penalty, out=penalty)
         penalty -= 1
         penalty = np.square(penalty, out=penalty)
-        penalty = scipy.integrate.simps(penalty, x=eval_points_normalized)
+        penalty = scipy.integrate.simpson(penalty, x=eval_points_normalized)
 
         distance = np.sqrt(distance**2 + lam * penalty)
 
@@ -322,7 +322,7 @@ def fisher_rao_phase_distance(
 
     derivative_warping = np.sqrt(derivative_warping, out=derivative_warping)
 
-    d = scipy.integrate.simps(derivative_warping, x=eval_points_normalized)
+    d = scipy.integrate.simpson(derivative_warping, x=eval_points_normalized)
     d = np.clip(d, -1, 1)
 
     return np.arccos(d)  # type: ignore[no-any-return]
@@ -394,7 +394,7 @@ def _fisher_rao_warping_distance(
 
     product = np.multiply(srsf_warping1, srsf_warping2, out=srsf_warping1)
 
-    d = scipy.integrate.simps(product, x=warping1.grid_points[0])
+    d = scipy.integrate.simpson(product, x=warping1.grid_points[0])
     d = np.clip(d, -1, 1)
 
     return np.arccos(d)  # type: ignore[no-any-return]

skfda/misc/metrics/_lp_norms.py

@@ -161,7 +161,7 @@ def __call__(self, vector: Union[NDArrayFloat, FData]) -> NDArrayFloat:
 
                 # Computes the norm, approximating the integral with Simpson's
                 # rule.
-                res = scipy.integrate.simps(
+                res = scipy.integrate.simpson(
                     data_matrix[..., 0] ** self.p,
                     x=vector.grid_points,
                 ) ** (1 / self.p)

skfda/ml/regression/_historical_linear_model.py

@@ -42,7 +42,7 @@ def _pairwise_fem_inner_product(
     eval_fd = fd(grid)
 
     prod = eval_fem * eval_fd
-    integral = scipy.integrate.simps(prod, grid, axis=1)
+    integral = scipy.integrate.simpson(prod, grid, axis=1)
     return np.sum(integral, axis=-1)  # type: ignore[no-any-return]
 
 

skfda/preprocessing/dim_reduction/_fpca.py

@@ -348,7 +348,7 @@ def _fit_grid(
         if self._weights is None:
             # grid_points is a list with one array in the 1D case
             identity = np.eye(len(X.grid_points[0]))
-            self._weights = scipy.integrate.simps(identity, X.grid_points[0])
+            self._weights = scipy.integrate.simpson(identity, X.grid_points[0])
         elif callable(self._weights):
             self._weights = self._weights(X.grid_points[0])
             # if its a FDataGrid then we need to reduce the dimension to 1-D

skfda/preprocessing/feature_construction/_functions.py

@@ -3,34 +3,23 @@
 from __future__ import annotations
 
 import itertools
-from typing import Optional, Sequence, Tuple, TypeVar, Union, cast, overload
+from typing import Optional, Sequence, Tuple, TypeVar, Union, cast
 
 import numpy as np
-from typing_extensions import Literal, Protocol, TypeGuard
+from typing_extensions import Literal, TypeGuard
 
-from ..._utils import nquad_vec
-from ...misc.validation import check_fdata_dimensions, validate_domain_range
+from ..._utils.ndfunction._functions import (
+    _average_function_ufunc,
+    _UnaryUfunc,
+)
+from ...misc.validation import check_fdata_dimensions
 from ...representation import FData, FDataBasis, FDataGrid
 from ...typing._base import DomainRangeLike
 from ...typing._numpy import ArrayLike, NDArrayBool, NDArrayFloat, NDArrayInt
 
 T = TypeVar("T", bound=Union[NDArrayFloat, FDataGrid])
 
 
-class _BasicUfuncProtocol(Protocol):
-
-    @overload
-    def __call__(self, __arg: FDataGrid) -> FDataGrid:  # noqa: WPS112
-        pass
-
-    @overload
-    def __call__(self, __arg: NDArrayFloat) -> NDArrayFloat:  # noqa: WPS112
-        pass
-
-    def __call__(self, __arg: T) -> T:  # noqa: WPS112
-        pass
-
-
 def _sequence_of_ints(data: Sequence[object]) -> TypeGuard[Sequence[int]]:
     """Check that every element is an integer."""
     return all(isinstance(d, int) for d in data)
@@ -332,7 +321,7 @@ def number_crossings(
 
         >>> from skfda.preprocessing.feature_construction import (
         ...     number_crossings,
-        ... ) 
+        ... )
         >>> from scipy.special import jv
         >>> import numpy as np
         >>> x_grid = np.linspace(0, 14, 14)
@@ -521,7 +510,7 @@ def unconditional_moment(
 
 def unconditional_expected_value(
     data: FData,
-    function: _BasicUfuncProtocol,
+    function: _UnaryUfunc,
     *,
     domain: DomainRangeLike | None = None,
 ) -> NDArrayFloat:
@@ -537,19 +526,19 @@ def unconditional_expected_value(
             f_p(x(t))=\frac{1}{\left( b-a\right)}\int_a^b g
             \left(x_p(t)\right) dt
 
-        Args:
-            data: FDataGrid or FDataBasis where we want to calculate
-                the expected value.
-            function: function that specifies how the expected value to is
-                calculated. It has to be a function of X(t).
-            domain: Integration domain. By default, the whole domain is used.
+    Args:
+        data: FDataGrid or FDataBasis where we want to calculate
+            the expected value.
+        function: function that specifies how the expected value to is
+            calculated. It has to be a function of X(t).
+        domain: Integration domain. By default, the whole domain is used.
 
-        Returns:
-            ndarray of shape (n_dimensions, n_samples) with the values of the
-            expectations.
+    Returns:
+        ndarray of shape (n_dimensions, n_samples) with the values of the
+        expectations.
 
     Examples:
-        We will use this funtion to calculate the logarithmic first moment
+        We will use this function to calculate the logarithmic first moment
         of the first 5 samples of the Berkeley Growth dataset.
         We will start by importing it.
 
@@ -574,26 +563,4 @@ def unconditional_expected_value(
                    [ 4.84]])
 
     """
-    if domain is None:
-        domain = data.domain_range
-    else:
-        domain = validate_domain_range(domain)
-
-    lebesgue_measure = np.prod(
-        [
-            (iterval[1] - iterval[0])
-            for iterval in domain
-        ],
-    )
-
-    if isinstance(data, FDataGrid):
-        return function(data).integrate(domain=domain) / lebesgue_measure
-
-    def integrand(*args: NDArrayFloat) -> NDArrayFloat:  # noqa: WPS430
-        f1 = data(args)[:, 0, :]
-        return function(f1)
-
-    return nquad_vec(
-        integrand,
-        domain,
-    ) / lebesgue_measure
+    return _average_function_ufunc(data, ufunc=function, domain=domain)

skfda/representation/grid.py

@@ -515,7 +515,7 @@ def integrate(
         integrand = data.data_matrix
 
         for g in data.grid_points[::-1]:
-            integrand = scipy.integrate.simps(
+            integrand = scipy.integrate.simpson(
                 integrand,
                 x=g,
                 axis=-2,