turn split and refine into coroutines and move them to _Interval

1 files changed, 62 insertions, 49 deletions

diff --git a/vquad/core.py b/vquad/core.py
index e7faef2..2f3cac0 100644
--- a/vquad/core.py
+++ b/vquad/core.py
@@ -134,6 +134,64 @@ class _Interval:
         x = (2 * x - (a + b)) / (b - a)
         return _eval_legendre(self.coeffs, x)
 
+    def split(self):
+        """Split this interval in the center into two children.
+
+        This is a coroutine that initially yields an array of x values
+        of points to be evaluated.  Once the corresponding values have
+        been sent back a tuple containing the child intervals is
+        yielded and execution ends.
+        """
+        m = (self.a + self.b) / 2
+        f_center = self.vals[(len(self.vals) - 1) // 2]
+
+        depth = self.depth + 1
+        children = [_Interval(self.a, m, min_level, depth),
+                    _Interval(m, self.b, min_level, depth)]
+        points = np.concatenate([child.points()[1:-1] for child in children])
+        valss = np.empty((2, tbls.sizes[min_level]))
+        valss[:, 0] = self.vals[0], f_center
+        valss[:, -1] = f_center, self.vals[-1]
+        valss[:, 1:-1] = (yield points).reshape((2, -1))
+
+        for child, vals, T in zip(children, valss, tbls.Ts):
+            child.interpolate(vals, T[:, :self.coeffs.shape[0]] @ self.coeffs)
+            child.ndiv = (self.ndiv
+                         + (self.c00 and child.c00 / self.c00 > 2))
+            if child.ndiv > ndiv_max and 2*child.ndiv > child.depth:
+                msg = ('Possibly divergent integral in the interval '
+                       '[{}, {}]! (h={})')
+                raise DivergentIntegralError(
+                    msg.format(child.a, child.b, child.b - child.a),
+                    child.igral * np.inf, None)
+        yield children
+
+    def refine(self):
+        """Increase degree of interval.
+
+        This is a coroutine that initially yields an array of x values
+        of points to be evaluated.  Once the corresponding values have
+        been sent back, a bool is yielded and execution ends.
+
+        It is "true" if further refinements/splits of the interval seem
+        promising, and "false" otherwise.  This is the case when
+        neigboring points can be resolved only barely by floating point
+        numbers, or when the estimated relative error is already at the
+        limit of numerical accuracy and cannot be reduced further.
+        """
+        self.level += 1
+        points = self.points()
+        vals = np.empty(points.shape)
+        vals[0::2] = self.vals
+        vals[1::2] = (yield points[1::2])
+        self.interpolate(vals, self.coeffs)
+
+        yield (points[1] - points[0] > points[0] * min_sep
+               and points[-1] - points[-2] > points[-2] * min_sep
+               and self.err > (abs(self.igral)
+                               * eps * tbls.V_cond_nums[self.level]))
+
+
 class Vquad:
     """Evaluate an integral using adaptive quadrature.
 
@@ -169,60 +227,14 @@ class Vquad:
         ival.next = self.end = _Terminator()
         self.end.prev = ival
 
-    def split(self, ival):
-        m = (ival.a + ival.b) / 2
-        f_center = ival.vals[(len(ival.vals) - 1) // 2]
-
-        depth = ival.depth + 1
-        children = [_Interval(ival.a, m, min_level, depth),
-                    _Interval(m, ival.b, min_level, depth)]
-        points = np.concatenate([child.points()[1:-1] for child in children])
-        valss = np.empty((2, tbls.sizes[min_level]))
-        valss[:, 0] = ival.vals[0], f_center
-        valss[:, -1] = f_center, ival.vals[-1]
-        valss[:, 1:-1] = self.f(points).reshape((2, -1))
-
-        for child, vals, T in zip(children, valss, tbls.Ts):
-            child.interpolate(vals, T[:, :ival.coeffs.shape[0]] @ ival.coeffs)
-            child.ndiv = (ival.ndiv
-                         + (ival.c00 and child.c00 / ival.c00 > 2))
-            if child.ndiv > ndiv_max and 2*child.ndiv > child.depth:
-                msg = ('Possibly divergent integral in the interval '
-                       '[{}, {}]! (h={})')
-                raise DivergentIntegralError(
-                    msg.format(child.a, child.b, child.b - child.a),
-                    child.igral * np.inf, None)
-        return children
-
-    def refine(self, ival):
-        """Increase degree of interval.
-
-        Returns True if the refined interval is OK, and False if it is
-        borderline and should not be refined/split further.  This
-        happens when neigboring points can be only barely resolved by
-        floating point numbers, or when the estimated relative error is
-        already at the limit of numerical accuracy and cannot be reduced
-        further.
-        """
-        ival.level += 1
-        points = ival.points()
-        vals = np.empty(points.shape)
-        vals[0::2] = ival.vals
-        vals[1::2] = self.f(points[1::2])
-        ival.interpolate(vals, ival.coeffs)
-
-        return (points[1] - points[0] > points[0] * min_sep
-                and points[-1] - points[-2] > points[-2] * min_sep
-                and ival.err > (abs(ival.igral)
-                                * eps * tbls.V_cond_nums[ival.level]))
-
     def improve(self):
         ival = self.ivals[-1]
 
         if ival.level == max_level:
             split = True
         else:
-            if not self.refine(ival):
+            refine = ival.refine()
+            if not refine.send(self.f(next(refine))):
                 # Remove the interval but remember the excess integral and
                 # error.
                 self.err_excess += ival.err
@@ -234,7 +246,8 @@ class Vquad:
         if split:
             # Replace current interval by its children.
             self.ivals.pop()
-            child0, child1 = self.split(ival)
+            split = ival.split()
+            child0, child1 = split.send(self.f(next(split)))
             bisect.insort(self.ivals, child0)
             bisect.insort(self.ivals, child1)