Merge pull request #311 from python-adaptive/coeff-cache

lazily evaluate the integrator coefficients

Author: basnijholt

adaptive/learner/integrator_coeffs.py

@@ -2,6 +2,7 @@
 
 from collections import defaultdict
 from fractions import Fraction
+from functools import lru_cache
 
 import numpy as np
 import scipy.linalg
@@ -141,39 +142,49 @@ def calc_V(x, n):
     return np.array(V).T
 
 
-eps = np.spacing(1)
+@lru_cache(maxsize=None)
+def _coefficients():
+    """Compute the coefficients on demand, in order to avoid doing linear algebra on import."""
+    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]
+    # 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]
+    # 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 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]]
+    # 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
+    # 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
+    # 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
+    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))])
+    # 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)
+    b_def = calc_bdef(ns)
+    return locals()
+
+
+def __getattr__(attr):
+    return _coefficients()[attr]

adaptive/learner/integrator_learner.py

@@ -10,37 +10,24 @@
 from scipy.linalg import norm
 from sortedcontainers import SortedSet
 
+import adaptive.learner.integrator_coeffs as coeff
 from adaptive.learner.base_learner import BaseLearner
 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,
-)
-
 
 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 = coeff.b_def[depth].copy()
+    m = coeff.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] /= coeff.alpha[m]
+        xii = coeff.xi[depth][i]
+        b[m] = (b[m] + xii * b[m + 1]) / coeff.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] - coeff.gamma[j + 1] * b[j + 2]
+            ) / coeff.alpha[j - 1]
         b = b[1:]
 
         c[:m] -= c[m] / b[m] * b[:m]
@@ -62,7 +49,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 = coeff.V_inv[depth] @ fx
     if nans:
         fx[nans] = np.nan
         c_new = _downdate(c_new, nans, depth)
@@ -168,11 +155,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 coeff.T_left if left else coeff.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) < coeff.ns[depth]:
             return False
         return all(p in self.data for p in self.points(depth))
 
@@ -181,7 +168,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) * coeff.xi[depth] / 2
 
     def refine(self):
         self.depth += 1
@@ -234,7 +221,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 > coeff.ndiv_max and 2 * self.ndiv > self.rdepth:
             raise DivergentIntegralError
 
         if div:
@@ -243,7 +230,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 > coeff.ndiv_max and 2 * self.ndiv > self.rdepth:
             raise DivergentIntegralError
 
         for child in self.children:
@@ -276,21 +263,23 @@ 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 > coeff.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[:, : coeff.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[:, : coeff.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 +309,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) * coeff.eps * coeff.Vcond[depth])
 
         return force_split, remove
 
@@ -496,8 +485,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] * coeff.min_sep
+            or points[-1] - points[-2] < points[-2] * coeff.min_sep
         ):
             self.ivals.remove(ival)
         elif ival.depth == 3 or force_split:

adaptive/tests/test_cquad.py

@@ -4,8 +4,8 @@
 import numpy as np
 import pytest
 
+import adaptive.learner.integrator_coeffs as coeff
 from adaptive.learner import IntegratorLearner
-from adaptive.learner.integrator_coeffs import ns
 from adaptive.learner.integrator_learner import DivergentIntegralError
 
 from .algorithm_4 import DivergentIntegralError as A4DivergentIntegralError
@@ -245,7 +245,7 @@ def test_removed_choose_mutiple_points_at_once():
     we should use the high-precision interval.
     """
     learner = IntegratorLearner(np.exp, bounds=(0, 1), tol=1e-15)
-    n = ns[-1] + 2 * (ns[0] - 2)  # first + two children (33+6=39)
+    n = coeff.ns[-1] + 2 * (coeff.ns[0] - 2)  # first + two children (33+6=39)
     xs, _ = learner.ask(n)
     for x in xs:
         learner.tell(x, learner.function(x))
@@ -257,7 +257,7 @@ def test_removed_ask_one_by_one():
         # This test should raise because integrating np.exp should be done
         # after the 33th point
         learner = IntegratorLearner(np.exp, bounds=(0, 1), tol=1e-15)
-        n = ns[-1] + 2 * (ns[0] - 2)  # first + two children (33+6=39)
+        n = coeff.ns[-1] + 2 * (coeff.ns[0] - 2)  # first + two children (33+6=39)
         for _ in range(n):
             xs, _ = learner.ask(1)
             for x in xs: