From dbb027b9fa442e9378d4213cfaf8719fa49166ca Mon Sep 17 00:00:00 2001
From: Adrien Morison <adrien.morison@gmail.com>
Date: Sun, 26 Feb 2023 02:26:07 +0000
Subject: [PATCH] New ChebyshevSampling class

---
 dmsuite/interp.py | 56 +++++++++++++++++++++++++++++++++--------------
 1 file changed, 39 insertions(+), 17 deletions(-)

diff --git a/dmsuite/interp.py b/dmsuite/interp.py
index d5b1a80..cd1bc68 100644
--- a/dmsuite/interp.py
+++ b/dmsuite/interp.py
@@ -1,11 +1,49 @@
 """Interpolation."""
 
+from __future__ import annotations
+
+from dataclasses import dataclass
+from functools import cached_property
+
 import numpy as np
 from numpy.typing import NDArray
 
 from .poly_diff import Chebyshev
 
 
+@dataclass(frozen=True)
+class ChebyshevSampling:
+    """Barycentric polynomial interpolation on Chebyshev nodes.
+
+    Attributes:
+        degree: degree of Chebyshev polynomials.
+        positions: positions of [-1, 1] to sample.
+    """
+
+    degree: int
+    positions: NDArray
+
+    @cached_property
+    def _wgt_dif(self) -> tuple[NDArray, NDArray]:
+        # Chebyshev points
+        xxk = Chebyshev(self.degree).nodes
+        # weights for Chebyshev formula
+        wgt = (-1.0) ** np.arange(xxk.size)
+        wgt[0] /= 2
+        wgt[-1] /= 2
+        # Compute quantities xxx-xxk
+        nnx = self.positions.size
+        dif = np.tile(self.positions, (xxk.size, 1)).T - np.tile(xxk, (nnx, 1))
+        dif = 1 / (dif + np.where(dif == 0, np.finfo(float).eps, 0))
+        return wgt, dif
+
+    def apply_on(self, values: NDArray) -> NDArray:
+        """Apply desired sampling on values known at Chebyshev nodes."""
+        assert values.size == self.degree + 1
+        wgt, dif = self._wgt_dif
+        return np.dot(dif, wgt * values) / np.dot(dif, wgt)
+
+
 def chebint(ffk: NDArray, xxx: NDArray) -> NDArray:
     """Barycentric polynomial interpolation on Chebyshev nodes.
 
@@ -29,20 +67,4 @@ def chebint(ffk: NDArray, xxx: NDArray) -> NDArray:
     J.A.C. Weideman, S.C. Reddy 1998
     """
     ncheb = ffk.shape[0] - 1
-    if ncheb <= 1:
-        raise Exception("At least two data points are necessary in chebint")
-
-    nnx = xxx.shape[0]
-
-    # Chebyshev points
-    xxk = Chebyshev(degree=ncheb).nodes
-    # weights for Chebyshev formula
-    wgt = (-1.0) ** np.arange(ncheb + 1)
-    wgt[0] /= 2
-    wgt[ncheb] /= 2
-
-    # Compute quantities xxx-xxk
-    dif = np.tile(xxx, (ncheb + 1, 1)).T - np.tile(xxk, (nnx, 1))
-    dif = 1 / (dif + np.where(dif == 0, np.finfo(float).eps, 0))
-
-    return np.dot(dif, wgt * ffk) / np.dot(dif, wgt)
+    return ChebyshevSampling(degree=ncheb, positions=xxx).apply_on(ffk)