TL: refined design for QDelta generators

docs/devdoc/roadmap.md

@@ -20,6 +20,7 @@ Detailed description of all specific versions and their associated changes is av
 - ✅ integration of all RK schemes from [pySDC](https://github.com/Parallel-in-Time/pySDC)
 - ✅ citation file and zenodo releases
 - integration of `qmat` into [pySDC](https://github.com/Parallel-in-Time/pySDC)
+- ✅ refined design for $Q_\Delta$ generators
 - full documentation of classes and functions
 - finalization of extended usage tutorials ($S$-matrix, `dTau` coefficient for initial sweep, prolongation)
 - ✅ full definition and documentation of the version update pipeline

qmat/qdelta/__init__.py

@@ -7,37 +7,58 @@
 
 from qmat.utils import checkOverriding, storeClass, importAll, checkGenericConstr
 
+
 class QDeltaGenerator(object):
+    _K_DEP = False
 
     def __init__(self, Q, **kwargs):
         self.Q = np.asarray(Q, dtype=float)
-        self.QDelta = np.zeros_like(self.Q)
 
-    def storeAndReturn(self, QDelta):
-        np.copyto(self.QDelta, QDelta)
-        return self.QDelta
+    @property
+    def size(self):
+        return self.Q.shape[0]
 
-    def getQDelta(self, k=None):
+    @property
+    def zeros(self):
+        M = self.size
+        return np.zeros((M, M), dtype=float)
+
+    def computeQDelta(self, k=None) -> np.ndarray:
+        """Compute and returns the QDelta matrix"""
         raise NotImplementedError("mouahahah")
 
+    def getQDelta(self, k=None, copy=True):
+        try:
+            QDelta = self._QDelta[k] if self._K_DEP else self._QDelta
+        except Exception as e:
+            QDelta = self.computeQDelta(k)
+            if type(e) == AttributeError:
+                self._QDelta = {k: QDelta} if self._K_DEP else QDelta
+            elif type(e) == KeyError:
+                self._QDelta[k] = QDelta
+            else:
+                raise Exception("some very weird bug happened ... did you do fishy stuff ?")
+        return QDelta.copy() if copy else QDelta
+
     @property
     def dTau(self):
-        return self.QDelta[0]*0
+        return np.zeros(self.size, dtype=float)
 
     def genCoeffs(self, k=None, dTau=False):
         if isinstance(k, list):
-            out = [np.array([self.getQDelta(_k) for _k in k])]
+            out = [np.array([self.getQDelta(_k, copy=False) for _k in k])]
         else:
             out = [self.getQDelta(k)]
         if dTau:
             out += [self.dTau]
         return out if len(out) > 1 else out[0]
 
+
 QDELTA_GENERATORS = {}
 
 def register(cls:QDeltaGenerator)->QDeltaGenerator:
     checkGenericConstr(cls)
-    checkOverriding(cls, "getQDelta", isProperty=False)
+    checkOverriding(cls, "computeQDelta", isProperty=False)
     storeClass(cls, QDELTA_GENERATORS)
     return cls
 

qmat/qdelta/algebraic.py

@@ -14,50 +14,50 @@ class PIC(QDeltaGenerator):
     """Picard approximation (zeros)"""
     aliases = ["Picard"]
 
-    def getQDelta(self, k=None):
-        return self.QDelta
+    def computeQDelta(self, k=None):
+        return self.zeros
 
 
 @register
 class Exact(QDeltaGenerator):
     """Takes Q (exact approximation)"""
     aliases = ["EXACT"]
 
-    def getQDelta(self, k=None):
-        return self.storeAndReturn(self.Q)
+    def computeQDelta(self, k=None):
+        return self.Q.copy()  # TODO: do we really want a copy here ... ?
 
 
 @register
 class LU(QDeltaGenerator):
     """LU approximation from Weiser"""
 
-    def getQDelta(self, k=None):
+    def computeQDelta(self, k=None):
         _, _, U = spl.lu(self.Q.T)
-        return self.storeAndReturn(U.T)
+        return U.T
 
 
 @register
 class LU2(QDeltaGenerator):
     """LU approximation from Weiser multiplied by 2"""
 
-    def getQDelta(self, k=None):
+    def computeQDelta(self, k=None):
         _, _, U = spl.lu(self.Q.T)
-        return self.storeAndReturn(2*U.T)
+        return 2*U.T
 
 
 @register
 class QPar(QDeltaGenerator):
     """Approximation using diagonal of Q"""
     aliases = ["Qpar", "Qdiag"]
 
-    def getQDelta(self, k=None):
-        return self.storeAndReturn(np.diag(np.diag(self.Q)))
+    def computeQDelta(self, k=None):
+        return np.diag(np.diag(self.Q))
 
 
 @register
 class GS(QDeltaGenerator):
     """Approximation using lower triangular part of Q"""
     aliases = ["GaussSeidel"]
 
-    def getQDelta(self, k=None):
-        return self.storeAndReturn(np.tril(self.Q))
+    def computeQDelta(self, k=None):
+        return np.tril(self.Q)

qmat/qdelta/min.py

@@ -17,15 +17,14 @@ class MIN(QDeltaGenerator):
     aliases = ["MIN-Speck"]
 
     def rho(self, x):
-        M = self.QDelta.shape[0]
+        M = self.size
         return max(abs(
             np.linalg.eigvals(np.eye(M) - np.diag(x).dot(self.Q))))
 
-    def getQDelta(self, k=None):
-        x0 = 10 * np.ones(self.Q.shape[0])
+    def computeQDelta(self, k=None):
+        x0 = 10 * np.ones(self.size)
         d = spo.minimize(self.rho, x0, method='Nelder-Mead')
-        np.copyto(self.QDelta, np.linalg.inv(np.diag(d.x)))
-        return self.QDelta
+        return np.linalg.inv(np.diag(d.x))
 
 
 class FromTable(QDeltaGenerator):
@@ -34,19 +33,30 @@ def __init__(self, nNodes, nodeType, quadType, **kwargs):
         self.nNodes = nNodes
         self.nodeType = nodeType
         self.quadType = quadType
-        self.QDelta = np.zeros((nNodes, nNodes), dtype=float)
+
+    @property
+    def size(self):
+        return self.nNodes
 
-    def getQDelta(self, k=None):
+    def computeQDelta(self, k=None):
         name = self.__class__.__name__
         try:
             table = getattr(tables, name)
-            coeffs = table[self.nodeType][self.quadType][self.nNodes]
+            coeffs = table[self.nodeType][self.quadType][self.size]
+            coeffs = np.asarray(coeffs, dtype=float)
+            assert coeffs.ndim == 1
+            assert coeffs.size == self.size
         except KeyError:
             raise NotImplementedError(
-                "no {name} MIN coefficients for"
+                f"no {name} MIN coefficients for"
                 f"{self.nNodes} {self.nodeType} {self.quadType} nodes")
-        np.copyto(self.QDelta, np.diag(coeffs))
-        return self.QDelta
+        except AssertionError:
+            raise ValueError(
+                f"MIN coefficients stored for {name} are inconsistent : {coeffs}")
+        except:
+            raise ValueError(
+                f"could not convert {name} MIN coefficients to numpy array")
+        return np.diag(coeffs)
 
     def check(cls):
         try:
@@ -85,25 +95,28 @@ class MIN_SR_NS(QDeltaGenerator):
 
     def __init__(self, nodes, **kwargs):
         self.nodes = np.asarray(nodes)
-        M = self.nodes.size
-        self.QDelta = np.zeros((M, M), dtype=float)
 
-    def getQDelta(self, k=None):
-        np.copyto(self.QDelta, np.diag(self.nodes)/self.nodes.size)
-        return self.QDelta
+    @property
+    def size(self):
+        return self.nodes.size
+
+    def computeQDelta(self, k=None):
+        return np.diag(self.nodes)/self.size
 
 
 @register
 class MIN_SR_S(QDeltaGenerator):
     aliases = ["MIN-SR-S"]
 
     def __init__(self, nNodes, nodeType, quadType, **kwargs):
-        self.nNodes = nNodes
         self.nodeType = nodeType
         self.quadType = quadType
-        self.QDelta = np.zeros((nNodes, nNodes), dtype=float)
         self.coll = Collocation(nNodes, nodeType, quadType)
 
+    @property
+    def size(self):
+        return self.coll.nodes.size
+
     def computeCoeffs(self, M, a=None, b=None):
         """
         Compute diagonal coefficients for a given number of nodes M.
@@ -132,7 +145,7 @@ def computeCoeffs(self, M, a=None, b=None):
             The nodes associated to the current coefficients.
         """
         nodeType, quadType = self.nodeType, self.quadType
-        if M == self.nNodes:
+        if M == self.size:
             collM = self.coll
         else:
             collM = Collocation(nNodes=M, nodeType=nodeType, quadType=quadType)
@@ -182,35 +195,32 @@ def lawDiff(ab):
         sol = spo.minimize(lawDiff, [1.0, 1.0], method="nelder-mead")
         return sol.x
 
-    def getQDelta(self, k=None):
+    def computeQDelta(self, k=None):
         # Compute coefficients incrementally
         a, b = None, None
         m0 = 2 if self.quadType in ['LOBATTO', 'RADAU-LEFT'] else 1
-        for m in range(m0, self.nNodes + 1):
+        for m in range(m0, self.size + 1):
             coeffs, nodes = self.computeCoeffs(m, a, b)
             if m > 1:
                 a, b = self.fit(coeffs * m, nodes)
-
-        np.copyto(self.QDelta, np.diag(coeffs))
-        return self.QDelta
+        return np.diag(coeffs)
 
 
 @register
 class MIN_SR_FLEX(MIN_SR_S):
+    _K_DEP = True
     aliases = ["MIN-SR-FLEX"]
 
-    def getQDelta(self, k=None):
+    def computeQDelta(self, k=None):
         if k is None:
             k = 1
         if k < 1:
             raise ValueError(f"k must be greater than 0 ({k})")
-        nodes = self.coll.nodes
-        if k <= nodes.size:
-            np.copyto(self.QDelta, np.diag(nodes)/k)
+        if k <= self.size:
+            return np.diag(self.coll.nodes/k)
         else:
             try:
-                self.QDelta_MIN_SR_S
+                self._QDelta_MIN_SR_S
             except AttributeError:
-                self.QDelta_MIN_SR_S = super().getQDelta().copy()
-            np.copyto(self.QDelta, self.QDelta_MIN_SR_S)
-        return self.QDelta
+                self._QDelta_MIN_SR_S = super().computeQDelta()
+            return self._QDelta_MIN_SR_S

qmat/qdelta/timestepping.py

@@ -7,53 +7,60 @@
 
 from qmat.qdelta import QDeltaGenerator, register
 
+
 class TimeStepping(QDeltaGenerator):
 
     def __init__(self, nodes, tLeft=0, **kwargs):
         nodes = np.asarray(nodes)
+
         deltas = nodes.copy()
         deltas[0] = nodes[0] - tLeft
         deltas[1:] = np.ediff1d(nodes)
+
         self.deltas = deltas
-        M = nodes.size
-        self.QDelta = np.zeros((M, M), dtype=float)
         self.nodes = nodes
+        self.tLeft = tLeft
+
+    @property
+    def size(self):
+        return self.nodes.size
 
 
 @register
 class BE(TimeStepping):
     """Approximation based on Backward Euler steps between the nodes"""
     aliases = ["IE"]
 
-    def getQDelta(self, k=None):
-        QDelta, M, deltas = self.QDelta, self.nodes.size, self.deltas
+    def computeQDelta(self, k=None):
+        QDelta, M, deltas = self.zeros, self.nodes.size, self.deltas
         for i in range(M):
             QDelta[i:, :M-i] += np.diag(deltas[:M-i])
-        return self.storeAndReturn(QDelta)
+        return QDelta
 
 
 @register
 class FE(TimeStepping):
     """Approximation based on Forward Euler steps between the nodes"""
     aliases = ["EE"]
 
-    def getQDelta(self, k=None):
-        QDelta, M, deltas = self.QDelta, self.nodes.size, self.deltas
+    def computeQDelta(self, k=None):
+        QDelta, M, deltas = self.zeros, self.nodes.size, self.deltas
         for i in range(1, M):
             QDelta[i:, :M-i] += np.diag(deltas[1:M-i+1])
-        return self.storeAndReturn(QDelta)
+        return QDelta
 
     @property
     def dTau(self):
         return self.nodes*0 + self.deltas[0]
 
+
 @register
 class TRAP(TimeStepping):
     """Approximation based on Trapezoidal Rule between the nodes"""
     aliases = ["CN"]
 
-    def getQDelta(self, k=None):
-        QDelta, M, deltas = self.QDelta, self.nodes.size, self.deltas
+    def computeQDelta(self, k=None):
+        QDelta, M, deltas = self.zeros, self.nodes.size, self.deltas
         for i in range(0, M):
             QDelta[i:, :M-i] += np.diag(deltas[:M-i])
         for i in range(1, M):
@@ -71,19 +78,17 @@ class BEPAR(TimeStepping):
     """Approximation based on parallel Backward Euler steps from zero to nodes"""
     aliases = ["IEpar"]
 
-    def getQDelta(self, k=None):
-        self.QDelta[:] = np.diag(self.nodes)
-        return self.QDelta
+    def computeQDelta(self, k=None):
+        return np.diag(self.nodes) - self.tLeft
 
 
 @register
 class TRAPAR(TimeStepping):
     """Approximation based on parallel Trapezoidal Rule from zero to nodes"""
 
-    def getQDelta(self, k=None):
-        self.QDelta[:] = np.diag(self.nodes/2)
-        return self.QDelta
+    def computeQDelta(self, k=None):
+        return np.diag(self.nodes/2) - self.tLeft
 
     @property
     def dTau(self):
-        return self.nodes/2.0
+        return self.nodes/2.0 - self.tLeft

tests/test_qdelta/test_min.py

@@ -106,10 +106,6 @@ def testFlex(nNodes, nodeType, quadType):
     for k in range(nNodes):
         P = (I - np.linalg.solve(gen.getQDelta(k+1), Q)) @ P
 
-    assert np.allclose(np.tril(P), P), \
-        "QDelta product is not lower triangular"
-    assert np.allclose(P, np.diag(np.diag(P))), \
-        "QDelta product is not diagonal"
     assert np.linalg.norm(P, ord=np.inf) < 1e-13 * margin, \
         "nilpotency measure is to high"