[Fea]Add trapezoid integrate (PaddlePaddle#729)

* [Fea]Add trapezoid integrate

* fix doc

* update 1

docs/zh/api/experimental.md

@@ -14,5 +14,6 @@
         - bessel_i1e
         - fractional_diff
         - gaussian_integrate
+        - trapezoid_integrate
       show_root_heading: true
       heading_level: 3

ppsci/experimental/__init__.py

@@ -22,6 +22,7 @@
 from ppsci.experimental.math_module import bessel_i1e
 from ppsci.experimental.math_module import fractional_diff
 from ppsci.experimental.math_module import gaussian_integrate
+from ppsci.experimental.math_module import trapezoid_integrate
 
 __all__ = [
     "bessel_i0",
@@ -30,4 +31,5 @@
     "bessel_i1e",
     "fractional_diff",
     "gaussian_integrate",
+    "trapezoid_integrate",
 ]

ppsci/experimental/math_module.py

@@ -396,3 +396,69 @@ def int_func(s):
         int_func, dim=1, N=2**10 + 1, integration_domains=[[a, t]], dtype=dtype
     )
     return result
+
+
+def trapezoid_integrate(
+    y: paddle.Tensor,
+    x: paddle.Tensor = None,
+    dx: float = None,
+    axis: int = -1,
+    mode: Literal["sum", "cumsum"] = "sum",
+) -> paddle.Tensor:
+    """
+    Integrate along the given axis using the composite trapezoidal rule. Use the sum method.
+
+    Args:
+        y (paddle.Tensor): Input to be integrated.
+        x (paddle.Tensor, optional): The sample points corresponding to the input samples. its shape should be,
+        (1) input.shape; (2) the input.shape[axis] if axis is not default. Defaults to None.
+        dx (float, optional): The sample points are assumed to be evenly spaced and it is the spacing between sample points.
+        If 'x' and 'dx' are both default, 'dx' is set to 1 by default. Defaults to None.
+        axis (int, optional): The axis along which to integrate. Defaults to -1.
+        mode (Literal["sum", "cumsum"], optional): Which type cumulative sum function used. Defaults to "sum".
+
+    Returns:
+        paddle.Tensor: Integral result. If dim of input is N, return is N-1 dim.
+
+    Examples:
+        >>> import paddle
+        >>> import ppsci
+        >>> y = paddle.to_tensor([[0, 1, 2], [3, 4, 5]], dtype="float32")
+        >>> res = ppsci.experimental.trapezoid_integrate(y)
+        >>> print(res)
+        Tensor(shape=[2], dtype=float32, place=Place(gpu:0), stop_gradient=True,
+            [2., 8.])
+        >>> res = ppsci.experimental.trapezoid_integrate(y, mode="cumsum")
+        >>> print(res)
+        Tensor(shape=[2, 2], dtype=float32, place=Place(gpu:0), stop_gradient=True,
+            [[0.50000000, 2.        ],
+            [3.50000000, 8.        ]])
+        >>> res = ppsci.experimental.trapezoid_integrate(
+            y, x=paddle.to_tensor([[0, 1, 2], [3, 4, 5]], dtype="float32")
+        )
+        >>> print(res)
+        Tensor(shape=[2], dtype=float32, place=Place(gpu:0), stop_gradient=True,
+            [2., 8.])
+        >>> res = ppsci.experimental.trapezoid_integrate(
+            y, x=paddle.to_tensor([0, 1], dtype="float32"), axis=0
+        )
+        >>> print(res)
+        Tensor(shape=[3], dtype=float32, place=Place(gpu:0), stop_gradient=True,
+            [1.50000000, 2.50000000, 3.50000000])
+        >>> res = ppsci.experimental.trapezoid_integrate(
+            y, x=paddle.to_tensor([0, 1, 2], dtype="float32"), axis=1
+        )
+        >>> print(res)
+        Tensor(shape=[2], dtype=float32, place=Place(gpu:0), stop_gradient=True,
+            [2., 8.])
+        >>> res = ppsci.experimental.trapezoid_integrate(y, dx=2)
+        >>> print(res)
+        Tensor(shape=[2], dtype=float32, place=Place(gpu:0), stop_gradient=True,
+            [4. , 16.])
+    """
+    if mode == "sum":
+        return paddle.trapezoid(y, x, dx, axis)
+    elif mode == "cumsum":
+        return paddle.cumulative_trapezoid(y, x, dx, axis)
+    else:
+        raise ValueError(f'mode should be "sum" or "cumsum", but got {mode}')

test/experimental/test_trapezoid_integrate.py

@@ -0,0 +1,86 @@
+from typing import Callable
+
+import numpy as np
+import paddle
+import pytest
+from typing_extensions import Literal
+
+from ppsci.experimental import trapezoid_integrate
+
+paddle.seed(1024)
+
+
+def trapezoid_sum_test(y, x, dx):
+    dx = 1 if not dx else dx
+    y = y.numpy()
+    res = []
+    for i in range(len(y) - 1):
+        res.append((y[i] + y[i + 1]) * dx / 2)
+    return np.sum(np.array(res))
+
+
+def trapezoid_cum_test(y, x, dx):
+    dx = 1 if not dx else dx
+    y = y.numpy()
+    res = []
+    for i in range(len(y) - 1):
+        res.append((y[i] + y[i + 1]) * dx / 2)
+    return np.cumsum(np.array(res))
+
+
+def trapezoid_x_test(y, x, dx):
+    dx = 1 if not dx else dx
+    y = y.numpy()
+    res = []
+    for yi in y:
+        res_i = []
+        for i in range(len(yi) - 1):
+            res_i.append((yi[i] + yi[i + 1]) * (x[i + 1] - x[i]) / 2)
+        res.append(res_i)
+    return np.sum(np.array(res), axis=1)
+
+
+@pytest.mark.parametrize(
+    "y,x,dx,axis,mode,antideriv_func",
+    [
+        (
+            paddle.to_tensor([0, 1, 2, 3, 4, 5], dtype="float32"),
+            None,
+            None,
+            -1,
+            "sum",
+            trapezoid_sum_test,
+        ),
+        (
+            paddle.to_tensor([0, 1, 2, 3, 4, 5], dtype="float32"),
+            None,
+            2,
+            -1,
+            "cumsum",
+            trapezoid_cum_test,
+        ),
+        (
+            paddle.to_tensor([[0, 1, 2], [3, 4, 5]], dtype="float32"),
+            paddle.to_tensor([0, 1, 2], dtype="float32"),
+            None,
+            1,
+            "sum",
+            trapezoid_x_test,
+        ),
+    ],
+)
+def test_trapezoid_integrate(
+    y: paddle.Tensor,
+    x: paddle.Tensor,
+    dx: float,
+    axis: int,
+    mode: Literal["sum", "cumsum"],
+    antideriv_func: Callable,
+):
+    integrate_result = trapezoid_integrate(y, x, dx, axis, mode)
+    reference_result = antideriv_func(y, x, dx)
+    assert np.allclose(integrate_result.numpy(), reference_result)
+
+
+if __name__ == "__main__":
+    pytest.main()