.gitlab-ci.yml

 image: quantumtinkerer/research
+
 test:
   script:
     - pip install -r test-requirements.txt
@@ -7,3 +8,9 @@ test:
   artifacts:
     paths:
       - htmlcov
+
+authors check:
+  script:
+    - MISSING_AUTHORS=$(git shortlog -s HEAD | sed -e "s/^[0-9\t ]*//"| xargs -i sh -c 'grep -q "{}" AUTHORS.md || echo "{} missing from authors"')
+    - if [ ! -z "$MISSING_AUTHORS" ]; then { echo $MISSING_AUTHORS; exit 1; }; fi
+  allow_failure: true

RELEASE.md

+# Making a Adaptive release
+
+This document guides a contributor through creating a release of Adaptive.
+
+
+## Preflight checks
+
+The following checks should be made *before* tagging the release.
+
+
+#### Check that all issues are resolved
+
+Check that all the issues and merge requests for the appropriate
+[milestone](https://gitlab.kwant-project.org/qt/adaptive/issues)
+have been resolved. Any unresolved issues should have their milestone
+bumped.
+
+
+#### Ensure that all tests pass
+
+For major and minor releases we will be tagging the ``master`` branch.
+This should be as simple as verifying that the 
+[latest CI pipeline](https://gitlab.kwant-project.org/qt/adaptive/pipelines) 
+succeeded.
+
+
+#### Verify that `AUTHORS.md` is up-to-date
+
+The following command shows the number of commits per author since the last
+annotated tag:
+```
+t=$(git describe --abbrev=0); echo Commits since $t; git shortlog -s $t..
+```
+
+## Make a release, but do not publish it yet
+
+
+### Tag the release
+
+Make an **annotated, signed** tag for the release. The tag must have the name:
+```
+git tag -s v<version> -m "version <version>"
+```
+
+
+### Build a source tarball and wheels and test it
+
+```
+rm -fr build dist
+python setup.py sdist bdist_wheel
+```
+
+This creates the file `dist/adaptive-<version>.tar.gz`.  It is a good idea to unpack it 
+and check that the tests run:
+```
+tar xzf dist/adaptive*.tar.gz
+cd adaptive-*
+py.test .
+```
+
+## Upload to PyPI
+
+```
+twine upload dist/*
+```
+
+
+## Update the [conda-forge recipe](https://github.com/conda-forge/adaptive-feedstock)
+
+* Fork the [feedstock repo](https://github.com/conda-forge/adaptive-feedstock)
+* Change the version number and sha256 in `recipe/meta.yaml` and commit to your fork
+* Open a [Pull Request](https://github.com/conda-forge/adaptive-feedstock/compare)
+* Type `@conda-forge-admin, please rerender` as a comment
+* When the tests succeed, merge

adaptive/learner/average_learner.py

@@ -14,12 +14,19 @@ class AverageLearner(BaseLearner):
     The learned function must depend on an integer input variable that
     represents the source of randomness.
 
-    Parameters:
-    -----------
+    Parameters
+    ----------
     atol : float
         Desired absolute tolerance
     rtol : float
         Desired relative tolerance
+
+    Attributes
+    ----------
+    data : dict
+        Sampled points and values.
+    pending_points : set
+        Points that still have to be evaluated.
     """
 
     def __init__(self, function, atol=None, rtol=None):
@@ -31,6 +38,7 @@ class AverageLearner(BaseLearner):
             rtol = np.inf
 
         self.data = {}
+        self.pending_points = set()
         self.function = function
         self.atol = atol
         self.rtol = rtol
@@ -40,28 +48,36 @@ class AverageLearner(BaseLearner):
 
     @property
     def n_requested(self):
-        return len(self.data)
+        return len(self.data) + len(self.pending_points)
 
-    def ask(self, n, add_data=True):
+    def ask(self, n, tell_pending=True):
         points = list(range(self.n_requested, self.n_requested + n))
+
+        if any(p in self.data or p in self.pending_points for p in points):
+            # This means some of the points `< self.n_requested` do not exist.
+            points = list(set(range(self.n_requested + n))
+                          - set(self.data)
+                          - set(self.pending_points))[:n]
+
         loss_improvements = [self.loss_improvement(n) / n] * n
-        if add_data:
-            self.tell_many(points, itertools.repeat(None))
+        if tell_pending:
+            for p in points:
+                self.tell_pending(p)
         return points, loss_improvements
 
     def tell(self, n, value):
-        value_is_new = not (n in self.data and value == self.data[n])
-        if not value_is_new:
-            value_old = self.data[n]
+        if n in self.data:
+            # The point has already been added before.
+            return
+
         self.data[n] = value
-        if value is not None:
-            self.sum_f += value
-            self.sum_f_sq += value**2
-            if value_is_new:
-                self.npoints += 1
-            else:
-                self.sum_f -= value_old
-                self.sum_f_sq -= value_old**2
+        self.pending_points.discard(n)
+        self.sum_f += value
+        self.sum_f_sq += value**2
+        self.npoints += 1
+
+    def tell_pending(self, n):
+        self.pending_points.add(n)
 
     @property
     def mean(self):
@@ -94,7 +110,7 @@ class AverageLearner(BaseLearner):
 
     def remove_unfinished(self):
         """Remove uncomputed data from the learner."""
-        pass
+        self.pending_points = set()
 
     def plot(self):
         hv = ensure_holoviews()

adaptive/learner/balancing_learner.py

@@ -64,25 +64,28 @@ class BalancingLearner(BaseLearner):
 
     def _ask_and_tell(self, n):
         points = []
+        loss_improvements = []
         for _ in range(n):
-            loss_improvements = []
+            improvements_per_learner = []
             pairs = []
             for index, learner in enumerate(self.learners):
                 if index not in self._points:
                     self._points[index] = learner.ask(
-                        n=1, add_data=False)
+                        n=1, tell_pending=False)
                 point, loss_improvement = self._points[index]
-                loss_improvements.append(loss_improvement[0])
+                improvements_per_learner.append(loss_improvement[0])
                 pairs.append((index, point[0]))
-            x, _ = max(zip(pairs, loss_improvements), key=itemgetter(1))
+            x, l = max(zip(pairs, improvements_per_learner),
+                       key=itemgetter(1))
             points.append(x)
-            self.tell(x, None)
+            loss_improvements.append(l)
+            self.tell_pending(x)
 
-        return points, None
+        return points, loss_improvements
 
-    def ask(self, n, add_data=True):
+    def ask(self, n, tell_pending=True):
         """Chose points for learners."""
-        if not add_data:
+        if not tell_pending:
             with restore(*self.learners):
                 return self._ask_and_tell(n)
         else:
@@ -94,6 +97,12 @@ class BalancingLearner(BaseLearner):
         self._loss.pop(index, None)
         self.learners[index].tell(x, y)
 
+    def tell_pending(self, x):
+        index, x = x
+        self._points.pop(index, None)
+        self._loss.pop(index, None)
+        self.learners[index].tell_pending(x)
+
     def loss(self, real=True):
         losses = []
         for index, learner in enumerate(self.learners):

adaptive/learner/base_learner.py

@@ -44,6 +44,12 @@ class BaseLearner(metaclass=abc.ABCMeta):
         for x, y in zip(xs, ys):
             self.tell(x, y)
 
+    @abc.abstractmethod
+    def tell_pending(self, x):
+        """Tell the learner that 'x' has been requested such
+        that it's not suggested again."""
+        pass
+
     @abc.abstractmethod
     def remove_unfinished(self):
         """Remove uncomputed data from the learner."""
@@ -62,14 +68,14 @@ class BaseLearner(metaclass=abc.ABCMeta):
         """
 
     @abc.abstractmethod
-    def ask(self, n, add_data=True):
+    def ask(self, n, tell_pending=True):
         """Choose the next 'n' points to evaluate.
 
         Parameters
         ----------
         n : int
             The number of points to choose.
-        add_data : bool, default: True
+        tell_pending : bool, default: True
             If True, add the chosen points to this
             learner's 'data' with 'None' for the 'y'
             values. Set this to False if you do not

adaptive/learner/data_saver.py

@@ -32,10 +32,13 @@ class DataSaver:
         return getattr(self.learner, attr)
 
     def tell(self, x, result):
-        y = self.arg_picker(result) if result is not None else None
+        y = self.arg_picker(result)
         self.extra_data[x] = result
         self.learner.tell(x, y)
 
+    def tell_pending(self, x):
+        self.learner.tell_pending(x)
+
 
 def _ds(learner_type, arg_picker, *args, **kwargs):
     args = args[2:]  # functools.partial passes the first 2 arguments in 'args'!

adaptive/learner/integrator_learner.py

@@ -15,6 +15,7 @@ from .integrator_coeffs import (b_def, T_left, T_right, ns, hint,
                                 ndiv_max, min_sep, eps, xi, V_inv,
                                 Vcond, alpha, gamma)
 from ..notebook_integration import ensure_holoviews
+from ..utils import restore
 
 
 def _downdate(c, nans, depth):
@@ -393,6 +394,9 @@ class IntegratorLearner(BaseLearner):
                         assert ival in self.ivals
                         self.priority_split.append(ival)
 
+    def tell_pending(self):
+        pass
+
     def propagate_removed(self, ival):
         def _propagate_removed_down(ival):
             ival.removed = True
@@ -414,10 +418,16 @@ class IntegratorLearner(BaseLearner):
                 self._stack.append(x)
         self.ivals.add(ival)
 
-    def ask(self, n, add_data=True):
-        if not add_data:
-            raise NotImplementedError(
-                "Asking points irreversibly changes the learner's data structure.")
+
+    def ask(self, n, tell_pending=True):
+        """Choose points for learners."""
+        if not tell_pending:
+            with restore(self):
+                return self._ask_and_tell_pending(n)
+        else:
+            return self._ask_and_tell_pending(n)
+
+    def _ask_and_tell_pending(self, n):
         points, loss_improvements = self.pop_from_stack(n)
         n_left = n - len(points)
         while n_left > 0:

adaptive/learner/learner1D.py

@@ -55,6 +55,18 @@ def default_loss(interval, scale, function_values):
     return loss
 
 
+def linspace(x_left, x_right, n):
+    """This is equivalent to
+    'np.linspace(x_left, x_right, n, endpoint=False)[1:]',
+    but it is 15-30 times faster for small 'n'."""
+    if n == 1:
+        # This is just an optimization
+        return []
+    else:
+        step = (x_right - x_left) / n
+        return [x_left + step * i for i in range(1, n)]
+
+
 class Learner1D(BaseLearner):
     """Learns and predicts a function 'f:ℝ → ℝ^N'.
 
@@ -137,44 +149,64 @@ class Learner1D(BaseLearner):
                                               self._scale, self.data)
             self.losses[x_left, x_right] = loss
 
-            start = self.neighbors_combined.bisect_right(x_left)
-            end = self.neighbors_combined.bisect_left(x_right)
-            for i in range(start, end):
-                a, b = (self.neighbors_combined.iloc[i],
-                        self.neighbors_combined.iloc[i + 1])
+            # Iterate over all interpolated intervals in between
+            # x_left and x_right and set the newly interpolated loss.
+            a, b = x_left, None
+            while b != x_right:
+                b = self.neighbors_combined[a][1]
                 self.losses_combined[a, b] = (b - a) * loss / dx
-            if start == end:
-                self.losses_combined[x_left, x_right] = loss
+                a = b
 
     def update_losses(self, x, real=True):
+        # When we add a new point x, we should update the losses
+        # (x_left, x_right) are the "real" neighbors of 'x'.
+        x_left, x_right = self.find_neighbors(x, self.neighbors)
+        # (a, b) are the neighbors of the combined interpolated
+        # and "real" intervals.
+        a, b = self.find_neighbors(x, self.neighbors_combined)
+
+        # (a, b) is splitted into (a, x) and (x, b) so if (a, b) exists
+        self.losses_combined.pop((a, b), None)  # we get rid of (a, b).
+
         if real:
-            x_left, x_right = self.find_neighbors(x, self.neighbors)
+            # We need to update all interpolated losses in the interval
+            # (x_left, x) and (x, x_right). Since the addition of the point
+            # 'x' could change their loss.
             self.update_interpolated_loss_in_interval(x_left, x)
             self.update_interpolated_loss_in_interval(x, x_right)
-            self.losses.pop((x_left, x_right), None)
-        else:
-            losses_combined = self.losses_combined
-            x_left, x_right = self.find_neighbors(x, self.neighbors)
-            a, b = self.find_neighbors(x, self.neighbors_combined)
-            if x_left is not None and x_right is not None:
-                dx = x_right - x_left
-                loss = self.losses[x_left, x_right]
-                losses_combined[a, x] = (x - a) * loss / dx
-                losses_combined[x, b] = (b - x) * loss / dx
-            else:
-                if a is not None:
-                    losses_combined[a, x] = float('inf')
-                if b is not None:
-                    losses_combined[x, b] = float('inf')
-
-            losses_combined.pop((a, b), None)
 
-    def find_neighbors(self, x, neighbors):
+            # Since 'x' is in between (x_left, x_right),
+            # we get rid of the interval.
+            self.losses.pop((x_left, x_right), None)
+            self.losses_combined.pop((x_left, x_right), None)
+        elif x_left is not None and x_right is not None:
+            # 'x' happens to be in between two real points,
+            # so we can interpolate the losses.
+            dx = x_right - x_left
+            loss = self.losses[x_left, x_right]
+            self.losses_combined[a, x] = (x - a) * loss / dx
+            self.losses_combined[x, b] = (b - x) * loss / dx
+
+        # (no real point left of x) or (no real point right of a)
+        left_loss_is_unknown = ((x_left is None) or
+                                (not real and x_right is None))
+        if (a is not None) and left_loss_is_unknown:
+            self.losses_combined[a, x] = float('inf')
+
+        # (no real point right of x) or (no real point left of b)
+        right_loss_is_unknown = ((x_right is None) or
+                                 (not real and x_left is None))
+        if (b is not None) and right_loss_is_unknown:
+            self.losses_combined[x, b] = float('inf')
+
+    @staticmethod
+    def find_neighbors(x, neighbors):
         if x in neighbors:
             return neighbors[x]
         pos = neighbors.bisect_left(x)
-        x_left = neighbors.iloc[pos - 1] if pos != 0 else None
-        x_right = neighbors.iloc[pos] if pos != len(neighbors) else None
+        keys = neighbors.keys()
+        x_left = keys[pos - 1] if pos != 0 else None
+        x_right = keys[pos] if pos != len(neighbors) else None
         return x_left, x_right
 
     def update_neighbors(self, x, neighbors):
@@ -192,7 +224,7 @@ class Learner1D(BaseLearner):
 
         When the function returns a vector the learners y-scale is set by
         the level with the the largest peak-to-peak value.
-         """
+        """
         self._bbox[0][0] = min(self._bbox[0][0], x)
         self._bbox[0][1] = max(self._bbox[0][1], x)
         self._scale[0] = self._bbox[0][1] - self._bbox[0][0]
@@ -212,33 +244,33 @@ class Learner1D(BaseLearner):
                 self._scale[1] = self._bbox[1][1] - self._bbox[1][0]
 
     def tell(self, x, y):
-        real = y is not None
+        if x in self.data:
+            # The point is already evaluated before
+            return
 
-        if real:
-            # Add point to the real data dict
-            self.data[x] = y
-            # remove from set of pending points
-            self.pending_points.discard(x)
+        # either it is a float/int, if not, try casting to a np.array
+        if not isinstance(y, (float, int)):
+            y = np.asarray(y, dtype=float)
 
-            if self._vdim is None:
-                try:
-                    self._vdim = len(np.squeeze(y))
-                except TypeError:
-                    self._vdim = 1
-        else:
-            # The keys of pending_points are the unknown points
-            self.pending_points.add(x)
+        # Add point to the real data dict
+        self.data[x] = y
 
-        # Update the neighbors
-        self.update_neighbors(x, self.neighbors_combined)
-        if real:
-            self.update_neighbors(x, self.neighbors)
+        # remove from set of pending points
+        self.pending_points.discard(x)
 
-        # Update the scale
-        self.update_scale(x, y)
+        if self._vdim is None:
+            try:
+                self._vdim = len(np.squeeze(y))
+            except TypeError:
+                self._vdim = 1
 
-        # Update the losses
-        self.update_losses(x, real)
+        if not self.bounds[0] <= x <= self.bounds[1]:
+            return
+
+        self.update_neighbors(x, self.neighbors_combined)
+        self.update_neighbors(x, self.neighbors)
+        self.update_scale(x, y)
+        self.update_losses(x, real=True)
 
         # If the scale has increased enough, recompute all losses.
         if self._scale[1] > self._oldscale[1] * 2:
@@ -248,15 +280,35 @@ class Learner1D(BaseLearner):
 
             self._oldscale = deepcopy(self._scale)
 
-    def ask(self, n, add_data=True):
+    def tell_pending(self, x):
+        if x in self.data:
+            # The point is already evaluated before
+            return
+        self.pending_points.add(x)
+        self.update_neighbors(x, self.neighbors_combined)
+        self.update_losses(x, real=False)
+
+    def ask(self, n, tell_pending=True):
         """Return n points that are expected to maximally reduce the loss."""
+        points, loss_improvements = self._ask_points_without_adding(n)
+
+        if tell_pending:
+            for p in points:
+                self.tell_pending(p)
+
+        return points, loss_improvements
+
+    def _ask_points_without_adding(self, n):
+        """Return n points that are expected to maximally reduce the loss.
+        Without altering the state of the learner"""
         # Find out how to divide the n points over the intervals
         # by finding  positive integer n_i that minimize max(L_i / n_i) subject
         # to a constraint that sum(n_i) = n + N, with N the total number of
         # intervals.
-
         # Return equally spaced points within each interval to which points
         # will be added.
+
+        # XXX: when is this used and could we safely remove it without impacting performance?
         if n == 0:
             return [], []
 
@@ -264,43 +316,59 @@ class Learner1D(BaseLearner):
         missing_bounds = [b for b in self.bounds if b not in self.data
                           and b not in self.pending_points]
 
-        if missing_bounds:
-            loss_improvements = [np.inf] * n
-            # XXX: should check if points are present in self.data or self.pending_points
-            points = np.linspace(*self.bounds, n + 2 - len(missing_bounds)).tolist()
-            if len(missing_bounds) == 1:
-                points = points[1:] if missing_bounds[0] == self.bounds[1] else points[:-1]
-        else:
-            def xs(x_left, x_right, n):
-                if n == 1:
-                    # This is just an optimization
-                    return []
-                else:
-                    step = (x_right - x_left) / n
-                    return [x_left + step * i for i in range(1, n)]
-
-            # Calculate how many points belong to each interval.
-            x_scale = self._scale[0]
-            quals = [((-loss if not math.isinf(loss) else -(x[1] - x[0]) / x_scale, x, 1))
-                     for x, loss in self.losses_combined.items()]
-            heapq.heapify(quals)
-
-            for point_number in range(n):
-                quality, x, n = quals[0]
-                if abs(x[1] - x[0]) / (n + 1) <= self._dx_eps:
-                    # The interval is too small and should not be subdivided
-                    quality = np.inf
-                heapq.heapreplace(quals, (quality * n / (n + 1), x, n + 1))
-
-            points = list(itertools.chain.from_iterable(
-                xs(*x, n) for quality, x, n in quals))
-
-            loss_improvements = list(itertools.chain.from_iterable(
-                                     itertools.repeat(-quality, n - 1)
-                                     for quality, x, n in quals))
-
-        if add_data:
-            self.tell_many(points, itertools.repeat(None))
+        if len(missing_bounds) >= n:
+            return missing_bounds[:n], [np.inf] * n
+
+        def finite_loss(loss, xs):
+            # If the loss is infinite we return the
+            # distance between the two points.
+            return (loss if not math.isinf(loss)
+                else (xs[1] - xs[0]) / self._scale[0])
+
+        quals = [(-finite_loss(loss, x), x, 1)
+                 for x, loss in self.losses_combined.items()]
+
+        # Add bound intervals to quals if bounds were missing.
+        if len(self.data) + len(self.pending_points) == 0:
+            # We don't have any points, so return a linspace with 'n' points.
+            return np.linspace(*self.bounds, n).tolist(), [np.inf] * n
+        elif len(missing_bounds) > 0:
+            # There is at least one point in between the bounds.
+            all_points = list(self.data.keys()) + list(self.pending_points)
+            intervals = [(self.bounds[0], min(all_points)),
+                         (max(all_points), self.bounds[1])]
+            for interval, bound in zip(intervals, self.bounds):
+                if bound in missing_bounds:
+                    qual = (-finite_loss(np.inf, interval), interval, 1)
+                    quals.append(qual)
+
+        # Calculate how many points belong to each interval.
+        points, loss_improvements = self._subdivide_quals(
+            quals, n - len(missing_bounds))
+
+        points = missing_bounds + points
+        loss_improvements = [np.inf] * len(missing_bounds) + loss_improvements
+
+        return points, loss_improvements
+
+    def _subdivide_quals(self, quals, n):
+        # Calculate how many points belong to each interval.
+        heapq.heapify(quals)
+
+        for _ in range(n):
+            quality, x, n = quals[0]
+            if abs(x[1] - x[0]) / (n + 1) <= self._dx_eps:
+                # The interval is too small and should not be subdivided.
+                quality = np.inf
+                # XXX: see https://gitlab.kwant-project.org/qt/adaptive/issues/104
+            heapq.heapreplace(quals, (quality * n / (n + 1), x, n + 1))
+
+        points = list(itertools.chain.from_iterable(
+            linspace(*interval, n) for quality, interval, n in quals))
+
+        loss_improvements = list(itertools.chain.from_iterable(
+            itertools.repeat(-quality, n - 1)
+            for quality, interval, n in quals))
 
         return points, loss_improvements
 

adaptive/learner/learner2D.py

@@ -166,6 +166,9 @@ class Learner2D(BaseLearner):
     ----------
     data : dict
         Sampled points and values.
+    pending_points : set
+        Points that still have to be evaluated and are currently
+        interpolated, see `data_combined`.
     stack_size : int, default 10
         The size of the new candidate points stack. Set it to 1
         to recalculate the best points at each call to `ask`.
@@ -180,7 +183,7 @@ class Learner2D(BaseLearner):
     -------
     data_combined : dict
         Sampled points and values so far including
-        the unknown interpolated ones.
+        the unknown interpolated points in `pending_points`.
 
     Notes
     -----
@@ -217,7 +220,7 @@ class Learner2D(BaseLearner):
         self.bounds = tuple((float(a), float(b)) for a, b in bounds)
         self.data = OrderedDict()
         self._stack = OrderedDict()
-        self._interp = set()
+        self.pending_points = set()
 
         self.xy_mean = np.mean(self.bounds, axis=1)
         self._xy_scale = np.ptp(self.bounds, axis=1)
@@ -263,57 +266,83 @@ class Learner2D(BaseLearner):
 
     @property
     def bounds_are_done(self):
-        return not any((p in self._interp or p in self._stack)
+        return not any((p in self.pending_points or p in self._stack)
                        for p in self._bounds_points)
 
-    def data_combined(self):
-        # Interpolate the unfinished points
-        data_combined = copy(self.data)
-        if self._interp:
-            points_interp = list(self._interp)
+    def _data_in_bounds(self):
+        if self.data:
+            points = np.array(list(self.data.keys()))
+            values = np.array(list(self.data.values()), dtype=float)
+            ll, ur = np.reshape(self.bounds, (2, 2)).T
+            inds = np.all(np.logical_and(ll <= points, points <= ur), axis=1)
+            return points[inds], values[inds].reshape(-1, self.vdim)
+        return np.zeros((0, 2)), np.zeros((0, self.vdim), dtype=float)
+
+    def _data_interp(self):
+        if self.pending_points:
+            points = list(self.pending_points)
             if self.bounds_are_done:
-                values_interp = self.ip()(self._scale(points_interp))
+                values = self.ip()(self._scale(points))
             else:
                 # Without the bounds the interpolation cannot be done properly,
                 # so we just set everything to zero.
-                values_interp = np.zeros((len(points_interp), self.vdim))
+                values = np.zeros((len(points), self.vdim))
+            return points, values
+        return np.zeros((0, 2)), np.zeros((0, self.vdim), dtype=float)
+
+    def _data_combined(self):
+        points, values = self._data_in_bounds()
+        if not self.pending_points:
+            return points, values
+        points_interp, values_interp = self._data_interp()
+        points_combined = np.vstack([points, points_interp])
+        values_combined = np.vstack([values, values_interp])
+        return points_combined, values_combined
 
-            for point, value in zip(points_interp, values_interp):
-                data_combined[point] = value
-
-        return data_combined
+    def data_combined(self):
+        # Interpolate the unfinished points
+        points, values = self._data_combined()
+        return {tuple(k): v for k, v in zip(points, values)}
 
     def ip(self):
         if self._ip is None:
-            points = self._scale(list(self.data.keys()))
-            values = np.array(list(self.data.values()), dtype=float)
+            points, values = self._data_in_bounds()
+            points = self._scale(points)
             self._ip = interpolate.LinearNDInterpolator(points, values)
         return self._ip
 
     def ip_combined(self):
         if self._ip_combined is None:
-            data_combined = self.data_combined()
-            points = self._scale(list(data_combined.keys()))
-            values = np.array(list(data_combined.values()), dtype=float)
+            points, values = self._data_combined()
+            points = self._scale(points)
             self._ip_combined = interpolate.LinearNDInterpolator(points,
                                                                  values)
         return self._ip_combined
 
+    def inside_bounds(self, xy):
+        x, y = xy
+        (xmin, xmax), (ymin, ymax) = self.bounds
+        return xmin <= x <= xmax and ymin <= y <= ymax
+
     def tell(self, point, value):
         point = tuple(point)
+        self.data[point] = value
+        if not self.inside_bounds(point):
+            return
+        self.pending_points.discard(point)
+        self._ip = None
+        self._stack.pop(point, None)
 
-        if value is None:
-            self._interp.add(point)
-            self._ip_combined = None
-        else:
-            self.data[point] = value
-            self._interp.discard(point)
-            self._ip = None
-
+    def tell_pending(self, point):
+        point = tuple(point)
+        if not self.inside_bounds(point):
+            return
+        self.pending_points.add(point)
+        self._ip_combined = None
         self._stack.pop(point, None)
 
     def _fill_stack(self, stack_till=1):
-        if len(self.data) + len(self._interp) < self.ndim + 1:
+        if len(self.data) + len(self.pending_points) < self.ndim + 1:
             raise ValueError("too few points...")
 
         # Interpolate
@@ -342,13 +371,14 @@ class Learner2D(BaseLearner):
 
         return points_new, losses_new
 
-    def ask(self, n, add_data=True):
-        # Even if add_data is False we add the point such that _fill_stack
+    def ask(self, n, tell_pending=True):
+        # Even if tell_pending is False we add the point such that _fill_stack
         # will return new points, later we remove these points if needed.
         points = list(self._stack.keys())
         loss_improvements = list(self._stack.values())
         n_left = n - len(points)
-        self.tell_many(points[:n], itertools.repeat(None))
+        for p in points[:n]:
+            self.tell_pending(p)
 
         while n_left > 0:
             # The while loop is needed because `stack_till` could be larger
@@ -356,17 +386,18 @@ class Learner2D(BaseLearner):
             # it could fill up till a length smaller than `stack_till`.
             new_points, new_loss_improvements = self._fill_stack(
                 stack_till=max(n_left, self.stack_size))
-            self.tell_many(new_points[:n_left], itertools.repeat(None))
+            for p in points[:n_left]:
+                self.tell_pending(p)
             n_left -= len(new_points)
 
             points += new_points
             loss_improvements += new_loss_improvements
 
-        if not add_data:
+        if not tell_pending:
             self._stack = OrderedDict(zip(points[:self.stack_size],
                                           loss_improvements))
             for point in points[:n]:
-                self._interp.discard(point)
+                self.pending_points.discard(point)
 
         return points[:n], loss_improvements[:n]
 
@@ -379,7 +410,7 @@ class Learner2D(BaseLearner):
         return self._loss
 
     def remove_unfinished(self):
-        self._interp = set()
+        self.pending_points = set()
         for p in self._bounds_points:
             if p not in self.data:
                 self._stack[p] = np.inf

adaptive/learner/learnerND.py

@@ -8,11 +8,12 @@ import numpy as np
 from scipy import interpolate
 import scipy.spatial
 
-from ..notebook_integration import ensure_holoviews
 from .base_learner import BaseLearner
 
-from .triangulation import Triangulation, point_in_simplex, \
-                           circumsphere, simplex_volume_in_embedding
+from ..notebook_integration import ensure_holoviews
+from .triangulation import (Triangulation, point_in_simplex,
+                            circumsphere, simplex_volume_in_embedding)
+from ..utils import restore
 
 
 def volume(simplex, ys=None):
@@ -245,7 +246,7 @@ class LearnerND(BaseLearner):
             return  # we already know about the point
 
         if value is None:
-            return self._tell_pending(point)
+            return self.tell_pending(point)
 
         self._pending.discard(point)
         tri = self.tri
@@ -263,7 +264,7 @@ class LearnerND(BaseLearner):
         simplex = tuple(sorted(simplex))
         return simplex in self.tri.simplices
 
-    def _tell_pending(self, point, simplex=None):
+    def tell_pending(self, point, *, simplex=None):
         point = tuple(point)
         self._pending.add(point)
 
@@ -309,15 +310,23 @@ class LearnerND(BaseLearner):
             heapq.heappush(self._simplex_queue,
                            (-subloss, simplex, subsimplex))
 
-    def ask(self, n=1):
+    def _ask_and_tell_pending(self, n=1):
         xs, losses = zip(*(self._ask() for _ in range(n)))
         return list(xs), list(losses)
 
+    def ask(self, n, tell_pending=True):
+        """Chose points for learners."""
+        if not tell_pending:
+            with restore(self):
+                return self._ask_and_tell_pending(n)
+        else:
+            return self._ask_and_tell_pending(n)
+
     def _ask_bound_point(self):
         # get the next bound point that is still available
         new_point = next(p for p in self._bounds_points
                          if p not in self.data and p not in self._pending)
-        self._tell_pending(new_point)
+        self.tell_pending(new_point)
         return new_point, np.inf
 
     def _ask_point_without_known_simplices(self):
@@ -330,7 +339,7 @@ class LearnerND(BaseLearner):
         p = r * a + b
         p = tuple(p)
 
-        self._tell_pending(p)
+        self.tell_pending(p)
         return p, np.inf
 
     def _pop_highest_existing_simplex(self):
@@ -350,8 +359,8 @@ class LearnerND(BaseLearner):
         # Could not find a simplex, this code should never be reached
         assert self.tri is not None
         raise AssertionError(
-            """Could not find a simplex to. Yet there should always be a simplex 
-            available if LearnerND.tri() is not None"""
+            "Could not find a simplex to. Yet there should always be a simplex "
+            "available if LearnerND.tri() is not None"
         )
 
     def _ask_best_point(self):
@@ -371,7 +380,7 @@ class LearnerND(BaseLearner):
                                                   transform=self._transform))
 
         self._pending_to_simplex[point_new] = simplex
-        self._tell_pending(point_new, simplex)  # O(??)
+        self.tell_pending(point_new, simplex=simplex)  # O(??)
 
         return point_new, loss
 

adaptive/learner/skopt_learner.py

@@ -27,11 +27,13 @@ class SKOptLearner(Optimizer, BaseLearner):
         self.function = function
         super().__init__(**kwargs)
 
-    def tell(self, x, y, fit=True):
-        if y is not None:
-            # 'skopt.Optimizer' takes care of points we
-            # have not got results for.
-            super().tell([x], y, fit)
+    def tell(self, x, y, fit=True):    
+        super().tell([x], y, fit)
+
+    def tell_pending(self, x):
+        # 'skopt.Optimizer' takes care of points we
+        # have not got results for.
+        pass
 
     def remove_unfinished(self):
         pass
@@ -46,7 +48,10 @@ class SKOptLearner(Optimizer, BaseLearner):
             # estimator of loss, but it is the cheapest.
             return 1 - model.score(self.Xi, self.yi)
 
-    def ask(self, n, add_data=True):
+    def ask(self, n, tell_pending=True):
+        if not tell_pending:
+            raise NotImplementedError('Asking points is an irreversible '
+                'action, so use `ask(n, tell_pending=True`.')
         points = super().ask(n)
         # TODO: Choose a better estimate for the loss improvement.
         if self.space.n_dims > 1:

adaptive/tests/test_average_learner.py

+# -*- coding: utf-8 -*-
+
+from ..learner import AverageLearner
+
+
+def test_only_returns_new_points():
+    learner = AverageLearner(lambda x: x, atol=None, rtol=0.01)
+
+    # Only tell it n = 5...10
+    for i in range(5, 10):
+        learner.tell(i, 1)
+
+    learner.tell_pending(0)  # This means it shouldn't return 0 anymore
+
+    assert learner.ask(1)[0][0] == 1
+    assert learner.ask(1)[0][0] == 2
+    assert learner.ask(1)[0][0] == 3
+    assert learner.ask(1)[0][0] == 4
+    assert learner.ask(1)[0][0] == 10

adaptive/tests/test_learner1d.py

+# -*- coding: utf-8 -*-
+
+import random
+import numpy as np
+
+from ..learner import Learner1D
+from ..runner import simple, replay_log
+
+
+def test_pending_loss_intervals():
+    # https://gitlab.kwant-project.org/qt/adaptive/issues/99
+    l = Learner1D(lambda x: x, (0, 4))
+
+    l.tell(0, 0)
+    l.tell(1, 0)
+    l.tell(2, 0)
+    assert set(l.losses_combined.keys()) == {(0, 1), (1, 2)}
+    l.ask(1)
+    assert set(l.losses_combined.keys()) == {(0, 1), (1, 2), (2, 4)}
+    l.tell(3.5, 0)
+    assert set(l.losses_combined.keys()) == {
+        (0, 1), (1, 2), (2, 3.5), (3.5, 4.0)}
+
+
+def test_loss_interpolation_for_unasked_point():
+    # https://gitlab.kwant-project.org/qt/adaptive/issues/99
+    l = Learner1D(lambda x: x, (0, 4))
+
+    l.tell(0, 0)
+    l.tell(1, 0)
+    l.tell(2, 0)
+
+    assert l.ask(1) == ([4], [np.inf])
+    assert l.losses == {(0, 1): 0.25, (1, 2): 0.25}
+    assert l.losses_combined == {(0, 1): 0.25, (1, 2): 0.25, (2, 4.0): np.inf}
+
+    # assert l.ask(1) == ([3], [np.inf])  # XXX: This doesn't return np.inf as loss_improvement...
+    l.ask(1)
+    assert l.losses == {(0, 1): 0.25, (1, 2): 0.25}
+    assert l.losses_combined == {
+        (0, 1): 0.25, (1, 2): 0.25, (2, 3.0): np.inf, (3.0, 4.0): np.inf}
+
+    l.tell(4, 0)
+
+    assert l.losses_combined == {
+        (0, 1): 0.25, (1, 2): 0.25, (2, 3): 0.25, (3, 4): 0.25}
+
+
+def test_first_iteration():
+    """Edge cases where we ask for a few points at the start."""
+    learner = Learner1D(lambda x: None, (-1, 1))
+    points, loss_improvements = learner.ask(2)
+    assert set(points) == set(learner.bounds)
+
+    learner = Learner1D(lambda x: None, (-1, 1))
+    points, loss_improvements = learner.ask(3)
+    assert set(points) == set([-1, 0, 1])
+
+    learner = Learner1D(lambda x: None, (-1, 1))
+    points, loss_improvements = learner.ask(1)
+    assert len(points) == 1 and points[0] in learner.bounds
+    rest = set([-1, 0, 1]) - set(points)
+    points, loss_improvements = learner.ask(2)
+    assert set(points) == set(rest)
+
+    learner = Learner1D(lambda x: None, (-1, 1))
+    points, loss_improvements = learner.ask(1)
+    to_see = set(learner.bounds) - set(points)
+    points, loss_improvements = learner.ask(1)
+    assert set(points) == set(to_see)
+
+    learner = Learner1D(lambda x: None, (-1, 1))
+    learner.tell(1, 0)
+    points, loss_improvements = learner.ask(1)
+    assert points == [-1]
+
+    learner = Learner1D(lambda x: None, (-1, 1))
+    learner.tell(-1, 0)
+    points, loss_improvements = learner.ask(1)
+    assert points == [1]
+
+
+def test_loss_interpolation():
+    learner = Learner1D(lambda _: 0, bounds=(-1, 1))
+
+    learner.tell(-1, 0)
+    learner.tell(1, 0)
+    for i in range(100):
+        # Add a 100 points with either None or 0
+        if random.random() < 0.9:
+            learner.tell(random.uniform(-1, 1), None)
+        else:
+            learner.tell(random.uniform(-1, 1), 0)
+
+    for (x1, x2), loss in learner.losses_combined.items():
+        expected_loss = (x2 - x1) / 2
+        assert abs(expected_loss - loss) < 1e-15, (expected_loss, loss)
+
+
+def _run_on_discontinuity(x_0, bounds):
+
+    def f(x):
+        return -1 if x < x_0 else +1
+
+    learner = Learner1D(f, bounds)
+    while learner.loss() > 0.1:
+        (x,), _ = learner.ask(1)
+        learner.tell(x, learner.function(x))
+
+    return learner
+
+
+def test_termination_on_discontinuities():
+
+    learner = _run_on_discontinuity(0, (-1, 1))
+    smallest_interval = min(abs(a - b) for a, b in learner.losses.keys())
+    assert smallest_interval >= np.finfo(float).eps
+
+    learner = _run_on_discontinuity(1, (-2, 2))
+    smallest_interval = min(abs(a - b) for a, b in learner.losses.keys())
+    assert smallest_interval >= np.finfo(float).eps
+
+    learner = _run_on_discontinuity(0.5E3, (-1E3, 1E3))
+    smallest_interval = min(abs(a - b) for a, b in learner.losses.keys())
+    assert smallest_interval >= 0.5E3 * np.finfo(float).eps
+
+
+def test_order_adding_points():
+    # and https://gitlab.kwant-project.org/qt/adaptive/issues/98
+    l = Learner1D(lambda x: x, (0, 1))
+    l.tell_many([1, 0, 0.5], [0, 0, 0])
+    assert l.losses_combined == {(0, 0.5): 0.5, (0.5, 1): 0.5}
+    assert l.losses == {(0, 0.5): 0.5, (0.5, 1): 0.5}
+    l.ask(1)
+
+
+def test_adding_existing_point_passes_silently():
+    # See https://gitlab.kwant-project.org/qt/adaptive/issues/97
+    l = Learner1D(lambda x: x, (0, 4))
+    l.tell(0, 0)
+    l.tell(1, 0)
+    l.tell(2, 0)
+    l.tell(1, None)
+
+
+def test_loss_at_machine_precision_interval_is_zero():
+    """The loss of an interval smaller than _dx_eps
+    should be set to zero."""
+    def f(x):
+        return 1 if x == 0 else 0
+
+    def goal(l):
+        return l.loss() < 0.01 or l.npoints >= 1000
+
+    learner = Learner1D(f, bounds=(-1, 1))
+    simple(learner, goal=goal)
+
+    # this means loss < 0.01 was reached
+    assert learner.npoints != 1000
+
+
+def small_deviations(x):
+    return 0 if x <= 1 else 1 + 10**(-random.randint(12, 14))
+
+
+def test_small_deviations():
+    """This tests whether the Learner1D can handle small deviations.
+    See https://gitlab.kwant-project.org/qt/adaptive/merge_requests/73 and
+    https://gitlab.kwant-project.org/qt/adaptive/issues/61."""
+
+    eps = 5e-14
+    learner = Learner1D(small_deviations, bounds=(1 - eps, 1 + eps))
+
+    # Some non-determinism is needed to make this test fail so we keep
+    # a list of points that will be evaluated later to emulate
+    # parallel execution
+    stash = []
+
+    for i in range(100):
+        xs, _ = learner.ask(10)
+
+        # Save 5 random points out of `xs` for later
+        random.shuffle(xs)
+        for _ in range(5):
+            stash.append(xs.pop())
+
+        for x in xs:
+            learner.tell(x, learner.function(x))
+
+        # Evaluate and add 5 random points from `stash`
+        random.shuffle(stash)
+        for _ in range(5):
+            learner.tell(stash.pop(), learner.function(x))
+
+        if learner.loss() == 0:
+            # If this condition is met, the learner can't return any
+            # more points.
+            break
+
+
+def test_uniform_sampling1D_v2():
+    def check(known, expect):
+        def f(x): return x
+        learner = Learner1D(f, bounds=(-1, 1))
+        for x in known:
+            learner.tell(x, f(x))
+        pts, _ = learner.ask(len(expect))
+        assert set(pts) == expect
+
+    check([-1, 0, 1], {-0.5, 0.5})
+    check([-1, -0.5, 1], {0, 0.5})
+    check([-1, 0.5, 1], {-0.5, 0})
+    check([-1, 0], {1})
+    # although the following test might be unexpected, this is indeed correct
+    # given the default loss function
+    check([-1, 0], {-.5, 1})
+    check([-1, -.5], {-.75, 1})
+    check([-1, -.5], {-.75, .25, 1})
+
+
+def test_add_data_unordered():
+    # see https://gitlab.kwant-project.org/qt/adaptive/issues/95
+    learner = Learner1D(lambda x: x, bounds=(-1, 1))
+    xs = [-1, 1, 0]
+
+    ys = [learner.function(x) for x in xs]
+    for x, y in zip(xs, ys):
+        learner.tell(x, y)
+
+    learner.ask(3)
+
+
+def test_ask_does_not_return_known_points_when_returning_bounds():
+    learner = Learner1D(lambda x: None, (-1, 1))
+    learner.tell(0, 0)
+    points, _ = learner.ask(3)
+    assert 0 not in points

adaptive/tests/test_learnernd.py

+# -*- coding: utf-8 -*-
+
+from ..learner import LearnerND
+from ..runner import replay_log
+
+
+def test_faiure_case_LearnerND():
+    log = [
+        ('ask', 4),
+        ('tell', (-1, -1, -1), 1.607873907219222e-101),
+        ('tell', (-1, -1, 1), 1.607873907219222e-101),
+        ('ask', 2),
+        ('tell', (-1, 1, -1), 1.607873907219222e-101),
+        ('tell', (-1, 1, 1), 1.607873907219222e-101),
+        ('ask', 2),
+        ('tell', (1, -1, 1), 2.0),
+        ('tell', (1, -1, -1), 2.0),
+        ('ask', 2),
+        ('tell', (0.0, 0.0, 0.0), 4.288304431237686e-06),
+        ('tell', (1, 1, -1), 2.0)
+    ]
+    learner = LearnerND(lambda *x: x, bounds=[(-1, 1), (-1, 1), (-1, 1)])
+    replay_log(learner, log)

adaptive/tests/test_learner.py → adaptive/tests/test_learners.py

@@ -14,12 +14,6 @@ import pytest
 from ..learner import *
 from ..runner import simple, replay_log
 
-try:
-    import skopt
-    with_scikit_optimize = True
-except ModuleNotFoundError:
-    with_scikit_optimize = False
-
 
 def generate_random_parametrization(f):
     """Return a realization of 'f' with parameters bound to random values.
@@ -108,7 +102,8 @@ def run_with(*learner_types):
             # Check if learner was marked with our `xfail` decorator
             # XXX: doesn't work when feeding kwargs to xfail.
             if is_xfail:
-                pars.append(pytest.param(l, f, dict(k), marks=[pytest.mark.xfail]))
+                pars.append(pytest.param(l, f, dict(k),
+                                         marks=[pytest.mark.xfail]))
             else:
                 pars.append((l, f, dict(k)))
     return pytest.mark.parametrize('learner_type, f, learner_kwargs', pars)
@@ -128,23 +123,7 @@ def ask_randomly(learner, rounds, points):
     return xs, ls
 
 
-@pytest.mark.skipif(not with_scikit_optimize,
-                    reason='scikit-optimize is not installed')
-def test_skopt_learner_runs():
-    """The SKOptLearner provides very few guarantees about its
-       behaviour, so we only test the most basic usage
-    """
-
-    def g(x, noise_level=0.1):
-        return (np.sin(5 * x) * (1 - np.tanh(x ** 2))
-                + np.random.randn() * noise_level)
-
-    learner = SKOptLearner(g, dimensions=[(-2., 2.)])
-
-    for _ in range(11):
-        (x,), _ = learner.ask(1)
-        learner.tell(x, learner.function(x))
-
+# Tests
 
 @run_with(Learner1D)
 def test_uniform_sampling1D(learner_type, f, learner_kwargs):
@@ -188,6 +167,19 @@ def test_uniform_sampling2D(learner_type, f, learner_kwargs):
     assert max(distances) < math.sqrt(dx**2 + dy**2)
 
 
+@pytest.mark.parametrize('learner_type, bounds', [
+    (Learner1D, (-1, 1)),
+    (Learner2D, [(-1, 1), (-1, 1)]),
+    (LearnerND, [(-1, 1), (-1, 1), (-1, 1)]),
+])
+def test_learner_accepts_lists(learner_type, bounds):
+    def f(x):
+        return [0, 1]
+
+    learner = learner_type(f, bounds=bounds)
+    simple(learner, goal=lambda l: l.npoints > 10)
+
+
 @run_with(xfail(Learner1D), Learner2D, LearnerND)
 def test_adding_existing_data_is_idempotent(learner_type, f, learner_kwargs):
     """Adding already existing data is an idempotent operation.
@@ -359,122 +351,38 @@ def test_learner_performance_is_invariant_under_scaling(learner_type, f, learner
     assert abs(learner.loss() - control.loss()) / learner.loss() < 1e-11
 
 
-def test_learner1d_first_iteration():
-    """Edge cases where we ask for a few points at the start."""
-    learner = Learner1D(lambda x: None, (-1, 1))
-    points, loss_improvements = learner.ask(2)
-    assert set(points) == set(learner.bounds)
-
-    learner = Learner1D(lambda x: None, (-1, 1))
-    points, loss_improvements = learner.ask(3)
-    assert set(points) == set([-1, 0, 1])
-
-    learner = Learner1D(lambda x: None, (-1, 1))
-    points, loss_improvements = learner.ask(1)
-    assert len(points) == 1 and points[0] in learner.bounds
-    rest = set([-1, 0, 1]) - set(points)
-    points, loss_improvements = learner.ask(2)
-    assert set(points) == set(rest)
-
-    learner = Learner1D(lambda x: None, (-1, 1))
-    points, loss_improvements = learner.ask(1)
-    to_see = set(learner.bounds) - set(points)
-    points, loss_improvements = learner.ask(1)
-    assert set(points) == set(to_see)
-
-    learner = Learner1D(lambda x: None, (-1, 1))
-    learner.tell(1, 0)
-    points, loss_improvements = learner.ask(1)
-    assert points == [-1]
-
-    learner = Learner1D(lambda x: None, (-1, 1))
-    learner.tell(-1, 0)
-    points, loss_improvements = learner.ask(1)
-    assert points == [1]
-
-
-def _run_on_discontinuity(x_0, bounds):
-
-    def f(x):
-        return -1 if x < x_0 else +1
-
-    learner = Learner1D(f, bounds)
-    while learner.loss() > 0.1:
-        (x,), _ = learner.ask(1)
-        learner.tell(x, learner.function(x))
+# XXX: The LearnerND shouldn't fail, see https://gitlab.kwant-project.org/qt/adaptive/issues/105
+@run_with(Learner1D, Learner2D, xfail(LearnerND), AverageLearner)
+def test_balancing_learner(learner_type, f, learner_kwargs):
+    """Test if the BalancingLearner works with the different types of learners."""
+    learners = [learner_type(generate_random_parametrization(f), **learner_kwargs)
+                for i in range(4)]
 
-    return learner
+    learner = BalancingLearner(learners)
 
-
-def test_termination_on_discontinuities():
-
-    learner = _run_on_discontinuity(0, (-1, 1))
-    smallest_interval = min(abs(a - b) for a, b in learner.losses.keys())
-    assert smallest_interval >= np.finfo(float).eps
-
-    learner = _run_on_discontinuity(1, (-2, 2))
-    smallest_interval = min(abs(a - b) for a, b in learner.losses.keys())
-    assert smallest_interval >= np.finfo(float).eps
-
-    learner = _run_on_discontinuity(0.5E3, (-1E3, 1E3))
-    smallest_interval = min(abs(a - b) for a, b in learner.losses.keys())
-    assert smallest_interval >= 0.5E3 * np.finfo(float).eps
-
-
-def test_loss_at_machine_precision_interval_is_zero():
-    """The loss of an interval smaller than _dx_eps
-    should be set to zero."""
-    def f(x):
-        return 1 if x == 0 else 0
-
-    def goal(l):
-        return l.loss() < 0.01 or l.npoints >= 1000
-
-    learner = Learner1D(f, bounds=(-1, 1))
-    simple(learner, goal=goal)
-
-    # this means loss < 0.01 was reached
-    assert learner.npoints != 1000
-
-
-def small_deviations(x):
-    import random
-    return 0 if x <= 1 else 1 + 10**(-random.randint(12, 14))
-
-
-def test_small_deviations():
-    """This tests whether the Learner1D can handle small deviations.
-    See https://gitlab.kwant-project.org/qt/adaptive/merge_requests/73 and
-    https://gitlab.kwant-project.org/qt/adaptive/issues/61."""
-
-    eps = 5e-14
-    learner = Learner1D(small_deviations, bounds=(1 - eps, 1 + eps))
-
-    # Some non-determinism is needed to make this test fail so we keep
-    # a list of points that will be evaluated later to emulate
-    # parallel execution
+    # Emulate parallel execution
     stash = []
 
     for i in range(100):
-        xs, _ = learner.ask(10)
+        n = random.randint(1, 10)
+        m = random.randint(0, n)
+        xs, _ = learner.ask(n, tell_pending=False)
 
-        # Save 5 random points out of `xs` for later
+        # Save 'm' random points out of `xs` for later
         random.shuffle(xs)
-        for _ in range(5):
+        for _ in range(m):
             stash.append(xs.pop())
 
         for x in xs:
             learner.tell(x, learner.function(x))
 
-        # Evaluate and add 5 random points from `stash`
+        # Evaluate and add 'm' random points from `stash`
         random.shuffle(stash)
-        for _ in range(5):
-            learner.tell(stash.pop(), learner.function(x))
+        for _ in range(m):
+            x = stash.pop()
+            learner.tell(x, learner.function(x))
 
-        if learner.loss() == 0:
-            # If this condition is met, the learner can't return any
-            # more points.
-            break
+    assert all(l.npoints > 10 for l in learner.learners), [l.npoints for l in learner.learners]
 
 
 @pytest.mark.xfail
@@ -496,22 +404,3 @@ def test_learner_subdomain(learner_type, f, learner_kwargs):
        perform 'similarly' to learners defined on that subdomain only."""
     # XXX: not sure how to implement this. How do we measure "performance"?
     raise NotImplementedError()
-
-
-def test_faiure_case_LearnerND():
-    log = [
-        ('ask', 4),
-        ('tell', (-1, -1, -1), 1.607873907219222e-101),
-        ('tell', (-1, -1, 1), 1.607873907219222e-101),
-        ('ask', 2),
-        ('tell', (-1, 1, -1), 1.607873907219222e-101),
-        ('tell', (-1, 1, 1), 1.607873907219222e-101),
-        ('ask', 2),
-        ('tell', (1, -1, 1), 2.0),
-        ('tell', (1, -1, -1), 2.0),
-        ('ask', 2),
-        ('tell', (0.0, 0.0, 0.0), 4.288304431237686e-06),
-        ('tell', (1, 1, -1), 2.0)
-    ]
-    learner = LearnerND(lambda *x: x, bounds=[(-1, 1), (-1, 1), (-1, 1)])
-    replay_log(learner, log)

adaptive/tests/test_skopt_learner.py

+# -*- coding: utf-8 -*-
+
+import random
+import numpy as np
+
+import pytest
+
+
+try:
+    import skopt
+    with_scikit_optimize = True
+    from ..learner import SKOptLearner
+except ModuleNotFoundError:
+    with_scikit_optimize = False
+
+
+@pytest.mark.skipif(not with_scikit_optimize,
+                    reason='scikit-optimize is not installed')
+def test_skopt_learner_runs():
+    """The SKOptLearner provides very few guarantees about its
+       behaviour, so we only test the most basic usage
+    """
+
+    def g(x, noise_level=0.1):
+        return (np.sin(5 * x) * (1 - np.tanh(x ** 2))
+                + np.random.randn() * noise_level)
+
+    learner = SKOptLearner(g, dimensions=[(-2., 2.)])
+
+    for _ in range(11):
+        (x,), _ = learner.ask(1)
+        learner.tell(x, learner.function(x))

learner.ipynb

 %% Cell type:markdown id: tags:
 
 # Adaptive
 
 %% Cell type:markdown id: tags:
 
 [`adaptive`](https://gitlab.kwant-project.org/qt/adaptive-evaluation) is a package for adaptively sampling functions with support for parallel evaluation.
 
 This is an introductory notebook that shows some basic use cases.
 
 `adaptive` needs at least Python 3.6, and the following packages:
 
 + `scipy`
 + `sortedcontainers`
 
 Additionally `adaptive` has lots of extra functionality that makes it simple to use from Jupyter notebooks.
 This extra functionality depends on the following packages
 
 + `ipykernel>=4.8.0`
 + `jupyter_client>=5.2.2`
 + `holoviews`
 + `bokeh`
 + `ipywidgets`
 
 %% Cell type:code id: tags:
 
 ``` 
 import adaptive
 adaptive.notebook_extension()
 
 # Import modules that are used in multiple cells
 import holoviews as hv
 import numpy as np
 from functools import partial
 import random
 ```
 
 %% Cell type:markdown id: tags:
 
 # 1D function learner
 
 %% Cell type:markdown id: tags:
 
 We start with the most common use-case: sampling a 1D function $\ f: ℝ → ℝ$.
 
 We will use the following function, which is a smooth (linear) background with a sharp peak at a random location:
 
 %% Cell type:code id: tags:
 
 ``` 
 offset = random.uniform(-0.5, 0.5)
 
 def f(x, offset=offset, wait=True):
     from time import sleep
     from random import random
 
     a = 0.01
     if wait:
         sleep(random())
     return x + a**2 / (a**2 + (x - offset)**2)
 ```
 
 %% Cell type:markdown id: tags:
 
 We start by initializing a 1D "learner", which will suggest points to evaluate, and adapt its suggestions as more and more points are evaluated.
 
 %% Cell type:code id: tags:
 
 ``` 
 learner = adaptive.Learner1D(f, bounds=(-1, 1))
 ```
 
 %% Cell type:markdown id: tags:
 
 Next we create a "runner" that will request points from the learner and evaluate 'f' on them.
 
 By default on Unix-like systems the runner will evaluate the points in parallel using local processes ([`concurrent.futures.ProcessPoolExecutor`](https://docs.python.org/3/library/concurrent.futures.html#processpoolexecutor)).
 
 On Windows systems the runner will try to use a [`distributed.Client`](https://distributed.readthedocs.io/en/latest/client.html) if [`distributed`](https://distributed.readthedocs.io/en/latest/index.html) is installed. A `ProcessPoolExecutor` cannot be used on Windows for reasons.
 
 %% Cell type:code id: tags:
 
 ``` 
 # The end condition is when the "loss" is less than 0.1. In the context of the
 # 1D learner this means that we will resolve features in 'func' with width 0.1 or wider.
 runner = adaptive.Runner(learner, goal=lambda l: l.loss() < 0.05)
 runner.live_info()
 ```
 
 %% Cell type:markdown id: tags:
 
 When instantiated in a Jupyter notebook the runner does its job in the background and does not block the IPython kernel.
 We can use this to create a plot that updates as new data arrives:
 
 %% Cell type:code id: tags:
 
 ``` 
 runner.live_plot(update_interval=0.1)
 ```
 
 %% Cell type:markdown id: tags:
 
 We can now compare the adaptive sampling to a homogeneous sampling with the same number of points:
 
 %% Cell type:code id: tags:
 
 ``` 
 if not runner.task.done():
     raise RuntimeError('Wait for the runner to finish before executing the cells below!')
 ```
 
 %% Cell type:code id: tags:
 
 ``` 
 learner2 = adaptive.Learner1D(f, bounds=learner.bounds)
 
 xs = np.linspace(*learner.bounds, len(learner.data))
 learner2.tell_many(xs, map(partial(f, wait=False), xs))
 
 learner.plot() + learner2.plot()
 ```
 
 %% Cell type:markdown id: tags:
 
 # 2D function learner
 
 %% Cell type:markdown id: tags:
 
 Besides 1D functions, we can also learn 2D functions: $\ f: ℝ^2 → ℝ$
 
 %% Cell type:code id: tags:
 
 ``` 
 def ring(xy, wait=True):
     import numpy as np
     from time import sleep
     from random import random
     if wait:
         sleep(random()/10)
     x, y = xy
     a = 0.2
     return x + np.exp(-(x**2 + y**2 - 0.75**2)**2/a**4)
 
 learner = adaptive.Learner2D(ring, bounds=[(-1, 1), (-1, 1)])
 ```
 
 %% Cell type:code id: tags:
 
 ``` 
 runner = adaptive.Runner(learner, goal=lambda l: l.loss() < 0.01)
 runner.live_info()
 ```
 
 %% Cell type:code id: tags:
 
 ``` 
 def plot(learner):
     plot = learner.plot(tri_alpha=0.2)
     title = f'loss={learner._loss:.3f}, n_points={learner.npoints}'
     return (plot.Image
             + plot.EdgePaths.I.opts(plot=dict(title_format=title))
             + plot)
 
 runner.live_plot(plotter=plot, update_interval=0.1)
 ```
 
 %% Cell type:code id: tags:
 
 ``` 
 %%opts EdgePaths (color='w')
 
 import itertools
 
 # Create a learner and add data on homogeneous grid, so that we can plot it
 learner2 = adaptive.Learner2D(ring, bounds=learner.bounds)
 n = int(learner.npoints**0.5)
 xs, ys = [np.linspace(*bounds, n) for bounds in learner.bounds]
 xys = list(itertools.product(xs, ys))
 learner2.tell_many(xys, map(partial(ring, wait=False), xys))
 
 (learner2.plot(n).relabel('Homogeneous grid') + learner.plot().relabel('With adaptive') +
  learner2.plot(n, tri_alpha=0.4) + learner.plot(tri_alpha=0.4)).cols(2)
 ```
 
 %% Cell type:markdown id: tags:
 
-# N-dimensional function learner
-Besides 1 and 2 dimensional functions, we can also learn N-D functions: $\ f: ℝ^N → ℝ, N \ge 2$
-
-Do keep in mind the speed and [effectiveness](https://en.wikipedia.org/wiki/Curse_of_dimensionality) of the learner drops quickly with increasing number of dimensions.
-
-%% Cell type:code id: tags:
-
-``` 
-# this step takes a lot of time, it will finish at about 3300 points, which can take up to 6 minutes
-def sphere(xyz):
-    x, y, z = xyz
-    a = 0.4
-    return x + z**2 + np.exp(-(x**2 + y**2 + z**2 - 0.75**2)**2/a**4)
-
-learner = adaptive.LearnerND(sphere, bounds=[(-1, 1), (-1, 1), (-1, 1)])
-runner = adaptive.Runner(learner, goal=lambda l: l.loss() < 0.01)
-runner.live_info()
-```
-
-%% Cell type:markdown id: tags:
-
-Let's plot 2D slices of the 3D function
-
-%% Cell type:code id: tags:
-
-``` 
-def plot_cut(x, direction, learner=learner):
-    cut_mapping = {'xyz'.index(direction): x}
-    return learner.plot_slice(cut_mapping, n=100)
-
-dm = hv.DynamicMap(plot_cut, kdims=['value', 'direction'])
-dm.redim.values(value=np.linspace(-1, 1), direction=list('xyz'))
-```
-
-%% Cell type:markdown id: tags:
-
-Or we can plot 1D slices
-
-%% Cell type:code id: tags:
-
-``` 
-%%opts Path {+framewise}
-def plot_cut(x1, x2, directions, learner=learner):
-    cut_mapping = {'xyz'.index(d): x for d, x in zip(directions, [x1, x2])}
-    return learner.plot_slice(cut_mapping)
-
-dm = hv.DynamicMap(plot_cut, kdims=['v1', 'v2', 'directions'])
-dm.redim.values(v1=np.linspace(-1, 1),
-                v2=np.linspace(-1, 1),
-                directions=['xy', 'xz', 'yz'])
-```
-
-%% Cell type:markdown id: tags:
-
-The plots show some wobbles while the original function was smooth, this is a result of the fact that the learner chooses points in 3 dimensions and the simplices are not in the same face as we try to interpolate our lines. However, as always, when you sample more points the graph will become gradually smoother.
-
-%% Cell type:markdown id: tags:
-
 # Averaging learner
 
 %% Cell type:markdown id: tags:
 
 The next type of learner averages a function until the uncertainty in the average meets some condition.
 
 This is useful for sampling a random variable. The function passed to the learner must formally take a single parameter,
 which should be used like a "seed" for the (pseudo-) random variable (although in the current implementation the seed parameter can be ignored by the function).
 
 %% Cell type:code id: tags:
 
 ``` 
 def g(n):
     import random
     from time import sleep
     sleep(random.random() / 1000)
     # Properly save and restore the RNG state
     state = random.getstate()
     random.seed(n)
     val = random.gauss(0.5, 1)
     random.setstate(state)
     return val
 ```
 
 %% Cell type:code id: tags:
 
 ``` 
 learner = adaptive.AverageLearner(g, atol=None, rtol=0.01)
 runner = adaptive.Runner(learner, goal=lambda l: l.loss() < 2)
 runner.live_info()
 ```
 
 %% Cell type:code id: tags:
 
 ``` 
 runner.live_plot(update_interval=0.1)
 ```
 
 %% Cell type:markdown id: tags:
 
 # 1D integration learner with `cquad`
 
 %% Cell type:markdown id: tags:
 
 This learner learns a 1D function and calculates the integral and error of the integral with it. It is based on Pedro Gonnet's [implementation](https://www.academia.edu/1976055/Adaptive_quadrature_re-revisited).
 
 Let's try the following function with cusps (that is difficult to integrate):
 
 %% Cell type:code id: tags:
 
 ``` 
 def f24(x):
     return np.floor(np.exp(x))
 
 xs = np.linspace(0, 3, 200)
 hv.Scatter((xs, [f24(x) for x in xs]))
 ```
 
 %% Cell type:markdown id: tags:
 
 Just to prove that this really is a difficult to integrate function, let's try a familiar function integrator `scipy.integrate.quad`, which will give us warnings that it encounters difficulties.
 
 %% Cell type:code id: tags:
 
 ``` 
 import scipy.integrate
 scipy.integrate.quad(f24, 0, 3)
 ```
 
 %% Cell type:markdown id: tags:
 
 We initialize a learner again and pass the bounds and relative tolerance we want to reach. Then in the `Runner` we pass `goal=lambda l: l.done()` where `learner.done()` is `True` when the relative tolerance has been reached.
 
 %% Cell type:code id: tags:
 
 ``` 
 from adaptive.runner import SequentialExecutor
 
 learner = adaptive.IntegratorLearner(f24, bounds=(0, 3), tol=1e-8)
 
 # We use a SequentialExecutor, which runs the function to be learned in *this* process only. This means we don't pay
 # the overhead of evaluating the function in another process.
 runner = adaptive.Runner(learner, executor=SequentialExecutor(), goal=lambda l: l.done())
 runner.live_info()
 ```
 
 %% Cell type:markdown id: tags:
 
 Now we could do the live plotting again, but lets just wait untill the runner is done.
 
 %% Cell type:code id: tags:
 
 ``` 
 if not runner.task.done():
     raise RuntimeError('Wait for the runner to finish before executing the cells below!')
 ```
 
 %% Cell type:code id: tags:
 
 ``` 
 print('The integral value is {} with the corresponding error of {}'.format(learner.igral, learner.err))
 learner.plot()
 ```
 
 %% Cell type:markdown id: tags:
 
 # 1D learner with vector output: `f:ℝ → ℝ^N`
 
 %% Cell type:markdown id: tags:
 
 Sometimes you may want to learn a function with vector output:
 
 %% Cell type:code id: tags:
 
 ``` 
 random.seed(0)
 offsets = [random.uniform(-0.8, 0.8) for _ in range(3)]
 
 # sharp peaks at random locations in the domain
 def f_levels(x, offsets=offsets):
     a = 0.01
     return np.array([offset + x + a**2 / (a**2 + (x - offset)**2)
                      for offset in offsets])
 ```
 
 %% Cell type:markdown id: tags:
 
 `adaptive` has you covered! The `Learner1D` can be used for such functions:
 
 %% Cell type:code id: tags:
 
 ``` 
 learner = adaptive.Learner1D(f_levels, bounds=(-1, 1))
 runner = adaptive.Runner(learner, goal=lambda l: l.loss() < 0.01)
 runner.live_info()
 ```
 
 %% Cell type:code id: tags:
 
 ``` 
 runner.live_plot(update_interval=0.1)
 ```
 
 %% Cell type:markdown id: tags:
 
+# N-dimensional function learner (beta)
+Besides 1 and 2 dimensional functions, we can also learn N-D functions: $\ f: ℝ^N → ℝ, N \ge 2$
+
+Do keep in mind the speed and [effectiveness](https://en.wikipedia.org/wiki/Curse_of_dimensionality) of the learner drops quickly with increasing number of dimensions.
+
+%% Cell type:code id: tags:
+
+``` 
+# this step takes a lot of time, it will finish at about 3300 points, which can take up to 6 minutes
+def sphere(xyz):
+    x, y, z = xyz
+    a = 0.4
+    return x + z**2 + np.exp(-(x**2 + y**2 + z**2 - 0.75**2)**2/a**4)
+
+learner = adaptive.LearnerND(sphere, bounds=[(-1, 1), (-1, 1), (-1, 1)])
+runner = adaptive.Runner(learner, goal=lambda l: l.loss() < 0.01)
+runner.live_info()
+```
+
+%% Cell type:markdown id: tags:
+
+Let's plot 2D slices of the 3D function
+
+%% Cell type:code id: tags:
+
+``` 
+def plot_cut(x, direction, learner=learner):
+    cut_mapping = {'xyz'.index(direction): x}
+    return learner.plot_slice(cut_mapping, n=100)
+
+dm = hv.DynamicMap(plot_cut, kdims=['value', 'direction'])
+dm.redim.values(value=np.linspace(-1, 1), direction=list('xyz'))
+```
+
+%% Cell type:markdown id: tags:
+
+Or we can plot 1D slices
+
+%% Cell type:code id: tags:
+
+``` 
+%%opts Path {+framewise}
+def plot_cut(x1, x2, directions, learner=learner):
+    cut_mapping = {'xyz'.index(d): x for d, x in zip(directions, [x1, x2])}
+    return learner.plot_slice(cut_mapping)
+
+dm = hv.DynamicMap(plot_cut, kdims=['v1', 'v2', 'directions'])
+dm.redim.values(v1=np.linspace(-1, 1),
+                v2=np.linspace(-1, 1),
+                directions=['xy', 'xz', 'yz'])
+```
+
+%% Cell type:markdown id: tags:
+
+The plots show some wobbles while the original function was smooth, this is a result of the fact that the learner chooses points in 3 dimensions and the simplices are not in the same face as we try to interpolate our lines. However, as always, when you sample more points the graph will become gradually smoother.
+
+%% Cell type:markdown id: tags:
+
 # Custom adaptive logic for 1D and 2D
 
 %% Cell type:markdown id: tags:
 
 `Learner1D` and `Learner2D` both work on the principle of subdividing their domain into subdomains, and assigning a property to each subdomain, which we call the *loss*. The algorithm for choosing the best place to evaluate our function is then simply *take the subdomain with the largest loss and add a point in the center, creating new subdomains around this point*.
 
 The *loss function* that defines the loss per subdomain is the canonical place to define what regions of the domain are "interesting".
 The default loss function for `Learner1D` and `Learner2D` is sufficient for a wide range of common cases, but it is by no means a panacea. For example, the default loss function will tend to get stuck on divergences.
 
 Both the `Learner1D` and `Learner2D` allow you to specify a *custom loss function*. Below we illustrate how you would go about writing your own loss function. The documentation for `Learner1D` and `Learner2D` specifies the signature that your loss function needs to have in order for it to work with `adaptive`.
 
 
 Say we want to properly sample a function that contains divergences. A simple (but naive) strategy is to *uniformly* sample the domain:
 
 %% Cell type:code id: tags:
 
 ``` 
 def uniform_sampling_1d(interval, scale, function_values):
     # Note that we never use 'function_values'; the loss is just the size of the subdomain
     x_left, x_right = interval
     x_scale, _ = scale
     dx = (x_right - x_left) / x_scale
     return dx
 
 def f_divergent_1d(x):
     return 1 / x**2
 
 learner = adaptive.Learner1D(f_divergent_1d, (-1, 1), loss_per_interval=uniform_sampling_1d)
 runner = adaptive.BlockingRunner(learner, goal=lambda l: l.loss() < 0.01)
 learner.plot().select(y=(0, 10000))
 ```
 
 %% Cell type:code id: tags:
 
 ``` 
 %%opts EdgePaths (color='w') Image [logz=True]
 
 from adaptive.runner import SequentialExecutor
 
 def uniform_sampling_2d(ip):
     from adaptive.learner.learner2D import areas
     A = areas(ip)
     return np.sqrt(A)
 
 def f_divergent_2d(xy):
     x, y = xy
     return 1 / (x**2 + y**2)
 
 learner = adaptive.Learner2D(f_divergent_2d, [(-1, 1), (-1, 1)], loss_per_triangle=uniform_sampling_2d)
 
 # this takes a while, so use the async Runner so we know *something* is happening
 runner = adaptive.Runner(learner, goal=lambda l: l.loss() < 0.02)
 runner.live_info()
 runner.live_plot(update_interval=0.2,
                  plotter=lambda l: l.plot(tri_alpha=0.3).relabel('1 / (x^2 + y^2) in log scale'))
 ```
 
 %% Cell type:markdown id: tags:
 
 The uniform sampling strategy is a common case to benchmark against, so the 1D and 2D versions are included in `adaptive` as `adaptive.learner.learner1D.uniform_sampling` and `adaptive.learner.learner2D.uniform_sampling`.
 
 %% Cell type:markdown id: tags:
 
 ### Doing better
 
 Of course, using `adaptive` for uniform sampling is a bit of a waste!
 
 Let's see if we can do a bit better. Below we define a loss per subdomain that scales with the degree of nonlinearity of the function (this is very similar to the default loss function for `Learner2D`), but which is 0 for subdomains smaller than a certain area, and infinite for subdomains larger than a certain area.
 
 A loss defined in this way means that the adaptive algorithm will first prioritise subdomains that are too large (infinite loss). After all subdomains are appropriately small it will prioritise places where the function is very nonlinear, but will ignore subdomains that are too small (0 loss).
 
 %% Cell type:code id: tags:
 
 ``` 
 %%opts EdgePaths (color='w') Image [logz=True]
 
 def resolution_loss(ip, min_distance=0, max_distance=1):
     """min_distance and max_distance should be in between 0 and 1
     because the total area is normalized to 1."""
 
     from adaptive.learner.learner2D import areas, deviations
 
     A = areas(ip)
 
     # 'deviations' returns an array of shape '(n, len(ip))', where
     # 'n' is the  is the dimension of the output of the learned function
     # In this case we know that the learned function returns a scalar,
     # so 'deviations' returns an array of shape '(1, len(ip))'.
     # It represents the deviation of the function value from a linear estimate
     # over each triangular subdomain.
     dev = deviations(ip)[0]
 
     # we add terms of the same dimension: dev == [distance], A == [distance**2]
     loss = np.sqrt(A) * dev + A
 
     # Setting areas with a small area to zero such that they won't be chosen again
     loss[A < min_distance**2] = 0
 
     # Setting triangles that have a size larger than max_distance to infinite loss
     loss[A > max_distance**2] = np.inf
 
     return loss
 
 loss = partial(resolution_loss, min_distance=0.01)
 
 learner = adaptive.Learner2D(f_divergent_2d, [(-1, 1), (-1, 1)], loss_per_triangle=loss)
 runner = adaptive.BlockingRunner(learner, goal=lambda l: l.loss() < 0.02)
 learner.plot(tri_alpha=0.3).relabel('1 / (x^2 + y^2) in log scale')
 ```
 
 %% Cell type:markdown id: tags:
 
 Awesome! We zoom in on the singularity, but not at the expense of sampling the rest of the domain a reasonable amount.
 
 The above strategy is available as `adaptive.learner.learner2D.resolution_loss`.
 
 %% Cell type:markdown id: tags:
 
 # Balancing learner
 
 %% Cell type:markdown id: tags:
 
 The balancing learner is a "meta-learner" that takes a list of learners. When you request a point from the balancing learner, it will query all of its "children" to figure out which one will give the most improvement.
 
 The balancing learner can for example be used to implement a poor-man's 2D learner by using the `Learner1D`.
 
 %% Cell type:code id: tags:
 
 ``` 
 def h(x, offset=0):
     a = 0.01
     return x + a**2 / (a**2 + (x - offset)**2)
 
 learners = [adaptive.Learner1D(partial(h, offset=random.uniform(-1, 1)),
             bounds=(-1, 1)) for i in range(10)]
 
 bal_learner = adaptive.BalancingLearner(learners)
 runner = adaptive.Runner(bal_learner, goal=lambda l: l.loss() < 0.01)
 runner.live_info()
 ```
 
 %% Cell type:code id: tags:
 
 ``` 
 plotter = lambda learner: hv.Overlay([L.plot() for L in learner.learners])
 runner.live_plot(plotter=plotter, update_interval=0.1)
 ```
 
 %% Cell type:markdown id: tags:
 
 Often one wants to create a set of `learner`s for a cartesian product of parameters. For that particular case we've added a `classmethod` called `from_product`. See how it works below
 
 %% Cell type:code id: tags:
 
 ``` 
 from scipy.special import eval_jacobi
 
 def jacobi(x, n, alpha, beta): return eval_jacobi(n, alpha, beta, x)
 
 combos = {
     'n': [1, 2, 4, 8],
     'alpha': np.linspace(0, 2, 3),
     'beta': np.linspace(0, 1, 5),
 }
 
 learner = adaptive.BalancingLearner.from_product(
     jacobi, adaptive.Learner1D, dict(bounds=(0, 1)), combos)
 
 runner = adaptive.BlockingRunner(learner, goal=lambda l: l.loss() < 0.01)
 
 # The `cdims` will automatically be set when using `from_product`, so
 # `plot()` will return a HoloMap with correctly labeled sliders.
 learner.plot().overlay('beta').grid()
 ```
 
 %% Cell type:markdown id: tags:
 
 # DataSaver
 
 %% Cell type:markdown id: tags:
 
 If the function that you want to learn returns a value along with some metadata, you can wrap your learner in an `adaptive.DataSaver`.
 
 In the following example the function to be learned returns its result and the execution time in a dictionary:
 
 %% Cell type:code id: tags:
 
 ``` 
 from operator import itemgetter
 
 def f_dict(x):
     """The function evaluation takes roughly the time we `sleep`."""
     import random
     from time import sleep
 
     waiting_time = random.random()
     sleep(waiting_time)
     a = 0.01
     y = x + a**2 / (a**2 + x**2)
     return {'y': y, 'waiting_time': waiting_time}
 
 # Create the learner with the function that returns a 'dict'
 # This learner cannot be run directly, as Learner1D does not know what to do with the 'dict'
 _learner = adaptive.Learner1D(f_dict, bounds=(-1, 1))
 
 # Wrapping the learner with 'adaptive.DataSaver' and tell it which key it needs to learn
 learner = adaptive.DataSaver(_learner, arg_picker=itemgetter('y'))
 ```
 
 %% Cell type:markdown id: tags:
 
 `learner.learner` is the original learner, so `learner.learner.loss()` will call the correct loss method.
 
 %% Cell type:code id: tags:
 
 ``` 
 runner = adaptive.Runner(learner, goal=lambda l: l.learner.loss() < 0.05)
 runner.live_info()
 ```
 
 %% Cell type:code id: tags:
 
 ``` 
 runner.live_plot(plotter=lambda l: l.learner.plot(), update_interval=0.1)
 ```
 
 %% Cell type:markdown id: tags:
 
 Now the `DataSavingLearner` will have an dictionary attribute `extra_data` that has `x` as key and the data that was returned by `learner.function` as values.
 
 %% Cell type:code id: tags:
 
 ``` 
 learner.extra_data
 ```
 
 %% Cell type:markdown id: tags:
 
 # `Scikit-Optimize`
 
 %% Cell type:markdown id: tags:
 
 We have wrapped the `Optimizer` class from [`scikit-optimize`](https://github.com/scikit-optimize/scikit-optimize), to show how existing libraries can be integrated with `adaptive`.
 
 The `SKOptLearner` attempts to "optimize" the given function `g` (i.e. find the global minimum of `g` in the window of interest).
 
 Here we use the same example as in the `scikit-optimize` [tutorial](https://github.com/scikit-optimize/scikit-optimize/blob/master/examples/ask-and-tell.ipynb). Although `SKOptLearner` can optimize functions of arbitrary dimensionality, we can only plot the learner if a 1D function is being learned.
 
 %% Cell type:code id: tags:
 
 ``` 
 def g(x, noise_level=0.1):
     return (np.sin(5 * x) * (1 - np.tanh(x ** 2))
             + np.random.randn() * noise_level)
 ```
 
 %% Cell type:code id: tags:
 
 ``` 
 learner = adaptive.SKOptLearner(g, dimensions=[(-2., 2.)],
                                 base_estimator="GP",
                                 acq_func="gp_hedge",
                                 acq_optimizer="lbfgs",
                                )
 runner = adaptive.Runner(learner, ntasks=1, goal=lambda l: l.npoints > 40)
 runner.live_info()
 ```
 
 %% Cell type:code id: tags:
 
 ``` 
 %%opts Overlay [legend_position='top']
 xs = np.linspace(*learner.space.bounds[0])
 to_learn = hv.Curve((xs, [g(x, 0) for x in xs]), label='to learn')
 
 runner.live_plot().relabel('prediction', depth=2) * to_learn
 ```
 
 %% Cell type:markdown id: tags:
 
 # Using multiple cores
 
 %% Cell type:markdown id: tags:
 
 Often you will want to evaluate the function on some remote computing resources. `adaptive` works out of the box with any framework that implements a [PEP 3148](https://www.python.org/dev/peps/pep-3148/) compliant executor that returns `concurrent.futures.Future` objects.
 
 %% Cell type:markdown id: tags:
 
 ### [`concurrent.futures`](https://docs.python.org/3/library/concurrent.futures.html)
 
 %% Cell type:markdown id: tags:
 
 On Unix-like systems by default `adaptive.Runner` creates a `ProcessPoolExecutor`, but you can also pass one explicitly e.g. to limit the number of workers:
 
 %% Cell type:code id: tags:
 
 ``` 
 from concurrent.futures import ProcessPoolExecutor
 
 executor = ProcessPoolExecutor(max_workers=4)
 
 learner = adaptive.Learner1D(f, bounds=(-1, 1))
 runner = adaptive.Runner(learner, executor=executor, goal=lambda l: l.loss() < 0.05)
 runner.live_info()
 runner.live_plot(update_interval=0.1)
 ```
 
 %% Cell type:markdown id: tags:
 
 ### [`ipyparallel`](https://ipyparallel.readthedocs.io/en/latest/intro.html)
 
 %% Cell type:code id: tags:
 
 ``` 
 import ipyparallel
 
 client = ipyparallel.Client()  # You will need to start an `ipcluster` to make this work
 
 learner = adaptive.Learner1D(f, bounds=(-1, 1))
 runner = adaptive.Runner(learner, executor=client, goal=lambda l: l.loss() < 0.01)
 runner.live_info()
 runner.live_plot()
 ```
 
 %% Cell type:markdown id: tags:
 
 ### [`distributed`](https://distributed.readthedocs.io/en/latest/)
 
 On Windows by default `adaptive.Runner` uses a `distributed.Client`.
 
 %% Cell type:code id: tags:
 
 ``` 
 import distributed
 
 client = distributed.Client()
 
 learner = adaptive.Learner1D(f, bounds=(-1, 1))
 runner = adaptive.Runner(learner, executor=client, goal=lambda l: l.loss() < 0.01)
 runner.live_info()
 runner.live_plot(update_interval=0.1)
 ```
 
 %% Cell type:markdown id: tags:
 
 ---
 
 %% Cell type:markdown id: tags:
 
 # Advanced Topics
 
 %% Cell type:markdown id: tags:
 
 ## A watched pot never boils!
 
 %% Cell type:markdown id: tags:
 
 `adaptive.Runner` does its work in an `asyncio` task that runs concurrently with the IPython kernel, when using `adaptive` from a Jupyter notebook. This is advantageous because it allows us to do things like live-updating plots, however it can trip you up if you're not careful.
 
 Notably: **if you block the IPython kernel, the runner will not do any work**.
 
 For example if you wanted to wait for a runner to complete, **do not wait in a busy loop**:
 ```python
 while not runner.task.done():
     pass
 ```
 
 If you do this then **the runner will never finish**.
 
 %% Cell type:markdown id: tags:
 
 What to do if you don't care about live plotting, and just want to run something until its done?
 
 The simplest way to accomplish this is to use `adaptive.BlockingRunner`:
 
 %% Cell type:code id: tags:
 
 ``` 
 learner = adaptive.Learner1D(partial(f, wait=False), bounds=(-1, 1))
 adaptive.BlockingRunner(learner, goal=lambda l: l.loss() < 0.005)
 # This will only get run after the runner has finished
 learner.plot()
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Reproducibility
 
 %% Cell type:markdown id: tags:
 
 By default `adaptive` runners evaluate the learned function in parallel across several cores. The runners are also opportunistic, in that as soon as a result is available they will feed it to the learner and request another point to replace the one that just finished.
 
 Because the order in which computations complete is non-deterministic, this means that the runner behaves in a non-deterministic way. Adaptive makes this choice because in many cases the speedup from parallel execution is worth sacrificing the "purity" of exactly reproducible computations.
 
 Nevertheless it is still possible to run a learner in a deterministic way with adaptive.
 
 The simplest way is to use `adaptive.runner.simple` to run your learner:
 
 %% Cell type:code id: tags:
 
 ``` 
 learner = adaptive.Learner1D(partial(f, wait=False), bounds=(-1, 1))
 
 # blocks until completion
 adaptive.runner.simple(learner, goal=lambda l: l.loss() < 0.002)
 
 learner.plot()
 ```
 
 %% Cell type:markdown id: tags:
 
 Note that unlike `adaptive.Runner`, `adaptive.runner.simple` *blocks* until it is finished.
 
 If you want to enable determinism, want to continue using the non-blocking `adaptive.Runner`, you can use the `adaptive.runner.SequentialExecutor`:
 
 %% Cell type:code id: tags:
 
 ``` 
 from adaptive.runner import SequentialExecutor
 
 learner = adaptive.Learner1D(f, bounds=(-1, 1))
 
 # blocks until completion
 runner = adaptive.Runner(learner, executor=SequentialExecutor(), goal=lambda l: l.loss() < 0.002)
 runner.live_info()
 runner.live_plot(update_interval=0.1)
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Cancelling a runner
 
 %% Cell type:markdown id: tags:
 
 Sometimes you want to interactively explore a parameter space, and want the function to be evaluated at finer and finer resolution and manually control when the calculation stops.
 
 If no `goal` is provided to a runner then the runner will run until cancelled.
 
 `runner.live_info()` will provide a button that can be clicked to stop the runner. You can also stop the runner programatically using `runner.cancel()`.
 
 %% Cell type:code id: tags:
 
 ``` 
 learner = adaptive.Learner1D(f, bounds=(-1, 1))
 runner = adaptive.Runner(learner)
 runner.live_info()
 runner.live_plot(update_interval=0.1)
 ```
 
 %% Cell type:code id: tags:
 
 ``` 
 runner.cancel()
 ```
 
 %% Cell type:code id: tags:
 
 ``` 
 print(runner.status())
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Debugging Problems
 
 %% Cell type:markdown id: tags:
 
 Runners work in the background with respect to the IPython kernel, which makes it convenient, but also means that inspecting errors is more difficult because exceptions will not be raised directly in the notebook. Often the only indication you will have that something has gone wrong is that nothing will be happening.
 
 Let's look at the following example, where the function to be learned will raise an exception 10% of the time.
 
 %% Cell type:code id: tags:
 
 ``` 
 def will_raise(x):
     from random import random
     from time import sleep
 
     sleep(random())
     if random() < 0.1:
         raise RuntimeError('something went wrong!')
     return x**2
 
 learner = adaptive.Learner1D(will_raise, (-1, 1))
 runner = adaptive.Runner(learner)  # without 'goal' the runner will run forever unless cancelled
 runner.live_info()
 runner.live_plot()
 ```
 
 %% Cell type:markdown id: tags:
 
 The above runner should continue forever, but we notice that it stops after a few points are evaluated.
 
 First we should check that the runner has really finished:
 
 %% Cell type:code id: tags:
 
 ``` 
 runner.task.done()
 ```
 
 %% Cell type:markdown id: tags:
 
 If it has indeed finished then we should check the `result` of the runner. This should be `None` if the runner stopped successfully. If the runner stopped due to an exception then asking for the result will raise the exception with the stack trace:
 
 %% Cell type:code id: tags:
 
 ``` 
 runner.task.result()
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Logging runners
 
 %% Cell type:markdown id: tags:
 
 Runners do their job in the background, which makes introspection quite cumbersome. One way to inspect runners is to instantiate one with `log=True`:
 
 %% Cell type:code id: tags:
 
 ``` 
 learner = adaptive.Learner1D(f, bounds=(-1, 1))
 runner = adaptive.Runner(learner, goal=lambda l: l.loss() < 0.1,
                          log=True)
 runner.live_info()
 ```
 
 %% Cell type:markdown id: tags:
 
 This gives a the runner a `log` attribute, which is a list of the `learner` methods that were called, as well as their arguments. This is useful because executors typically execute their tasks in a non-deterministic order.
 
 This can be used with `adaptive.runner.replay_log` to perfom the same set of operations on another runner:
 
 %% Cell type:code id: tags:
 
 ``` 
 reconstructed_learner = adaptive.Learner1D(f, bounds=learner.bounds)
 adaptive.runner.replay_log(reconstructed_learner, runner.log)
 ```
 
 %% Cell type:code id: tags:
 
 ``` 
 learner.plot().Scatter.I.opts(style=dict(size=6)) * reconstructed_learner.plot()
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Timing functions
 
 %% Cell type:markdown id: tags:
 
 To time the runner you **cannot** simply use
 ```python
 now = datetime.now()
 runner = adaptive.Runner(...)
 print(datetime.now() - now)
 ```
 because this will be done immediately. Also blocking the kernel with `while not runner.task.done()` will not work because the runner will not do anything when the kernel is blocked.
 
 Therefore you need to create an `async` function and hook it into the `ioloop` like so:
 
 %% Cell type:code id: tags:
 
 ``` 
 import asyncio
 
 async def time(runner):
     from datetime import datetime
     now = datetime.now()
     await runner.task
     return datetime.now() - now
 
 ioloop = asyncio.get_event_loop()
 
 learner = adaptive.Learner1D(f, bounds=(-1, 1))
 runner = adaptive.Runner(learner, goal=lambda l: l.loss() < 0.01)
 
 timer = ioloop.create_task(time(runner))
 ```
 
 %% Cell type:code id: tags:
 
 ``` 
 # The result will only be set when the runner is done.
 timer.result()
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Using Runners from a script
 
 %% Cell type:markdown id: tags:
 
 Runners can also be used from a Python script independently of the notebook.
 
 The simplest way to accomplish this is simply to use the `BlockingRunner`:
 
 ```python
 import adaptive
 
 def f(x):
     return x
 
 learner = adaptive.Learner1D(f, (-1, 1))
 
 adaptive.BlockingRunner(learner, goal=lambda: l: l.loss() < 0.1)
 ```
 
 If you use `asyncio` already in your script and want to integrate `adaptive` into it, then you can use the default `Runner` as you would from a notebook. If you want to wait for the runner to finish, then you can simply
 ```python
     await runner.task
 ```
 from within a coroutine.