TB: added option for periodic grids to sparse interpolation matrices (#19)

* Added option for periodic grids to sparse interpolation matrices

* Using camel case

Author: brownbaerchen

qmat/lagrange.py

@@ -537,7 +537,7 @@ def getDerivationMatrix(self, *args, **kwargs):
         return self.getDerivativeMatrix(*args, **kwargs)
 
 
-def getSparseInterpolationMatrix(inPoints, outPoints, order):
+def getSparseInterpolationMatrix(inPoints, outPoints, order, gridPeriod=-1):
     """
     Get a sparse interpolation matrix from `inPoints` to `outPoints` of order
     `order` using barycentric Lagrange interpolation.
@@ -553,6 +553,8 @@ def getSparseInterpolationMatrix(inPoints, outPoints, order):
             The points you want to interpolate to
         order : int
             Order of the interpolation
+        grid_period : float
+            Period of the grid. Negative values indicate non-periodic grids
 
     Returns
     -------
@@ -568,10 +570,21 @@ def getSparseInterpolationMatrix(inPoints, outPoints, order):
     lastClosestPoints = None
 
     for i in range(len(outPoints)):
-        closestPointsIdx = np.sort(np.argsort(np.abs(inPoints - outPoints[i]))[:order])
-        closestPoints = inPoints[closestPointsIdx] - outPoints[i]
+        if gridPeriod > 0:
+            pathL = (inPoints - gridPeriod - outPoints[i] % gridPeriod)
+            pathR = (inPoints + gridPeriod - outPoints[i] % gridPeriod)
+            pathC = (inPoints - outPoints[i] % gridPeriod)
+            path = np.append(np.append(pathR, pathL), pathC)
+            dist = np.abs(path)
+            _closestPointsIdx = np.sort(np.argsort(dist)[:order])
+            closestPointsIdx = _closestPointsIdx % len(inPoints)
+            closestPoints, sorting = np.unique(path[_closestPointsIdx], return_index=True)
+            closestPointsIdx = closestPointsIdx[sorting]
+        else:
+            closestPointsIdx = np.sort(np.argsort(np.abs(inPoints - outPoints[i]))[:order])
+            closestPoints = inPoints[closestPointsIdx] - outPoints[i]
 
-        if lastClosestPoints is not None and np.allclose(closestPoints, lastClosestPoints):
+        if lastClosestPoints is not None and len(closestPoints) == len(lastClosestPoints) and np.allclose(closestPoints, lastClosestPoints):
             interpolationLine = lastInterpolationLine
         else:
             interpolator = LagrangeApproximation(points = closestPoints)

tests/test_2_lagrange.py

@@ -203,3 +203,31 @@ def testSparseInterpolation(inPoints, outPoints, order):
         inPolynomial = np.polyval(polyCoeffs,inPoints)
         interpolated = interpolationMatrix @ inPolynomial
         assert not np.allclose(np.polyval(polyCoeffs, outPoints), interpolated), f'Interpolation of order {order+1} polynomial is unexpectedly exact'
+
+
+@pytest.mark.parametrize('inPoints', [np.linspace(0, 1, 64, endpoint=False)])
+@pytest.mark.parametrize('outPoints', [np.linspace(0.123, 4.1214, 256, endpoint=False), np.linspace(0, 1, 128, endpoint=False)])
+@pytest.mark.parametrize('order', [2, 3, 4])
+def testSparseInterpolationPeriodic(inPoints, outPoints, order, gridPeriod=1):
+    """
+    In this test, we do "extrapolation", which is really interpolation because of the periodicity of the grid.
+    """
+    from qmat.lagrange import getSparseInterpolationMatrix
+    import scipy.sparse as sp
+
+    data = np.sin(inPoints * 2 * np.pi / gridPeriod)
+
+    interpolationMatrix = getSparseInterpolationMatrix(inPoints, outPoints, order, gridPeriod=gridPeriod)
+    assert isinstance(interpolationMatrix, sp.csc_matrix)
+
+
+    interpolated = interpolationMatrix @ data
+    ref = np.sin(outPoints * 2* np.pi / gridPeriod)
+    error = np.linalg.norm(interpolated - ref)
+
+    max_error = {
+        4: 2e-5,
+        3: 4e-4,
+        2: 8e-3,
+    }
+    assert error < max_error[order], error