evaluate coeffs lazily

Author: basnijholt

adaptive/learner/integrator_coeffs.py

@@ -141,39 +141,105 @@ def calc_V(x, n):
     return np.array(V).T
 
 
-eps = np.spacing(1)
-
-# the nodes and Newton polynomials
-ns = (5, 9, 17, 33)
-xi = [-np.cos(np.linspace(0, np.pi, n)) for n in ns]
-
-# Make `xi` perfectly anti-symmetric, important for splitting the intervals
-xi = [(row - row[::-1]) / 2 for row in xi]
-
-# Compute the Vandermonde-like matrix and its inverse.
-V = [calc_V(x, n) for x, n in zip(xi, ns)]
-V_inv = list(map(scipy.linalg.inv, V))
-Vcond = [scipy.linalg.norm(a, 2) * scipy.linalg.norm(b, 2) for a, b in zip(V, V_inv)]
-
-# Compute the shift matrices.
-T_left, T_right = [V_inv[3] @ calc_V((xi[3] + a) / 2, ns[3]) for a in [-1, 1]]
-
-# If the relative difference between two consecutive approximations is
-# lower than this value, the error estimate is considered reliable.
-# See section 6.2 of Pedro Gonnet's thesis.
-hint = 0.1
-
-# Smallest acceptable relative difference of points in a rule.  This was chosen
-# such that no artifacts are apparent in plots of (i, log(a_i)), where a_i is
-# the sequence of estimates of the integral value of an interval and all its
-# ancestors..
-min_sep = 16 * eps
-
-ndiv_max = 20
-
-# set-up the downdate matrix
-k = np.arange(ns[3])
-alpha = np.sqrt((k + 1) ** 2 / (2 * k + 1) / (2 * k + 3))
-gamma = np.concatenate([[0, 0], np.sqrt(k[2:] ** 2 / (4 * k[2:] ** 2 - 1))])
-
-b_def = calc_bdef(ns)
+class Coefficients:
+    def __init__(self) -> None:
+        self.is_set = False
+        # The nodes
+        self.ns = (5, 9, 17, 33)
+        self.eps = np.spacing(1)
+
+        # If the relative difference between two consecutive approximations is
+        # lower than this value, the error estimate is considered reliable.
+        # See section 6.2 of Pedro Gonnet's thesis.
+        self.hint = 0.1
+
+        # Smallest acceptable relative difference of points in a rule.  This was chosen
+        # such that no artifacts are apparent in plots of (i, log(a_i)), where a_i is
+        # the sequence of estimates of the integral value of an interval and all its
+        # ancestors..
+        self.min_sep = 16 * self.eps
+
+        # Maximum amount of subdivisions
+        self.ndiv_max = 20
+
+    def set(self):
+        self.is_set = True
+        # The Newton polynomials
+        xi = [-np.cos(np.linspace(0, np.pi, n)) for n in self.ns]
+        # Make `xi` perfectly anti-symmetric, important for splitting the intervals
+        self._xi = [(row - row[::-1]) / 2 for row in xi]
+
+        # Compute the Vandermonde-like matrix and its inverse.
+        self._V = [calc_V(x, n) for x, n in zip(xi, self.ns)]
+        self._V_inv = list(map(scipy.linalg.inv, self._V))
+        self._Vcond = [
+            scipy.linalg.norm(a, 2) * scipy.linalg.norm(b, 2)
+            for a, b in zip(self._V, self._V_inv)
+        ]
+
+        # Compute the shift matrices.
+        self._T_left, self._T_right = [
+            self._V_inv[3] @ calc_V((xi[3] + a) / 2, self.ns[3]) for a in [-1, 1]
+        ]
+
+        # set-up the downdate matrix
+        k = np.arange(self.ns[3])
+        self._alpha = np.sqrt((k + 1) ** 2 / (2 * k + 1) / (2 * k + 3))
+        self._gamma = np.concatenate(
+            [[0, 0], np.sqrt(k[2:] ** 2 / (4 * k[2:] ** 2 - 1))]
+        )
+        self._b_def = calc_bdef(self.ns)
+
+    @property
+    def xi(self):
+        if not self.is_set:
+            self.set()
+        return self._xi
+
+    @property
+    def V(self):
+        if not self.is_set:
+            self.set()
+        return self._V
+
+    @property
+    def V_inv(self):
+        if not self.is_set:
+            self.set()
+        return self._V_inv
+
+    @property
+    def Vcond(self):
+        if not self.is_set:
+            self.set()
+        return self._Vcond
+
+    @property
+    def T_right(self):
+        if not self.is_set:
+            self.set()
+        return self._T_right
+
+    @property
+    def T_left(self):
+        if not self.is_set:
+            self.set()
+        return self._T_left
+
+    @property
+    def alpha(self):
+        if not self.is_set:
+            self.set()
+        return self._alpha
+
+    @property
+    def gamma(self):
+        if not self.is_set:
+            self.set()
+        return self._gamma
+
+    @property
+    def b_def(self):
+        if not self.is_set:
+            self.set()
+        return self._b_def

adaptive/learner/integrator_learner.py

@@ -14,33 +14,21 @@
 from adaptive.notebook_integration import ensure_holoviews
 from adaptive.utils import cache_latest, restore
 
-from .integrator_coeffs import (
-    T_left,
-    T_right,
-    V_inv,
-    Vcond,
-    alpha,
-    b_def,
-    eps,
-    gamma,
-    hint,
-    min_sep,
-    ndiv_max,
-    ns,
-    xi,
-)
+from .integrator_coeffs import Coefficients
+
+c = Coefficients()
 
 
 def _downdate(c, nans, depth):
     # This is algorithm 5 from the thesis of Pedro Gonnet.
