Support Lp norms for surfaces.

Remove more Simpson quadrature calls.

skfda/exploratory/stats/_stats.py

@@ -8,6 +8,8 @@
 from scipy import integrate
 from scipy.stats import rankdata
 
+from skfda._utils.ndfunction import average_function_value
+
 from ...misc.metrics._lp_distances import l2_distance
 from ...representation import FData, FDataGrid
 from ...typing._metric import Metric
@@ -108,33 +110,20 @@ def modified_epigraph_index(X: FDataGrid) -> NDArrayFloat:
     with all the other curves of our dataset.
 
     """
-    interval_len = (
-        X.domain_range[0][1]
-        - X.domain_range[0][0]
-    )
-
-    # Array containing at each point the number of curves
+    # Functions containing at each point the number of curves
     # are above it.
-    num_functions_above: NDArrayFloat = rankdata(
-        -X.data_matrix,
-        method='max',
-        axis=0,
-    ) - 1
-
-    integrand = num_functions_above
-
-    for d, s in zip(X.domain_range, X.grid_points):
-        integrand = integrate.simpson(
-            integrand,
-            x=s,
-            axis=1,
-        )
-        interval_len = d[1] - d[0]
-        integrand /= interval_len
-
-    integrand /= X.n_samples
+    num_functions_above = X.copy(
+        data_matrix=rankdata(
+            -X.data_matrix,
+            method='max',
+            axis=0,
+        ) - 1,
+    )
 
-    return integrand.flatten()
+    return (
+        average_function_value(num_functions_above)
+        / num_functions_above.n_samples
+    ).ravel()
 
 
 def depth_based_median(

skfda/misc/metrics/_lp_norms.py

@@ -7,7 +7,7 @@
 import scipy.integrate
 from typing_extensions import Final
 
-from ...representation import FData, FDataBasis
+from ...representation import FData, FDataBasis, FDataGrid
 from ...typing._metric import Norm
 from ...typing._numpy import NDArrayFloat
 
@@ -135,20 +135,21 @@ def __call__(self, vector: Union[NDArrayFloat, FData]) -> NDArrayFloat:
             )
             res = np.sqrt(integral[0]).flatten()
 
-        else:
+        elif isinstance(vector, FDataGrid):
             data_matrix = vector.data_matrix
-            original_shape = data_matrix.shape
-            data_matrix = data_matrix.reshape(-1, original_shape[-1])
 
-            data_matrix = (np.linalg.norm(
-                vector.data_matrix,
-                ord=vector_norm,
-                axis=-1,
-                keepdims=True,
-            ) if isinstance(vector_norm, (float, int))
-                else vector_norm(data_matrix)
-            )
-            data_matrix = data_matrix.reshape(original_shape[:-1] + (1,))
+            if isinstance(vector_norm, (float, int)):
+                data_matrix = np.linalg.norm(
+                    vector.data_matrix,
+                    ord=vector_norm,
+                    axis=-1,
+                    keepdims=True,
+                )
+            else:
+                original_shape = data_matrix.shape
+                data_matrix = data_matrix.reshape(-1, original_shape[-1])
+                data_matrix = vector_norm(data_matrix)
+                data_matrix = data_matrix.reshape(original_shape[:-1] + (1,))
 
             if np.isinf(self.p):
 
@@ -157,18 +158,19 @@ def __call__(self, vector: Union[NDArrayFloat, FData]) -> NDArrayFloat:
                     axis=tuple(range(1, data_matrix.ndim)),
                 )
 
-            elif vector.dim_domain == 1:
+            else:
 
+                integrand = vector.copy(
+                    data_matrix=data_matrix ** self.p,
+                    coordinate_names=(None,),
+                )
                 # Computes the norm, approximating the integral with Simpson's
                 # rule.
-                res = scipy.integrate.simpson(
-                    data_matrix[..., 0] ** self.p,
-                    x=vector.grid_points,
-                ) ** (1 / self.p)
-
-            else:
-                # Needed to perform surface integration
-                return NotImplemented
+                res = integrand.integrate().ravel() ** (1 / self.p)
+        else:
+            raise NotImplementedError(
+                f"LpNorm not implemented for type {type(vector)}",
+            )
 
         if len(res) == 1:
             return res[0]  # type: ignore[no-any-return]

skfda/tests/test_metrics.py

@@ -89,10 +89,11 @@ def test_lp_norm_surface_inf(self) -> None:
         )
 
     def test_lp_norm_surface(self) -> None:
-        """Test that integration of surfaces has not been implemented."""
-        # Integration of surfaces not implemented, add test case after
-        # implementation
-        self.assertEqual(lp_norm(self.fd_surface, p=1), NotImplemented)
+        """Test the integration of surfaces."""
+        np.testing.assert_allclose(
+            lp_norm(self.fd_surface, p=1),
+            [0.12566256, 0.12563755, 0.12566133],
+        )
 
     def test_lp_error_dimensions(self) -> None:
         """Test error on metric between different kind of objects."""