-    b = b_def[depth].copy()
-    m = ns[depth] - 1
+    b = c.b_def[depth].copy()
+    m = c.ns[depth] - 1
     for i in nans:
-        b[m + 1] /= alpha[m]
-        xii = xi[depth][i]
-        b[m] = (b[m] + xii * b[m + 1]) / alpha[m - 1]
+        b[m + 1] /= c.alpha[m]
+        xii = c.xi[depth][i]
+        b[m] = (b[m] + xii * b[m + 1]) / c.alpha[m - 1]
         for j in range(m - 1, 0, -1):
-            b[j] = (b[j] + xii * b[j + 1] - gamma[j + 1] * b[j + 2]) / alpha[j - 1]
+            b[j] = (b[j] + xii * b[j + 1] - c.gamma[j + 1] * b[j + 2]) / c.alpha[j - 1]
         b = b[1:]
 
         c[:m] -= c[m] / b[m] * b[:m]
@@ -62,7 +50,7 @@ def _zero_nans(fx):
 def _calc_coeffs(fx, depth):
     """Caution: this function modifies fx."""
     nans = _zero_nans(fx)
-    c_new = V_inv[depth] @ fx
+    c_new = c.V_inv[depth] @ fx
     if nans:
         fx[nans] = np.nan
         c_new = _downdate(c_new, nans, depth)
@@ -168,11 +156,11 @@ def T(self):
         left = self.a == self.parent.a
         right = self.b == self.parent.b
         assert left != right
-        return T_left if left else T_right
+        return c.T_left if left else c.T_right
 
     def refinement_complete(self, depth):
         """The interval has all the y-values to calculate the intergral."""
-        if len(self.data) < ns[depth]:
+        if len(self.data) < c.ns[depth]:
             return False
         return all(p in self.data for p in self.points(depth))
 
@@ -181,7 +169,7 @@ def points(self, depth=None):
             depth = self.depth
         a = self.a
         b = self.b
-        return (a + b) / 2 + (b - a) * xi[depth] / 2
+        return (a + b) / 2 + (b - a) * c.xi[depth] / 2
 
     def refine(self):
         self.depth += 1
@@ -234,7 +222,7 @@ def calc_ndiv(self):
         div = self.parent.c00 and self.c00 / self.parent.c00 > 2
         self.ndiv += div
 
-        if self.ndiv > ndiv_max and 2 * self.ndiv > self.rdepth:
+        if self.ndiv > c.ndiv_max and 2 * self.ndiv > self.rdepth:
             raise DivergentIntegralError
 
         if div:
@@ -243,7 +231,7 @@ def calc_ndiv(self):
 
     def update_ndiv_recursively(self):
         self.ndiv += 1
-        if self.ndiv > ndiv_max and 2 * self.ndiv > self.rdepth:
+        if self.ndiv > c.ndiv_max and 2 * self.ndiv > self.rdepth:
             raise DivergentIntegralError
 
         for child in self.children:
@@ -276,21 +264,21 @@ def complete_process(self, depth):
         if depth:
             # Refine
             c_diff = self.calc_err(c_old)
-            force_split = c_diff > hint * norm(self.c)
+            force_split = c_diff > c.hint * norm(self.c)
         else:
             # Split
             self.c00 = self.c[0]
 
             if self.parent.depth_complete is not None:
-                c_old = self.T[:, : ns[self.parent.depth_complete]] @ self.parent.c
+                c_old = self.T[:, : c.ns[self.parent.depth_complete]] @ self.parent.c
                 self.calc_err(c_old)
                 self.calc_ndiv()
 
             for child in self.children:
                 if child.depth_complete is not None:
                     child.calc_ndiv()
                 if child.depth_complete == 0:
-                    c_old = child.T[:, : ns[self.depth_complete]] @ self.c
+                    c_old = child.T[:, : c.ns[self.depth_complete]] @ self.c
                     child.calc_err(c_old)
 
         if self.done_leaves is not None and not len(self.done_leaves):
@@ -320,7 +308,7 @@ def complete_process(self, depth):
                 ival.done_leaves -= old_leaves
                 ival = ival.parent
 
-        remove = self.err < (abs(self.igral) * eps * Vcond[depth])
+        remove = self.err < (abs(self.igral) * c.eps * c.Vcond[depth])
 
         return force_split, remove
 
@@ -496,8 +484,8 @@ def _fill_stack(self):
         points = ival.points()
 
         if (
-            points[1] - points[0] < points[0] * min_sep
-            or points[-1] - points[-2] < points[-2] * min_sep
+            points[1] - points[0] < points[0] * c.min_sep
+            or points[-1] - points[-2] < points[-2] * c.min_sep
         ):
             self.ivals.remove(ival)
         elif ival.depth == 3 or force_split:

adaptive/tests/test_cquad.py

@@ -5,14 +5,16 @@
 import pytest
 
 from adaptive.learner import IntegratorLearner
-from adaptive.learner.integrator_coeffs import ns
+from adaptive.learner.integrator_coeffs import Coefficients
 from adaptive.learner.integrator_learner import DivergentIntegralError
 
 from .algorithm_4 import DivergentIntegralError as A4DivergentIntegralError
 from .algorithm_4 import algorithm_4, f0, f7, f21, f24, f63, fdiv
 
 eps = np.spacing(1)
 
+ns = Coefficients().ns
+
 
 def run_integrator_learner(f, a, b, tol, n):
     learner = IntegratorLearner(f, bounds=(a, b), tol=tol)