Browse code

Merge branch 'cquad' into 'master'

cquad

Closes #5

See merge request qt/adaptive!8

Joseph Weston authored on 31/10/2017 15:21:22
Showing 14 changed files
1 1
deleted file mode 100644
... ...
@@ -1,774 +0,0 @@
1
-# -*- coding: utf-8 -*-
2
-import abc
3
-from contextlib import contextmanager
4
-from copy import deepcopy as copy
5
-import functools
6
-import heapq
7
-import itertools
8
-from math import sqrt, hypot
9
-from operator import itemgetter
10
-
11
-import holoviews as hv
12
-import numpy as np
13
-from scipy import interpolate, optimize, special
14
-import sortedcontainers
15
-
16
-
17
-class BaseLearner(metaclass=abc.ABCMeta):
18
-    """Base class for algorithms for learning a function 'f: X → Y'.
19
-
20
-    Attributes
21
-    ----------
22
-    function : callable: X → Y
23
-        The function to learn.
24
-    data : dict: X → Y
25
-        'function' evaluated at certain points.
26
-        The values can be 'None', which indicates that the point
27
-        will be evaluated, but that we do not have the result yet.
28
-
29
-    Subclasses may define a 'plot' method that takes no parameters
30
-    and returns a holoviews plot.
31
-    """
32
-
33
-    def add_data(self, xvalues, yvalues):
34
-        """Add data to the learner.
35
-
36
-        Parameters
37
-        ----------
38
-        xvalues : value from the function domain, or iterable of such
39
-            Values from the domain of the learned function.
40
-        yvalues : value from the function image, or iterable of such
41
-            Values from the range of the learned function, or None.
42
-            If 'None', then it indicates that the value has not yet
43
-            been computed.
44
-        """
45
-        try:
46
-            for x, y in zip(xvalues, yvalues):
47
-                self.add_point(x, y)
48
-        except TypeError:
49
-            self.add_point(xvalues, yvalues)
50
-
51
-    @abc.abstractmethod
52
-    def add_point(self, x, y):
53
-        """Add a single datapoint to the learner."""
54
-        pass
55
-
56
-    @abc.abstractmethod
57
-    def remove_unfinished(self):
58
-        """Remove uncomputed data from the learner."""
59
-        pass
60
-
61
-    @abc.abstractmethod
62
-    def loss(self, real=True):
63
-        """Return the loss for the current state of the learner.
64
-
65
-        Parameters
66
-        ----------
67
-        real : bool, default: True
68
-            If False, return the "expected" loss, i.e. the
69
-            loss including the as-yet unevaluated points
70
-            (possibly by interpolation).
71
-        """
72
-
73
-    @abc.abstractmethod
74
-    def choose_points(self, n, add_data=True):
75
-        """Choose the next 'n' points to evaluate.
76
-
77
-        Parameters
78
-        ----------
79
-        n : int
80
-            The number of points to choose.
81
-        add_data : bool, default: True
82
-            If True, add the chosen points to this
83
-            learner's 'data' with 'None' for the 'y'
84
-            values. Set this to False if you do not
85
-            want to modify the state of the learner.
86
-        """
87
-        pass
88
-
89
-    def __getstate__(self):
90
-        return copy(self.__dict__)
91
-
92
-    def __setstate__(self, state):
93
-        self.__dict__ = state
94
-
95
-
96
-class AverageLearner(BaseLearner):
97
-    """A naive implementation of adaptive computing of averages.
98
-
99
-    The learned function must depend on an integer input variable that
100
-    represents the source of randomness.
101
-
102
-    Parameters:
103
-    -----------
104
-    atol : float
105
-        Desired absolute tolerance
106
-    rtol : float
107
-        Desired relative tolerance
108
-    """
109
-
110
-    def __init__(self, function, atol=None, rtol=None):
111
-        if atol is None and rtol is None:
112
-            raise Exception('At least one of `atol` and `rtol` should be set.')
113
-        if atol is None:
114
-            atol = np.inf
115
-        if rtol is None:
116
-            rtol = np.inf
117
-
118
-        self.data = {}
119
-        self.function = function
120
-        self.atol = atol
121
-        self.rtol = rtol
122
-        self.n = 0
123
-        self.n_requested = 0
124
-        self.sum_f = 0
125
-        self.sum_f_sq = 0
126
-
127
-    def choose_points(self, n, add_data=True):
128
-        points = list(range(self.n_requested, self.n_requested + n))
129
-        loss_improvements = [self.loss()] * n
130
-        if add_data:
131
-            self.add_data(points, itertools.repeat(None))
132
-        return points, loss_improvements
133
-
134
-    def add_point(self, n, value):
135
-        self.data[n] = value
136
-        if value is None:
137
-            self.n_requested += 1
138
-            return
139
-        else:
140
-            self.n += 1
141
-            self.sum_f += value
142
-            self.sum_f_sq += value**2
143
-
144
-    @property
145
-    def mean(self):
146
-        return self.sum_f / self.n
147
-
148
-    @property
149
-    def std(self):
150
-        n = self.n
151
-        if n < 2:
152
-            return np.inf
153
-        return sqrt((self.sum_f_sq - n * self.mean**2) / (n - 1))
154
-
155
-    def loss(self, real=True):
156
-        n = self.n
157
-        if n < 2:
158
-            return np.inf
159
-        standard_error = self.std / sqrt(n if real else self.n_requested)
160
-        return max(standard_error / self.atol,
161
-                   standard_error / abs(self.mean) / self.rtol)
162
-
163
-    def remove_unfinished(self):
164
-        """Remove uncomputed data from the learner."""
165
-        pass
166
-
167
-    def plot(self):
168
-        vals = [v for v in self.data.values() if v is not None]
169
-        if not vals:
170
-            return hv.Histogram([[], []])
171
-        num_bins = int(max(5, sqrt(self.n)))
172
-        vals = hv.Points(vals)
173
-        return hv.operation.histogram(vals, num_bins=num_bins, dimension=1)
174
-
175
-
176
-class Learner1D(BaseLearner):
177
-    """Learns and predicts a function 'f:ℝ → ℝ'.
178
-
179
-    Parameters
180
-    ----------
181
-    function : callable
182
-        The function to learn. Must take a single real parameter and
183
-        return a real number.
184
-    bounds : pair of reals
185
-        The bounds of the interval on which to learn 'function'.
186
-    """
187
-
188
-    def __init__(self, function, bounds):
189
-        self.function = function
190
-
191
-        # A dict storing the loss function for each interval x_n.
192
-        self.losses = {}
193
-        self.losses_combined = {}
194
-
195
-        self.data = sortedcontainers.SortedDict()
196
-        self.data_interp = {}
197
-
198
-        # A dict {x_n: [x_{n-1}, x_{n+1}]} for quick checking of local
199
-        # properties.
200
-        self.neighbors = sortedcontainers.SortedDict()
201
-        self.neighbors_combined = sortedcontainers.SortedDict()
202
-
203
-        # Bounding box [[minx, maxx], [miny, maxy]].
204
-        self._bbox = [list(bounds), [np.inf, -np.inf]]
205
-
206
-        # Data scale (maxx - minx), (maxy - miny)
207
-        self._scale = [bounds[1] - bounds[0], 0]
208
-        self._oldscale = copy(self._scale)
209
-
210
-        self.bounds = list(bounds)
211
-
212
-    @property
213
-    def data_combined(self):
214
-        return {**self.data, **self.data_interp}
215
-
216
-    def interval_loss(self, x_left, x_right, data):
217
-        """Calculate loss in the interval x_left, x_right.
218
-
219
-        Currently returns the rescaled length of the interval. If one of the
220
-        y-values is missing, returns 0 (so the intervals with missing data are
221
-        never touched. This behavior should be improved later.
222
-        """
223
-        y_right, y_left = data[x_right], data[x_left]
224
-        if self._scale[1] == 0:
225
-            return sqrt(((x_right - x_left) / self._scale[0])**2)
226
-        else:
227
-            return sqrt(((x_right - x_left) / self._scale[0])**2 +
228
-                        ((y_right - y_left) / self._scale[1])**2)
229
-
230
-    def loss(self, real=True):
231
-        losses = self.losses if real else self.losses_combined
232
-        if len(losses) == 0:
233
-            return float('inf')
234
-        else:
235
-            return max(losses.values())
236
-
237
-    def update_losses(self, x, data, neighbors, losses):
238
-        x_lower, x_upper = neighbors[x]
239
-        if x_lower is not None:
240
-            losses[x_lower, x] = self.interval_loss(x_lower, x, data)
241
-        if x_upper is not None:
242
-            losses[x, x_upper] = self.interval_loss(x, x_upper, data)
243
-        try:
244
-            del losses[x_lower, x_upper]
245
-        except KeyError:
246
-            pass
247
-
248
-    def find_neighbors(self, x, neighbors):
249
-        pos = neighbors.bisect_left(x)
250
-        x_lower = neighbors.iloc[pos-1] if pos != 0 else None
251
-        x_upper = neighbors.iloc[pos] if pos != len(neighbors) else None
252
-        return x_lower, x_upper
253
-
254
-    def update_neighbors(self, x, neighbors):
255
-        if x not in neighbors:  # The point is new
256
-            x_lower, x_upper = self.find_neighbors(x, neighbors)
257
-            neighbors[x] = [x_lower, x_upper]
258
-            neighbors.get(x_lower, [None, None])[1] = x
259
-            neighbors.get(x_upper, [None, None])[0] = x
260
-
261
-    def update_scale(self, x, y):
262
-        self._bbox[0][0] = min(self._bbox[0][0], x)
263
-        self._bbox[0][1] = max(self._bbox[0][1], x)
264
-        if y is not None:
265
-            self._bbox[1][0] = min(self._bbox[1][0], y)
266
-            self._bbox[1][1] = max(self._bbox[1][1], y)
267
-
268
-        self._scale = [self._bbox[0][1] - self._bbox[0][0],
269
-                       self._bbox[1][1] - self._bbox[1][0]]
270
-
271
-    def add_point(self, x, y):
272
-        real = y is not None
273
-
274
-        if real:
275
-            # Add point to the real data dict and pop from the unfinished
276
-            # data_interp dict.
277
-            self.data[x] = y
278
-            try:
279
-                del self.data_interp[x]
280
-            except KeyError:
281
-                pass
282
-        else:
283
-            # The keys of data_interp are the unknown points
284
-            self.data_interp[x] = None
285
-
286
-        # Update the neighbors
287
-        self.update_neighbors(x, self.neighbors_combined)
288
-        if real:
289
-            self.update_neighbors(x, self.neighbors)
290
-
291
-        # Update the scale
292
-        self.update_scale(x, y)
293
-
294
-        # Interpolate
295
-        if not real:
296
-            self.data_interp = self.interpolate()
297
-
298
-        # Update the losses
299
-        self.update_losses(x, self.data_combined, self.neighbors_combined,
300
-                           self.losses_combined)
301
-        if real:
302
-            self.update_losses(x, self.data, self.neighbors, self.losses)
303
-
304
-        # If the scale has doubled, recompute all losses.
305
-        if self._scale > self._oldscale * 2:
306
-            self.losses = {xs: self.interval_loss(*xs, self.data)
307
-                           for xs in self.losses}
308
-            self.losses_combined = {x: self.interval_loss(*x,
309
-                                                          self.data_combined)
310
-                                    for x in self.losses_combined}
311
-            self._oldscale = self._scale
312
-
313
-    def choose_points(self, n, add_data=True):
314
-        """Return n points that are expected to maximally reduce the loss."""
315
-        # Find out how to divide the n points over the intervals
316
-        # by finding  positive integer n_i that minimize max(L_i / n_i) subject
317
-        # to a constraint that sum(n_i) = n + N, with N the total number of
318
-        # intervals.
319
-
320
-        # Return equally spaced points within each interval to which points
321
-        # will be added.
322
-        if n == 0:
323
-            return []
324
-
325
-        # If the bounds have not been chosen yet, we choose them first.
326
-        points = []
327
-        for bound in self.bounds:
328
-            if bound not in self.data and bound not in self.data_interp:
329
-                points.append(bound)
330
-
331
-        # Ensure we return exactly 'n' points.
332
-        if points:
333
-            loss_improvements = [float('inf')] * n
334
-            if n <= 2:
335
-                points = points[:n]
336
-            else:
337
-                points = np.linspace(*self.bounds, n)
338
-        else:
339
-            def xs(x, n):
340
-                if n == 1:
341
-                    return []
342
-                else:
343
-                    step = (x[1] - x[0]) / n
344
-                    return [x[0] + step * i for i in range(1, n)]
345
-
346
-            # Calculate how many points belong to each interval.
347
-            quals = [(-loss, x_range, 1) for (x_range, loss) in
348
-                     self.losses_combined.items()]
349
-
350
-            heapq.heapify(quals)
351
-
352
-            for point_number in range(n):
353
-                quality, x, n = quals[0]
354
-                heapq.heapreplace(quals, (quality * n / (n + 1), x, n + 1))
355
-
356
-            points = list(itertools.chain.from_iterable(xs(x, n)
357
-                          for quality, x, n in quals))
358
-
359
-            loss_improvements = list(itertools.chain.from_iterable(
360
-                                     itertools.repeat(-quality, n)
361
-                                     for quality, x, n in quals))
362
-
363
-        if add_data:
364
-            self.add_data(points, itertools.repeat(None))
365
-
366
-        return points, loss_improvements
367
-
368
-    def interpolate(self, extra_points=None):
369
-        xs = list(self.data.keys())
370
-        ys = list(self.data.values())
371
-        xs_unfinished = list(self.data_interp.keys())
372
-
373
-        if extra_points is not None:
374
-            xs_unfinished += extra_points
375
-
376
-        if len(ys) == 0:
377
-            interp_ys = (0,) * len(xs_unfinished)
378
-        else:
379
-            interp_ys = np.interp(xs_unfinished, xs, ys)
380
-
381
-        data_interp = {x: y for x, y in zip(xs_unfinished, interp_ys)}
382
-
383
-        return data_interp
384
-
385
-    def plot(self):
386
-            if self.data:
387
-                return hv.Scatter(self.data)
388
-            else:
389
-                return hv.Scatter([])
390
-
391
-    def remove_unfinished(self):
392
-        self.data_interp = {}
393
-        self.losses_combined = copy(self.losses)
394
-        self.neighbors_combined = copy(self.neighbors)
395
-
396
-
397
-def dispatch(child_functions, arg):
398
-    index, x = arg
399
-    return child_functions[index](x)
400
-
401
-
402
-class BalancingLearner(BaseLearner):
403
-    """Choose the optimal points from a set of learners.
404
-
405
-    Parameters
406
-    ----------
407
-    learners : sequence of BaseLearner
408
-        The learners from which to choose. These must all have the same type.
409
-
410
-    Notes
411
-    -----
412
-    This learner compares the 'loss' calculated from the "child" learners.
413
-    This requires that the 'loss' from different learners *can be meaningfully
414
-    compared*. For the moment we enforce this restriction by requiring that
415
-    all learners are the same type but (depending on the internals of the
416
-    learner) it may be that the loss cannot be compared *even between learners
417
-    of the same type*. In this case the BalancingLearner will behave in an
418
-    undefined way.
419
-    """
420
-
421
-    def __init__(self, learners):
422
-        self.learners = learners
423
-
424
-        # Naively we would make 'function' a method, but this causes problems
425
-        # when using executors from 'concurrent.futures' because we have to
426
-        # pickle the whole learner.
427
-        self.function = functools.partial(dispatch, [l.function for l
428
-                                                     in self.learners])
429
-
430
-        if len(set(learner.__class__ for learner in self.learners)) > 1:
431
-            raise TypeError('A BalacingLearner can handle only one type'
432
-                            'of learners.')
433
-
434
-    def _choose_and_add_points(self, n):
435
-        points = []
436
-        for _ in range(n):
437
-            loss_improvements = []
438
-            pairs = []
439
-            for index, learner in enumerate(self.learners):
440
-                point, loss_improvement = learner.choose_points(n=1,
441
-                                                                add_data=False)
442
-                loss_improvements.append(loss_improvement[0])
443
-                pairs.append((index, point[0]))
444
-            x, _ = max(zip(pairs, loss_improvements), key=itemgetter(1))
445
-            points.append(x)
446
-            self.add_point(x, None)
447
-        return points, None
448
-
449
-    def choose_points(self, n, add_data=True):
450
-        """Chose points for learners."""
451
-        if not add_data:
452
-            with restore(*self.learners):
453
-                return self._choose_and_add_points(n)
454
-        else:
455
-            return self._choose_and_add_points(n)
456
-
457
-    def add_point(self, x, y):
458
-        index, x = x
459
-        self.learners[index].add_point(x, y)
460
-
461
-    def loss(self, real=True):
462
-        return max(learner.loss(real) for learner in self.learners)
463
-
464
-    def plot(self, index):
465
-        return self.learners[index].plot()
466
-
467
-    def remove_unfinished(self):
468
-        """Remove uncomputed data from the learners."""
469
-        for learner in self.learners:
470
-            learner.remove_unfinished()
471
-
472
-
473
-# Learner2D and helper functions.
474
-
475
-def _losses_per_triangle(ip):
476
-    tri = ip.tri
477
-    vs = ip.values.ravel()
478
-
479
-    gradients = interpolate.interpnd.estimate_gradients_2d_global(
480
-        tri, vs, tol=1e-6)
481
-    p = tri.points[tri.vertices]
482
-    g = gradients[tri.vertices]
483
-    v = vs[tri.vertices]
484
-    n_points_per_triangle = p.shape[1]
485
-
486
-    dev = 0
487
-    for j in range(n_points_per_triangle):
488
-        vest = v[:, j, None] + ((p[:, :, :] - p[:, j, None, :]) *
489
-                                g[:, j, None, :]).sum(axis=-1)
490
-        dev += abs(vest - v).max(axis=1)
491
-
492
-    q = p[:, :-1, :] - p[:, -1, None, :]
493
-    areas = abs(q[:, 0, 0] * q[:, 1, 1] - q[:, 0, 1] * q[:, 1, 0])
494
-    areas /= special.gamma(n_points_per_triangle)
495
-    areas = np.sqrt(areas)
496
-
497
-    vs_scale = vs[tri.vertices].ptp()
498
-    if vs_scale != 0:
499
-        dev /= vs_scale
500
-
501
-    return dev * areas
502
-
503
-class Learner2D(BaseLearner):
504
-    """Learns and predicts a function 'f: ℝ^2 → ℝ'.
505
-
506
-    Parameters
507
-    ----------
508
-    function : callable
509
-        The function to learn. Must take a tuple of two real
510
-        parameters and return a real number.
511
-    bounds : list of 2-tuples
512
-        A list ``[(a1, b1), (a2, b2)]`` containing bounds,
513
-        one per dimension.
514
-
515
-    Attributes
516
-    ----------
517
-    points_combined
518
-        Sample points so far including the unknown interpolated ones.
519
-    values_combined
520
-        Sampled values so far including the unknown interpolated ones.
521
-    points
522
-        Sample points so far with real results.
523
-    values
524
-        Sampled values so far with real results.
525
-
526
-    Notes
527
-    -----
528
-    Adapted from an initial implementation by Pauli Virtanen.
529
-
530
-    The sample points are chosen by estimating the point where the
531
-    linear and cubic interpolants based on the existing points have
532
-    maximal disagreement. This point is then taken as the next point
533
-    to be sampled.
534
-
535
-    In practice, this sampling protocol results to sparser sampling of
536
-    smooth regions, and denser sampling of regions where the function
537
-    changes rapidly, which is useful if the function is expensive to
538
-    compute.
539
-
540
-    This sampling procedure is not extremely fast, so to benefit from
541
-    it, your function needs to be slow enough to compute.
542
-    """
543
-
544
-    def __init__(self, function, bounds):
545
-        self.ndim = len(bounds)
546
-        if self.ndim != 2:
547
-            raise ValueError("Only 2-D sampling supported.")
548
-        self.bounds = tuple((float(a), float(b)) for a, b in bounds)
549
-        self._points = np.zeros([100, self.ndim])
550
-        self._values = np.zeros([100], dtype=float)
551
-        self._stack = []
552
-        self._interp = {}
553
-
554
-        xy_mean = np.mean(self.bounds, axis=1)
555
-        xy_scale = np.ptp(self.bounds, axis=1)
556
-
557
-        def scale(points):
558
-            return (points - xy_mean) / xy_scale
559
-
560
-        def unscale(points):
561
-            return points * xy_scale + xy_mean
562
-
563
-        self.scale = scale
564
-        self.unscale = unscale
565
-
566
-        # Keeps track till which index _points and _values are filled
567
-        self.n = 0
568
-
569
-        self._bounds_points = list(itertools.product(*bounds))
570
-
571
-        # Add the loss improvement to the bounds in the stack
572
-        self._stack = [(*p, np.inf) for p in self._bounds_points]
573
-
574
-        self.function = function
575
-
576
-    @property
577
-    def points_combined(self):
578
-        return self._points[:self.n]
579
-
580
-    @property
581
-    def values_combined(self):
582
-        return self._values[:self.n]
583
-
584
-    @property
585
-    def points(self):
586
-        return np.delete(self.points_combined,
587
-                         list(self._interp.values()), axis=0)
588
-
589
-    @property
590
-    def values(self):
591
-        return np.delete(self.values_combined,
592
-                         list(self._interp.values()), axis=0)
593
-
594
-    def ip(self):
595
-        points = self.scale(self.points)
596
-        return interpolate.LinearNDInterpolator(points, self.values)
597
-
598
-    @property
599
-    def n_real(self):
600
-        return self.n - len(self._interp)
601
-
602
-    def ip_combined(self):
603
-        points = self.scale(self.points_combined)
604
-        values = self.values_combined
605
-
606
-        # Interpolate the unfinished points
607
-        if self._interp:
608
-            n_interp = list(self._interp.values())
609
-            bounds_are_done = not any(p in self._interp
610
-                                      for p in self._bounds_points)
611
-            if bounds_are_done:
612
-                values[n_interp] = self.ip()(points[n_interp])
613
-            else:
614
-                # It is important not to return exact zeros because
615
-                # otherwise the algo will try to add the same point
616
-                # to the stack each time.
617
-                values[n_interp] = np.random.rand(len(n_interp)) * 1e-15
618
-
619
-        return interpolate.LinearNDInterpolator(points, values)
620
-
621
-    def add_point(self, point, value):
622
-        nmax = self.values_combined.shape[0]
623
-        if self.n >= nmax:
624
-            self._values = np.resize(self._values, [2*nmax + 10])
625
-            self._points = np.resize(self._points, [2*nmax + 10, self.ndim])
626
-
627
-        point = tuple(point)
628
-
629
-        # When the point is not evaluated yet, add an entry to self._interp
630
-        # that saves the point and index.
631
-        if value is None:
632
-            self._interp[point] = self.n
633
-            old_point = False
634
-        else:
635
-            old_point = point in self._interp
636
-
637
-        # If the point is new add it a new value to _points and _values,
638
-        # otherwise get the index of the value that is being replaced.
639
-        if old_point:
640
-            n = self._interp.pop(point)
641
-        else:
642
-            n = self.n
643
-            self.n += 1
644
-
645
-        self._points[n] = point
646
-        self._values[n] = value
647
-
648
-        # Remove the point if in the stack.
649
-        for i, (*_point, _) in enumerate(self._stack):
650
-            if point == tuple(_point):
651
-                self._stack.pop(i)
652
-                break
653
-
654
-    def _fill_stack(self, stack_till=None):
655
-        if stack_till is None:
656
-            stack_till = 1
657
-
658
-        if self.values_combined.shape[0] < self.ndim + 1:
659
-            raise ValueError("too few points...")
660
-
661
-        # Interpolate
662
-        ip = self.ip_combined()
663
-        tri = ip.tri
664
-
665
-        losses = _losses_per_triangle(ip)
666
-
667
-        def point_exists(p):
668
-            eps = np.finfo(float).eps * self.points_combined.ptp() * 100
669
-            if abs(p - self.points_combined).sum(axis=1).min() < eps:
670
-                return True
671
-            if self._stack:
672
-                _stack_points, _ = self._split_stack()
673
-                if abs(p - np.asarray(_stack_points)).sum(axis=1).min() < eps:
674
-                    return True
675
-            return False
676
-
677
-        for j, _ in enumerate(losses):
678
-            # Estimate point of maximum curvature inside the simplex
679
-            jsimplex = np.argmax(losses)
680
-            p = tri.points[tri.vertices[jsimplex]]
681
-            point_new = self.unscale(p.mean(axis=-2))
682
-
683
-            # XXX: not sure whether this is necessary it was there
684
-            # originally.
685
-            point_new = np.clip(point_new, *zip(*self.bounds))
686
-
687
-            # Check if it is really new
688
-            if point_exists(point_new):
689
-                losses[jsimplex] = 0
690
-                continue
691
-
692
-            # Add to stack
693
-            self._stack.append((*point_new, losses[jsimplex]))
694
-
695
-            if len(self._stack) >= stack_till:
696
-                break
697
-            else:
698
-                losses[jsimplex] = 0
699
-
700
-    def _split_stack(self, n=None):
701
-        points = []
702
-        loss_improvements = []
703
-        for *point, loss_improvement in self._stack[:n]:
704
-            points.append(point)
705
-            loss_improvements.append(loss_improvement)
706
-        return points, loss_improvements
707
-
708
-    def _choose_and_add_points(self, n):
709
-        if n <= len(self._stack):
710
-            points, loss_improvements = self._split_stack(n)
711
-            self.add_data(points, itertools.repeat(None))
712
-        else:
713
-            points = []
714
-            loss_improvements = []
715
-            n_left = n
716
-            while n_left > 0:
717
-                # The while loop is needed because `stack_till` could be larger
718
-                # than the number of triangles between the points. Therefore
719
-                # it could fill up till a length smaller than `stack_till`.
720
-                if self.n >= 2**self.ndim:
721
-                    # Only fill the stack if no more bounds left in _stack
722
-                    self._fill_stack(stack_till=n_left)
723
-                new_points, new_loss_improvements = self._split_stack(n_left)
724
-                points += new_points
725
-                loss_improvements += new_loss_improvements
726
-                self.add_data(new_points, itertools.repeat(None))
727
-                n_left -= len(new_points)
728
-
729
-        return points, loss_improvements
730
-
731
-    def choose_points(self, n, add_data=True):
732
-        if not add_data:
733
-            with restore(self):
734
-                return self._choose_and_add_points(n)
735
-        else:
736
-            return self._choose_and_add_points(n)
737
-
738
-    def loss(self, real=True):
739
-        n = self.n_real if real else self.n
740
-        bounds_are_not_done = any(p in self._interp
741
-                                  for p in self._bounds_points)
742
-        if n <= 4 or bounds_are_not_done:
743
-            return np.inf
744
-        ip = self.ip() if real else self.ip_combined()
745
-        losses = _losses_per_triangle(ip)
746
-        return losses.max()
747
-
748
-    def remove_unfinished(self):
749
-        self._points = self.points.copy()
750
-        self._values = self.values.copy()
751
-        self.n -= len(self._interp)
752
-        self._interp = {}
753
-
754
-    def plot(self, n_x=201, n_y=201):
755
-        x, y = self.bounds
756
-        lbrt = x[0], y[0], x[1], y[1]
757
-        if self.n_real >= 4:
758
-            x = np.linspace(-0.5, 0.5, n_x)
759
-            y = np.linspace(-0.5, 0.5, n_y)
760
-            ip = self.ip()
761
-            z = ip(x[:, None], y[None, :])
762
-            return hv.Image(np.rot90(z), bounds=lbrt)
763
-        else:
764
-            return hv.Image(np.zeros((2, 2)), bounds=lbrt)
765
-
766
-
767
-@contextmanager
768
-def restore(*learners):
769
-    states = [learner.__getstate__() for learner in learners]
770
-    try:
771
-        yield
772
-    finally:
773
-        for state, learner in zip(states, learners):
774
-            learner.__setstate__(state)
775 0
new file mode 100644
... ...
@@ -0,0 +1,7 @@
1
+# -*- coding: utf-8 -*-
2
+from .average_learner import AverageLearner
3
+from .base_learner import BaseLearner
4
+from .balancing_learner import BalancingLearner
5
+from .learner1D import Learner1D
6
+from .learner2D import Learner2D
7
+from .integrator_learner import IntegratorLearner
0 8
new file mode 100644
... ...
@@ -0,0 +1,87 @@
1
+# -*- coding: utf-8 -*-
2
+import itertools
3
+from math import sqrt
4
+
5
+import holoviews as hv
6
+import numpy as np
7
+
8
+from .base_learner import BaseLearner
9
+
10
+class AverageLearner(BaseLearner):
11
+    """A naive implementation of adaptive computing of averages.
12
+
13
+    The learned function must depend on an integer input variable that
14
+    represents the source of randomness.
15
+
16
+    Parameters:
17
+    -----------
18
+    atol : float
19
+        Desired absolute tolerance
20
+    rtol : float
21
+        Desired relative tolerance
22
+    """
23
+
24
+    def __init__(self, function, atol=None, rtol=None):
25
+        if atol is None and rtol is None:
26
+            raise Exception('At least one of `atol` and `rtol` should be set.')
27
+        if atol is None:
28
+            atol = np.inf
29
+        if rtol is None:
30
+            rtol = np.inf
31
+
32
+        self.data = {}
33
+        self.function = function
34
+        self.atol = atol
35
+        self.rtol = rtol
36
+        self.n = 0
37
+        self.n_requested = 0
38
+        self.sum_f = 0
39
+        self.sum_f_sq = 0
40
+
41
+    def choose_points(self, n, add_data=True):
42
+        points = list(range(self.n_requested, self.n_requested + n))
43
+        loss_improvements = [self.loss()] * n
44
+        if add_data:
45
+            self.add_data(points, itertools.repeat(None))
46
+        return points, loss_improvements
47
+
48
+    def add_point(self, n, value):
49
+        self.data[n] = value
50
+        if value is None:
51
+            self.n_requested += 1
52
+            return
53
+        else:
54
+            self.n += 1
55
+            self.sum_f += value
56
+            self.sum_f_sq += value**2
57
+
58
+    @property
59
+    def mean(self):
60
+        return self.sum_f / self.n
61
+
62
+    @property
63
+    def std(self):
64
+        n = self.n
65
+        if n < 2:
66
+            return np.inf
67
+        return sqrt((self.sum_f_sq - n * self.mean**2) / (n - 1))
68
+
69
+    def loss(self, real=True):
70
+        n = self.n
71
+        if n < 2:
72
+            return np.inf
73
+        standard_error = self.std / sqrt(n if real else self.n_requested)
74
+        return max(standard_error / self.atol,
75
+                   standard_error / abs(self.mean) / self.rtol)
76
+
77
+    def remove_unfinished(self):
78
+        """Remove uncomputed data from the learner."""
79
+        pass
80
+
81
+    def plot(self):
82
+        vals = [v for v in self.data.values() if v is not None]
83
+        if not vals:
84
+            return hv.Histogram([[], []])
85
+        num_bins = int(max(5, sqrt(self.n)))
86
+        vals = hv.Points(vals)
87
+        return hv.operation.histogram(vals, num_bins=num_bins, dimension=1)
0 88
new file mode 100644
... ...
@@ -0,0 +1,82 @@
1
+# -*- coding: utf-8 -*-
2
+import functools
3
+from operator import itemgetter
4
+
5
+from .base_learner import BaseLearner
6
+from .utils import restore
7
+
8
+
9
+def dispatch(child_functions, arg):
10
+    index, x = arg
11
+    return child_functions[index](x)
12
+
13
+
14
+class BalancingLearner(BaseLearner):
15
+    """Choose the optimal points from a set of learners.
16
+
17
+    Parameters
18
+    ----------
19
+    learners : sequence of BaseLearner
20
+        The learners from which to choose. These must all have the same type.
21
+
22
+    Notes
23
+    -----
24
+    This learner compares the 'loss' calculated from the "child" learners.
25
+    This requires that the 'loss' from different learners *can be meaningfully
26
+    compared*. For the moment we enforce this restriction by requiring that
27
+    all learners are the same type but (depending on the internals of the
28
+    learner) it may be that the loss cannot be compared *even between learners
29
+    of the same type*. In this case the BalancingLearner will behave in an
30
+    undefined way.
31
+    """
32
+
33
+    def __init__(self, learners):
34
+        self.learners = learners
35
+
36
+        # Naively we would make 'function' a method, but this causes problems
37
+        # when using executors from 'concurrent.futures' because we have to
38
+        # pickle the whole learner.
39
+        self.function = functools.partial(dispatch, [l.function for l
40
+                                                     in self.learners])
41
+
42
+        if len(set(learner.__class__ for learner in self.learners)) > 1:
43
+            raise TypeError('A BalacingLearner can handle only one type'
44
+                            'of learners.')
45
+
46
+    def _choose_and_add_points(self, n):
47
+        points = []
48
+        for _ in range(n):
49
+            loss_improvements = []
50
+            pairs = []
51
+            for index, learner in enumerate(self.learners):
52
+                point, loss_improvement = learner.choose_points(n=1,
53
+                                                                add_data=False)
54
+                loss_improvements.append(loss_improvement[0])
55
+                pairs.append((index, point[0]))
56
+            x, _ = max(zip(pairs, loss_improvements), key=itemgetter(1))
57
+            points.append(x)
58
+            self.add_point(x, None)
59
+        return points, None
60
+
61
+    def choose_points(self, n, add_data=True):
62
+        """Chose points for learners."""
63
+        if not add_data:
64
+            with restore(*self.learners):
65
+                return self._choose_and_add_points(n)
66
+        else:
67
+            return self._choose_and_add_points(n)
68
+
69
+    def add_point(self, x, y):
70
+        index, x = x
71
+        self.learners[index].add_point(x, y)
72
+
73
+    def loss(self, real=True):
74
+        return max(learner.loss(real) for learner in self.learners)
75
+
76
+    def plot(self, index):
77
+        return self.learners[index].plot()
78
+
79
+    def remove_unfinished(self):
80
+        """Remove uncomputed data from the learners."""
81
+        for learner in self.learners:
82
+            learner.remove_unfinished()
0 83
new file mode 100644
... ...
@@ -0,0 +1,83 @@
1
+# -*- coding: utf-8 -*-
2
+import abc
3
+import collections
4
+from copy import deepcopy
5
+
6
+
7
+class BaseLearner(metaclass=abc.ABCMeta):
8
+    """Base class for algorithms for learning a function 'f: X → Y'.
9
+
10
+    Attributes
11
+    ----------
12
+    function : callable: X → Y
13
+        The function to learn.
14
+    data : dict: X → Y
15
+        'function' evaluated at certain points.
16
+        The values can be 'None', which indicates that the point
17
+        will be evaluated, but that we do not have the result yet.
18
+
19
+    Subclasses may define a 'plot' method that takes no parameters
20
+    and returns a holoviews plot.
21
+    """
22
+
23
+    def add_data(self, xvalues, yvalues):
24
+        """Add data to the learner.
25
+
26
+        Parameters
27
+        ----------
28
+        xvalues : value from the function domain, or iterable of such
29
+            Values from the domain of the learned function.
30
+        yvalues : value from the function image, or iterable of such
31
+            Values from the range of the learned function, or None.
32
+            If 'None', then it indicates that the value has not yet
33
+            been computed.
34
+        """
35
+        if all(isinstance(i, collections.Iterable) for i in [xvalues, yvalues]):
36
+            for x, y in zip(xvalues, yvalues):
37
+                self.add_point(x, y)
38
+        else:
39
+            self.add_point(xvalues, yvalues)
40
+
41
+    @abc.abstractmethod
42
+    def add_point(self, x, y):
43
+        """Add a single datapoint to the learner."""
44
+        pass
45
+
46
+    @abc.abstractmethod
47
+    def remove_unfinished(self):
48
+        """Remove uncomputed data from the learner."""
49
+        pass
50
+
51
+    @abc.abstractmethod
52
+    def loss(self, real=True):
53
+        """Return the loss for the current state of the learner.
54
+
55
+        Parameters
56
+        ----------
57
+        real : bool, default: True
58
+            If False, return the "expected" loss, i.e. the
59
+            loss including the as-yet unevaluated points
60
+            (possibly by interpolation).
61
+        """
62
+
63
+    @abc.abstractmethod
64
+    def choose_points(self, n, add_data=True):
65
+        """Choose the next 'n' points to evaluate.
66
+
67
+        Parameters
68
+        ----------
69
+        n : int
70
+            The number of points to choose.
71
+        add_data : bool, default: True
72
+            If True, add the chosen points to this
73
+            learner's 'data' with 'None' for the 'y'
74
+            values. Set this to False if you do not
75
+            want to modify the state of the learner.
76
+        """
77
+        pass
78
+
79
+    def __getstate__(self):
80
+        return deepcopy(self.__dict__)
81
+
82
+    def __setstate__(self, state):
83
+        self.__dict__ = state
0 84
new file mode 100644
... ...
@@ -0,0 +1,160 @@
1
+# -*- coding: utf-8 -*-
2
+from fractions import Fraction
3
+from collections import defaultdict
4
+import numpy as np
5
+import scipy.linalg
6
+
7
+def legendre(n):
8
+    """Return the first n Legendre polynomials.
9
+
10
+    The polynomials have *standard* normalization, i.e.
11
+    int_{-1}^1 dx L_n(x) L_m(x) = delta(m, n) * 2 / (2 * n + 1).
12
+
13
+    The return value is a list of list of fraction.Fraction instances.
14
+    """
15
+    result = [[Fraction(1)], [Fraction(0), Fraction(1)]]
16
+    if n <= 2:
17
+        return result[:n]
18
+    for i in range(2, n):
19
+        # Use Bonnet's recursion formula.
20
+        new = (i + 1) * [Fraction(0)]
21
+        new[1:] = (r * (2*i - 1) for r in result[-1])
22
+        new[:-2] = (n - r * (i - 1) for n, r in zip(new[:-2], result[-2]))
23
+        new[:] = (n / i for n in new)
24
+        result.append(new)
25
+    return result
26
+
27
+
28
+def newton(n):
29
+    """Compute the monomial coefficients of the Newton polynomial over the
30
+    nodes of the n-point Clenshaw-Curtis quadrature rule.
31
+    """
32
+    # The nodes of the Clenshaw-Curtis rule are x_i = -cos(i * Pi / (n-1)).
33
+    # Here, we calculate the coefficients c_i such that sum_i c_i * x^i
34
+    # = prod_i (x - x_i).  The coefficients are thus sums of products of
35
+    # cosines.
36
+    #
37
+    # This routine uses the relation
38
+    #   cos(a) cos(b) = (cos(a + b) + cos(a - b)) / 2
39
+    # to efficiently calculate the coefficients.
40
+    #
41
+    # The dictionary 'terms' descibes the terms that make up the
42
+    # monomial coefficients.  Each item ((d, a), m) corresponds to a
43
+    # term m * cos(a * Pi / n) to be added to prefactor of the
44
+    # monomial x^(n-d).
45
+
46
+    mod = 2 * (n-1)
47
+    terms = defaultdict(int)
48
+    terms[0, 0] += 1
49
+
50
+    for i in range(n):
51
+        newterms = []
52
+        for (d, a), m in terms.items():
53
+            for b in [i, -i]:
54
+                # In order to reduce the number of terms, cosine
55
+                # arguments are mapped back to the inteval [0, pi/2).
56
+                arg = (a + b) % mod
57
+                if arg > n-1:
58
+                    arg = mod - arg
59
+                if arg >= n // 2:
60
+                    if n % 2 and arg == n // 2:
61
+                        # Zero term: ignore
62
+                        continue
63
+                    newterms.append((d + 1, n - 1 - arg, -m))
64
+                else:
65
+                    newterms.append((d + 1, arg, m))
66
+        for d, s, m in newterms:
67
+            terms[d, s] += m
68
+
69
+    c = (n + 1) * [0]
70
+    for (d, a), m in terms.items():
71
+        if m and a != 0:
72
+            raise ValueError("Newton polynomial cannot be represented exactly.")
73
+        c[n - d] += m
74
+        # The check could be removed and the above line replaced by
75
+        # the following, but then the result would be no longer exact.
76
+        # c[n - d] += m * np.cos(a * np.pi / (n - 1))
77
+
78
+    cf = np.array(c, float)
79
+    assert all(int(cfe) == ce for cfe, ce in zip(cf, c)), 'Precision loss'
80
+
81
+    cf /= 2.**np.arange(n, -1, -1)
82
+    return cf
83
+
84
+
85
+def scalar_product(a, b):
86
+    """Compute the polynomial scalar product int_-1^1 dx a(x) b(x).
87
+
88
+    The args must be sequences of polynomial coefficients.  This
89
+    function is careful to use the input data type for calculations.
90
+    """
91
+    la = len(a)
92
+    lc = len(b) + la + 1
93
+
94
+    # Compute the even coefficients of the product of a and b.
95
+    c = lc * [a[0].__class__()]
96
+    for i, bi in enumerate(b):
97
+        if bi == 0:
98
+            continue
99
+        for j in range(i % 2, la, 2):
100
+            c[i + j] += a[j] * bi
101
+
102
+    # Calculate the definite integral from -1 to 1.
103
+    return 2 * sum(c[i] / (i + 1) for i in range(0, lc, 2))
104
+
105
+
106
+def calc_bdef(ns):
107
+    """Calculate the decompositions of Newton polynomials (over the nodes
108
+    of the n-point Clenshaw-Curtis quadrature rule) in terms of
109
+    Legandre polynomials.
110
+
111
+    The parameter 'ns' is a sequence of numers of points of the
112
+    quadrature rule.  The return value is a corresponding sequence of
113
+    normalized Legendre polynomial coefficients.
114
+    """
115
+    legs = legendre(max(ns) + 1)
116
+    result = []
117
+    for n in ns:
118
+        poly = []
119
+        a = list(map(Fraction, newton(n)))
120
+        for b in legs[:n + 1]:
121
+            igral = scalar_product(a, b)
122
+
123
+            # Normalize & store.  (The polynomials returned by
124
+            # legendre() have standard normalization that is not
125
+            # orthonormal.)
126
+            poly.append(np.sqrt((2*len(b) - 1) / 2) * igral)
127
+
128
+        result.append(np.array(poly))
129
+    return result
130
+
131
+
132
+def calc_V(x, n):
133
+    V = [np.ones(x.shape), x.copy()]
134
+    for i in range(2, n):
135
+        V.append((2*i-1) / i * x * V[-1] - (i-1) / i * V[-2])
136
+    for i in range(n):
137
+        V[i] *= np.sqrt(i + 0.5)
138
+    return np.array(V).T
139
+
140
+# the nodes and Newton polynomials
141
+ns = (5, 9, 17, 33)
142
+xi = [-np.cos(np.linspace(0, np.pi, n)) for n in ns]
143
+
144
+# Make `xi` perfectly anti-symmetric, important for splitting the intervals
145
+xi = [(row - row[::-1]) / 2 for row in xi]
146
+
147
+# compute the coefficients
148
+V = [calc_V(x, n) for x, n in zip(xi, ns)]
149
+V_inv = list(map(scipy.linalg.inv, V))
150
+Vcond = list(map(np.linalg.cond, V))
151
+
152
+# shift matrix
153
+T_left, T_right = [V_inv[3] @ calc_V((xi[3] + a) / 2, ns[3]) for a in [-1, 1]]
154
+
155
+# set-up the downdate matrix
156
+k = np.arange(ns[3])
157
+alpha = np.sqrt((k+1)**2 / (2*k+1) / (2*k+3))
158
+gamma = np.concatenate([[0, 0], np.sqrt(k[2:]**2 / (4*k[2:]**2-1))])
159
+
160
+b_def = calc_bdef(ns)
0 161
new file mode 100644
... ...
@@ -0,0 +1,575 @@
1
+# -*- coding: utf-8 -*-
2
+# Copyright 2010 Pedro Gonnet
3
+# Copyright 2017 Christoph Groth
4
+# Copyright 2017 `adaptive` authors
5
+
6
+from collections import defaultdict
7
+from math import sqrt
8
+from operator import attrgetter
9
+
10
+import holoviews as hv
11
+import numpy as np
12
+from scipy.linalg import norm
13
+from sortedcontainers import SortedDict, SortedSet
14
+
15
+from .base_learner import BaseLearner
16
+from .integrator_coeffs import (b_def, T_left, T_right, ns,
17
+                                xi, V_inv, Vcond, alpha, gamma)
18
+
19
+
20
+eps = np.spacing(1)
21
+
22
+
23
+def _downdate(c, nans, depth):
24
+    b = b_def[depth].copy()
25
+    m = ns[depth] - 1
26
+    for i in nans:
27
+        b[m + 1] /= alpha[m]
28
+        xii = xi[depth][i]
29
+        b[m] = (b[m] + xii * b[m + 1]) / alpha[m - 1]
30
+        for j in range(m - 1, 0, -1):
31
+            b[j] = ((b[j] + xii * b[j + 1] - gamma[j + 1] * b[j + 2])
32
+                    / alpha[j - 1])
33
+        b = b[1:]
34
+
35
+        c[:m] -= c[m] / b[m] * b[:m]
36
+        c[m] = 0
37
+        m -= 1
38
+
39
+
40
+def _zero_nans(fx):
41
+    """Caution: this function modifies fx."""
42
+    nans = []
43
+    for i in range(len(fx)):
44
+        if not np.isfinite(fx[i]):
45
+            nans.append(i)
46
+            fx[i] = 0.0
47
+    return nans
48
+
49
+
50
+def _calc_coeffs(fx, depth):
51
+    """Caution: this function modifies fx."""
52
+    nans = _zero_nans(fx)
53
+    c_new = V_inv[depth] @ fx
54
+    if nans:
55
+        fx[nans] = np.nan
56
+        _downdate(c_new, nans, depth)
57
+    return c_new
58
+
59
+
60
+class DivergentIntegralError(ValueError):
61
+    pass
62
+
63
+
64
+class Interval:
65
+
66
+    """
67
+    Attributes
68
+    ----------
69
+    (a, b) : (float, float)
70
+        The left and right boundary of the interval.
71
+    c : numpy array of shape (4, 33)
72
+        Coefficients of the fit.
73
+    c_old : numpy array of shape `len(fx)`
74
+        Coefficients of the fit.
75
+    depth : int
76
+        The level of refinement, `depth=0` means that it has 5 (the minimal number of) points and
77
+        `depth=3` means it has 33 (the maximal number of) points.
78
+    fx : numpy array of size `(5, 9, 17, 33)[self.depth]`.
79
+        The function values at the points `self.points(self.depth)`.
80
+    igral : float
81
+        The integral value of the interval.
82
+    err : float
83
+        The error associated with the integral value.
84
+    tol : float
85
+        The relative tolerance that needs to be reached in the precision of the integral.
86
+    rdepth : int
87
+        The number of splits that the interval has gone through, starting at 1.
88
+    ndiv : int
89
+        A number that is used to determine whether the interval is divergent.
90
+    parent : Interval
91
+        The parent interval. If the interval resulted from a refinement, it has one parent. If
92
+        it resulted from a split, it has two parents.
93
+    children : list of `Interval`s
94
+        The intervals resulting from a split or refinement.
95
+    done_points : dict
96
+        A dictionary with the x-values and y-values: `{x1: y1, x2: y2 ...}`.
97
+    est_err : float
98
+        The sum of the errors of the children, if one of the children is not ready yet,
99
+        the error is infinity.
100
+    discard : bool
101
+        If True, the interval and it's children are not participating in the
102
+        determination of the total integral anymore because its parent had a
103
+        refinement when the data of the interval was not known, and later it appears
104
+        that this interval has to be split.
105
+    complete : bool
106
+        All the function values in the interval are known. This does not necessarily mean
107
+        that the integral value has been calculated, see `self.done`.
108
+    done : bool
109
+        The integral and the error for the interval has been calculated.
110
+    branch_complete : bool
111
+        The interval can be used to determine the total integral, however if its children are
112
+        also `branch_complete`, they should be used.
113
+
114
+    """
115
+
116
+    __slots__ = ['a', 'b', 'c', 'c_old', 'depth', 'fx', 'igral', 'err', 'tol',
117
+                 'rdepth', 'ndiv', 'parent', 'children', 'done_points',
118
+                 'est_err', 'discard']
119
+
120
+    def __init__(self, a, b):
121
+        self.children = []
122
+        self.done_points = SortedDict()
123
+        self.a = a
124
+        self.b = b
125
+        self.c = np.zeros((len(ns), ns[-1]))
126
+        self.est_err = np.inf
127
+        self.discard = False
128
+        self.igral = None
129
+
130
+    @classmethod
131
+    def make_first(cls, a, b, tol):
132
+        ival = Interval(a, b)
133
+        ival.tol = tol
134
+        ival.ndiv = 0
135
+        ival.rdepth = 1
136
+        ival.parent = None
137
+        ival.depth = 3
138
+        ival.c_old = np.zeros(ns[ival.depth])
139
+        ival.err = np.inf
140
+        return ival, ival.points(ival.depth)
141
+
142
+    @property
143
+    def complete(self):
144
+        """The interval has all the y-values to calculate the intergral."""
145
+        return len(self.done_points) == ns[self.depth]
146
+
147
+    @property
148
+    def done(self):
149
+        """The interval is complete and has the intergral calculated."""
150
+        return hasattr(self, 'fx') and self.complete
151
+
152
+    @property
153
+    def branch_complete(self):
154
+        if not self.children and self.complete:
155
+            return True
156
+        else:
157
+            return np.isfinite(sum(i.est_err for i in self.children))
158
+
159
+    @property
160
+    def T(self):
161
+        """Get the correct shift matrix.
162
+
163
+        Should only be called on children of a split interval.
164
+        """
165
+        assert self.parent is not None
166
+        left = self.a == self.parent.a
167
+        right = self.b == self.parent.b
168
+        assert left != right
169
+        return T_left if left else T_right
170
+
171
+    def points(self, depth):
172
+        a = self.a
173
+        b = self.b
174
+        return (a+b)/2 + (b-a)*xi[depth]/2
175
+
176
+    def refine(self):
177
+        ival = Interval(self.a, self.b)
178
+        ival.tol = self.tol
179
+        ival.rdepth = self.rdepth
180
+        ival.ndiv = self.ndiv
181
+        ival.c = self.c.copy()
182
+        ival.c_old = self.c_old.copy()
183
+        ival.parent = self
184
+        self.children = [ival]
185
+        ival.err = self.err
186
+        ival.depth = self.depth + 1
187
+        points = ival.points(ival.depth)
188
+        return ival, points
189
+
190
+    def split(self):
191
+        points = self.points(self.depth)
192
+
193
+        a = self.a
194
+        b = self.b
195
+        m = points[len(points) // 2]
196
+
197
+        ivals = [Interval(a, m), Interval(m, b)]
198
+        self.children = ivals
199
+
200
+        for ival in ivals:
201
+            ival.depth = 0
202
+            ival.tol = self.tol / sqrt(2)
203
+            ival.c_old = self.c_old.copy()
204
+            ival.rdepth = self.rdepth + 1
205
+            ival.parent = self
206
+            ival.ndiv = self.ndiv
207
+            ival.err = self.err / sqrt(2)
208
+
209
+        return ivals
210
+
211
+    def complete_process(self):
212
+        """Calculate the integral contribution and error from this interval,
213
+        and update the estimated error of all ancestor intervals."""
214
+        force_split = False
215
+        if self.parent is None:
216
+            self.process_make_first()
217
+        elif self.rdepth > self.parent.rdepth:
218
+            self.process_split()
219
+        else:
220
+            force_split = self.process_refine()
221
+
222
+        # Set the estimated error
223
+        if np.isinf(self.est_err):
224
+            self.est_err = self.err
225
+        ival = self.parent
226
+        while ival is not None:
227
+            # update the error estimates on all ancestor intervals
228
+            children_err = sum(i.est_err for i in ival.children)
229
+            if np.isfinite(children_err):
230
+                ival.est_err = children_err
231
+                ival = ival.parent
232
+            else:
233
+                break
234
+
235
+        # Check whether the point spacing is smaller than machine precision
236
+        # and pop the interval with the largest error and do not split
237
+        remove = self.err < (abs(self.igral) * eps * Vcond[self.depth])
238
+        if remove:
239
+            # If this interval is discarded from ivals, there is no need
240
+            # to split it further.
241
+            force_split = False
242
+
243
+        return force_split, remove
244
+
245
+    def process_make_first(self):
246
+        fx = np.array(self.done_points.values())
247
+        nans = _zero_nans(fx)
248
+
249
+        self.c[3] = V_inv[3] @ fx
250
+        self.c[2, :ns[2]] = V_inv[2] @ fx[:ns[3]:2]
251
+        fx[nans] = np.nan
252
+        self.fx = fx
253
+
254
+        self.c_old = np.zeros(fx.shape)
255
+        c_diff = norm(self.c[self.depth] - self.c[2])
256
+
257
+        a, b = self.a, self.b
258
+        self.err = (b - a) * c_diff
259
+        self.igral = (b - a) * self.c[self.depth, 0] / sqrt(2)
260
+
261
+        if c_diff / norm(self.c[3]) > 0.1:
262
+            self.err = max(self.err, (b-a) * norm(self.c[3]))
263
+
264
+    def process_split(self, ndiv_max=20):
265
+        fx = np.array(self.done_points.values())
266
+        self.c[self.depth, :ns[self.depth]] = c_new = _calc_coeffs(fx, self.depth)
267
+        self.fx = fx
268
+
269
+        parent = self.parent
270
+        self.c_old = self.T @ parent.c[parent.depth]
271
+        c_diff = norm(self.c[self.depth] - self.c_old)
272
+
273
+        a, b = self.a, self.b
274
+        self.err = (b - a) * c_diff
275
+        self.igral = (b - a) * self.c[self.depth, 0] / sqrt(2)
276
+
277
+        self.ndiv = (parent.ndiv
278
+                     + (abs(parent.c[0, 0]) > 0
279
+                        and self.c[0, 0] / parent.c[0, 0] > 2))
280
+
281
+        if self.ndiv > ndiv_max and 2*self.ndiv > self.rdepth:
282
+            raise DivergentIntegralError(self)
283
+
284
+    def process_refine(self):
285
+        fx = np.array(self.done_points.values())
286
+        self.c[self.depth, :ns[self.depth]] = c_new = _calc_coeffs(fx, self.depth)
287
+        self.fx = fx
288
+
289
+        c_diff = norm(self.c[self.depth - 1] - self.c[self.depth])
290
+
291
+        a, b = self.a, self.b
292
+        self.err = (b - a) * c_diff
293
+        self.igral = (b - a) * c_new[0] / sqrt(2)
294
+        nc = norm(c_new)
295
+        force_split = nc > 0 and c_diff / nc > 0.1
296
+        return force_split
297
+
298
+    def __repr__(self):
299
+        lst = [
300
+            '(a, b)=({:.5f}, {:.5f})'.format(self.a, self.b),
301
+            'depth={}'.format(self.depth),
302
+            'rdepth={}'.format(self.rdepth),
303
+            'err={:.5E}'.format(self.err),
304
+            'igral={:.5E}'.format(self.igral if self.igral else 0),
305
+            'est_err={:.5E}'.format(self.est_err),
306
+            'discard={}'.format(self.discard),
307
+        ]
308
+        return ' '.join(lst)
309
+
310
+    def equal(self, other, *, verbose=False):
311
+        """Note: Implementing __eq__ breaks SortedContainers in some way."""
312
+        if not self.complete:
313
+            if verbose:
314
+                print('Interval {} is not complete.'.format(self))
315
+            return False
316
+
317
+        slots = set(self.__slots__).intersection(other.__slots__)
318
+        same_slots = []
319
+        for s in slots:
320
+            a = getattr(self, s)
321
+            b = getattr(other, s)
322
+            is_equal = np.allclose(a, b, rtol=0, atol=eps, equal_nan=True)
323
+            if verbose and not is_equal:
324
+                print('self.{} - other.{} = {}'.format(s, s, a - b))
325
+            same_slots.append(is_equal)
326
+
327
+        return all(same_slots)
328
+
329
+
330
+class IntegratorLearner(BaseLearner):
331
+
332
+    def __init__(self, function, bounds, tol):
333
+        """
334
+        Parameters
335
+        ----------
336
+        function : callable: X → Y
337
+            The function to learn.
338
+        bounds : pair of reals
339
+            The bounds of the interval on which to learn 'function'.
340
+        tol : float
341
+            Relative tolerance of the error to the integral, this means that the
342
+            learner is done when: `tol > err / abs(igral)`.
343
+
344
+        Attributes
345
+        ----------
346
+        complete_branches : list of intervals
347
+            The intervals that can be used in the determination of the integral.
348
+        nr_points : int
349
+            The total number of evaluated points.
350
+        igral : float
351
+            The integral value in `self.bounds`.
352
+        err : float
353
+            The absolute error associated with `self.igral`.
354
+
355
+        Methods
356
+        ------- 
357
+        done : bool
358
+            Returns whether the `tol` has been reached.
359
+        plot : hv.Scatter
360
+            Plots all the points that are evaluated.
361
+        """
362
+        self.function = function
363
+        self.bounds = bounds
364
+        self.tol = tol
365
+        self.priority_split = []
366
+        self.ivals = SortedSet([], key=attrgetter('err'))
367
+        self.done_points = {}
368
+        self.not_done_points = set()
369
+        self._stack = []
370
+        self._err_final = 0
371
+        self._igral_final = 0
372
+        self.x_mapping = defaultdict(lambda: SortedSet([], key=attrgetter('rdepth')))
373
+        ival, points = Interval.make_first(*self.bounds, self.tol)
374
+        self._update_ival(ival, points)
375
+        self.first_ival = ival
376
+        self._complete_branches = []
377
+
378
+    def add_point(self, point, value):
379
+        if point not in self.x_mapping:
380
+            raise ValueError("Point {} doesn't belong to any interval"
381
+                             .format(point))
382
+        self.done_points[point] = value
383
+        self.not_done_points.discard(point)
384
+
385
+        # Select the intervals that have this point
386
+        ivals = self.x_mapping[point]
387
+        for ival in ivals:
388
+            ival.done_points[point] = value
389
+            if ival.complete and not ival.done and not ival.discard:
390
+                in_ivals = ival in self.ivals
391
+                self.ivals.discard(ival)
392
+                force_split, remove = ival.complete_process()
393
+                if remove:
394
+                    self._err_final += ival.err
395
+                    self._igral_final += ival.igral
396
+                elif in_ivals:
397
+                    self.ivals.add(ival)
398
+
399
+                if force_split:
400
+                    # Make sure that at the next execution of _fill_stack(),
401
+                    # this ival will be split.
402
+                    self.priority_split.append(ival)
403
+
404
+    def _update_ival(self, ival, points):
405
+        assert not ival.discard
406
+        for x in points:
407
+            # Update the mappings
408
+            self.x_mapping[x].add(ival)
409
+            if x in self.done_points:
410
+                self.add_point(x, self.done_points[x])
411
+            elif x not in self.not_done_points:
412
+                self.not_done_points.add(x)
413
+                self._stack.append(x)
414
+
415
+        # Add the new interval to the err sorted set
416
+        self.ivals.add(ival)
417
+
418
+    def set_discard(self, ival):
419
+        def _discard(ival):
420
+            ival.discard = True
421
+            self.ivals.discard(ival)
422
+            for point in self._stack:
423
+                # XXX: is this check worth it?
424
+                if all(i.discard for i in self.x_mapping[point]):
425
+                    self._stack.remove(point)
426
+            for child in ival.children:
427
+                _discard(child)
428
+        _discard(ival)
429
+
430
+    def choose_points(self, n):
431
+        points, loss_improvements = self.pop_from_stack(n)
432
+        n_left = n - len(points)
433
+        while n_left > 0:
434
+            assert n_left >= 0
435
+            self._fill_stack()
436
+            new_points, new_loss_improvements = self.pop_from_stack(n_left)
437
+            points += new_points
438
+            loss_improvements += new_loss_improvements
439
+            n_left -= len(new_points)
440
+
441
+        return points, loss_improvements
442
+
443
+    def pop_from_stack(self, n):
444
+        points = self._stack[:n]
445
+        self._stack = self._stack[n:]
446
+        loss_improvements = [max(ival.err for ival in self.x_mapping[x])
447
+                             for x in points]
448
+        return points, loss_improvements
449
+
450
+    def remove_unfinished(self):
451
+        pass
452
+
453
+    def _fill_stack(self):
454
+        # XXX: to-do if all the ivals have err=inf, take the interval
455
+        # with the lowest rdepth and no children.
456
+        if self.priority_split:
457
+            ival = self.priority_split.pop()
458
+            force_split = True
459
+            if ival.children:
460
+                # If the interval already has children (which is the result of an
461
+                # earlier refinement when the data of the interval wasn't known
462
+                # yet,) then discard the children and propagate it down.
463
+                for child in ival.children:
464
+                    self.set_discard(child)
465
+        else:
466
+            ival = self.ivals[-1]
467
+            force_split = False
468
+            assert not ival.children
469
+
470
+        # Remove the interval from the err sorted set because we're going to
471
+        # split or refine this interval
472
+        self.ivals.discard(ival)
473
+
474
+        # If the interval points are smaller than machine precision, then
475
+        # don't continue with splitting or refining.
476
+        points = ival.points(ival.depth)
477
+        reached_machine_tol = points[1] <= points[0] or points[-1] <= points[-2]
478
+
479
+        if (not ival.discard) and (not reached_machine_tol):
480
+            if ival.depth == 3 or force_split:
481
+                # Always split when depth is maximal or if refining didn't help
482
+                ivals_new = ival.split()
483
+                for ival_new in ivals_new:
484
+                    points = ival_new.points(depth=0)
485
+                    self._update_ival(ival_new, points)
486
+            else:
487
+                # Refine
488
+                ival_new, points = ival.refine()
489
+                self._update_ival(ival_new, points)
490
+
491
+        # Remove the smallest element if number of intervals is larger than 1000
492
+        if len(self.ivals) > 1000:
493
+            self.ivals.pop(0)
494
+
495
+        return self._stack
496
+
497
+    @staticmethod
498
+    def deepest_complete_branches(ival):
499
+        """Finds the deepest complete set of intervals starting from `ival`."""
500
+        complete_branches = []
501
+        def _find_deepest(ival):
502
+            children_err = (sum(i.est_err for i in ival.children)
503
+                            if ival.children else np.inf)
504
+            if np.isfinite(ival.est_err) and np.isinf(children_err):
505
+                complete_branches.append(ival)
506
+            else:
507
+                for i in ival.children:
508
+                    _find_deepest(i)
509
+        _find_deepest(ival)
510
+        return complete_branches
511
+
512
+    @property
513
+    def complete_branches(self):
514
+        if not self.first_ival.done:
515
+            return []
516
+
517
+        if not self._complete_branches:
518
+            self._complete_branches.append(self.first_ival)
519
+
520
+        complete_branches = []
521
+        for ival in self._complete_branches:
522
+            if ival.discard:
523
+                complete_branches = self.deepest_complete_branches(self.first_ival)
524
+                break
525
+            if not ival.children:
526
+                # If the interval has no children, than is already is the deepest
527
+                # complete branch.
528
+                complete_branches.append(ival)
529
+            else:
530
+                complete_branches.extend(self.deepest_complete_branches(ival))
531
+        self._complete_branches = complete_branches
532
+        return self._complete_branches
533
+
534
+    @property
535
+    def nr_points(self):
536
+        return len(self.done_points)
537
+
538
+    @property
539
+    def igral(self):
540
+        return sum(i.igral for i in self.complete_branches)
541
+
542
+    @property
543
+    def err(self):
544
+        complete_branches = self.complete_branches
545
+        if not complete_branches:
546
+            return np.inf
547
+        else:
548
+            return sum(i.err for i in complete_branches)
549
+
550
+    def done(self):
551
+        err = self.err
552
+        igral = self.igral
553
+        return (err == 0
554
+                or err < abs(igral) * self.tol
555
+                or (self._err_final > abs(igral) * self.tol
556
+                    and err - self._err_final < abs(igral) * self.tol)
557
+                or not self.ivals)
558
+
559
+    def loss(self, real=True):
560
+        return abs(abs(self.igral) * self.tol - self.err)
561
+
562
+    def equal(self, other, *, verbose=False):
563
+        """Note: `other` is a list of ivals."""
564
+        if len(self.ivals) != len(other):
565
+            if verbose:
566
+                print('len(self.ivals)={} != len(other)={}'.format(
567
+                    len(self.ivals), len(other)))
568
+            return False
569
+
570
+        ivals = [sorted(i, key=attrgetter('a')) for i in [self.ivals, other]]
571
+        return all(ival.equal(other_ival, verbose=verbose)
572
+                   for ival, other_ival in zip(*ivals))
573
+
574
+    def plot(self):
575
+        return hv.Scatter(self.done_points)
0 576
new file mode 100644
... ...
@@ -0,0 +1,231 @@
1
+# -*- coding: utf-8 -*-
2
+from copy import deepcopy
3
+import heapq
4
+import itertools
5
+from math import sqrt
6
+
7
+import holoviews as hv
8
+import numpy as np
9
+import sortedcontainers
10
+
11
+from .base_learner import BaseLearner
12
+
13
+class Learner1D(BaseLearner):
14
+    """Learns and predicts a function 'f:ℝ → ℝ'.
15
+
16
+    Parameters
17
+    ----------
18
+    function : callable
19
+        The function to learn. Must take a single real parameter and
20
+        return a real number.
21
+    bounds : pair of reals
22
+        The bounds of the interval on which to learn 'function'.
23
+    """
24
+
25
+    def __init__(self, function, bounds):
26
+        self.function = function
27
+
28
+        # A dict storing the loss function for each interval x_n.
29
+        self.losses = {}
30
+        self.losses_combined = {}
31
+
32
+        self.data = sortedcontainers.SortedDict()
33
+        self.data_interp = {}
34
+
35
+        # A dict {x_n: [x_{n-1}, x_{n+1}]} for quick checking of local
36
+        # properties.
37
+        self.neighbors = sortedcontainers.SortedDict()
38
+        self.neighbors_combined = sortedcontainers.SortedDict()
39
+
40
+        # Bounding box [[minx, maxx], [miny, maxy]].
41
+        self._bbox = [list(bounds), [np.inf, -np.inf]]
42
+
43
+        # Data scale (maxx - minx), (maxy - miny)
44
+        self._scale = [bounds[1] - bounds[0], 0]
45
+        self._oldscale = deepcopy(self._scale)
46
+
47
+        self.bounds = list(bounds)
48
+
49
+    @property
50
+    def data_combined(self):
51
+        return {**self.data, **self.data_interp}
52
+
53
+    def interval_loss(self, x_left, x_right, data):
54
+        """Calculate loss in the interval x_left, x_right.
55
+
56
+        Currently returns the rescaled length of the interval. If one of the
57
+        y-values is missing, returns 0 (so the intervals with missing data are
58
+        never touched. This behavior should be improved later.
59
+        """
60
+        y_right, y_left = data[x_right], data[x_left]
61
+        if self._scale[1] == 0:
62
+            return sqrt(((x_right - x_left) / self._scale[0])**2)
63
+        else:
64
+            return sqrt(((x_right - x_left) / self._scale[0])**2 +
65
+                        ((y_right - y_left) / self._scale[1])**2)
66
+
67
+    def loss(self, real=True):
68
+        losses = self.losses if real else self.losses_combined
69
+        if len(losses) == 0:
70
+            return float('inf')
71
+        else:
72
+            return max(losses.values())
73
+
74
+    def update_losses(self, x, data, neighbors, losses):
75
+        x_lower, x_upper = neighbors[x]
76
+        if x_lower is not None:
77
+            losses[x_lower, x] = self.interval_loss(x_lower, x, data)
78
+        if x_upper is not None:
79
+            losses[x, x_upper] = self.interval_loss(x, x_upper, data)
80
+        try:
81
+            del losses[x_lower, x_upper]
82
+        except KeyError:
83
+            pass
84
+
85
+    def find_neighbors(self, x, neighbors):
86
+        pos = neighbors.bisect_left(x)
87
+        x_lower = neighbors.iloc[pos-1] if pos != 0 else None
88
+        x_upper = neighbors.iloc[pos] if pos != len(neighbors) else None
89
+        return x_lower, x_upper
90
+
91
+    def update_neighbors(self, x, neighbors):
92
+        if x not in neighbors:  # The point is new
93
+            x_lower, x_upper = self.find_neighbors(x, neighbors)
94
+            neighbors[x] = [x_lower, x_upper]
95
+            neighbors.get(x_lower, [None, None])[1] = x
96
+            neighbors.get(x_upper, [None, None])[0] = x
97
+
98
+    def update_scale(self, x, y):
99
+        self._bbox[0][0] = min(self._bbox[0][0], x)
100
+        self._bbox[0][1] = max(self._bbox[0][1], x)
101
+        if y is not None:
102
+            self._bbox[1][0] = min(self._bbox[1][0], y)
103
+            self._bbox[1][1] = max(self._bbox[1][1], y)
104
+
105
+        self._scale = [self._bbox[0][1] - self._bbox[0][0],
106
+                       self._bbox[1][1] - self._bbox[1][0]]
107
+
108
+    def add_point(self, x, y):
109
+        real = y is not None
110
+
111
+        if real:
112
+            # Add point to the real data dict and pop from the unfinished
113
+            # data_interp dict.
114
+            self.data[x] = y
115
+            try:
116
+                del self.data_interp[x]
117
+            except KeyError:
118
+                pass
119
+        else:
120
+            # The keys of data_interp are the unknown points
121
+            self.data_interp[x] = None
122
+
123
+        # Update the neighbors
124
+        self.update_neighbors(x, self.neighbors_combined)
125
+        if real:
126
+            self.update_neighbors(x, self.neighbors)
127
+
128
+        # Update the scale
129
+        self.update_scale(x, y)
130
+
131
+        # Interpolate
132
+        if not real:
133
+            self.data_interp = self.interpolate()
134
+
135
+        # Update the losses
136
+        self.update_losses(x, self.data_combined, self.neighbors_combined,
137
+                           self.losses_combined)
138
+        if real:
139
+            self.update_losses(x, self.data, self.neighbors, self.losses)
140
+
141
+        # If the scale has doubled, recompute all losses.
142
+        if self._scale > self._oldscale * 2:
143
+            self.losses = {xs: self.interval_loss(*xs, self.data)
144
+                           for xs in self.losses}
145
+            self.losses_combined = {x: self.interval_loss(*x,
146
+                                                          self.data_combined)
147
+                                    for x in self.losses_combined}
148
+            self._oldscale = self._scale
149
+
150
+    def choose_points(self, n, add_data=True):
151
+        """Return n points that are expected to maximally reduce the loss."""
152
+        # Find out how to divide the n points over the intervals
153
+        # by finding  positive integer n_i that minimize max(L_i / n_i) subject
154
+        # to a constraint that sum(n_i) = n + N, with N the total number of
155
+        # intervals.
156
+
157
+        # Return equally spaced points within each interval to which points
158
+        # will be added.
159
+        if n == 0:
160
+            return []
161
+
162
+        # If the bounds have not been chosen yet, we choose them first.
163
+        points = []
164
+        for bound in self.bounds:
165
+            if bound not in self.data and bound not in self.data_interp:
166
+                points.append(bound)
167
+
168
+        # Ensure we return exactly 'n' points.
169
+        if points:
170
+            loss_improvements = [float('inf')] * n
171
+            if n <= 2:
172
+                points = points[:n]
173
+            else:
174
+                points = np.linspace(*self.bounds, n)
175
+        else:
176
+            def xs(x, n):
177
+                if n == 1:
178
+                    return []
179
+                else:
180
+                    step = (x[1] - x[0]) / n
181
+                    return [x[0] + step * i for i in range(1, n)]
182
+
183
+            # Calculate how many points belong to each interval.
184
+            quals = [(-loss, x_range, 1) for (x_range, loss) in
185
+                     self.losses_combined.items()]
186
+
187
+            heapq.heapify(quals)
188
+
189
+            for point_number in range(n):
190
+                quality, x, n = quals[0]
191
+                heapq.heapreplace(quals, (quality * n / (n + 1), x, n + 1))
192
+
193
+            points = list(itertools.chain.from_iterable(xs(x, n)
194
+                          for quality, x, n in quals))
195
+
196
+            loss_improvements = list(itertools.chain.from_iterable(
197
+                                     itertools.repeat(-quality, n)
198
+                                     for quality, x, n in quals))
199
+
200
+        if add_data:
201
+            self.add_data(points, itertools.repeat(None))
202
+
203
+        return points, loss_improvements
204
+
205
+    def interpolate(self, extra_points=None):
206
+        xs = list(self.data.keys())
207
+        ys = list(self.data.values())
208
+        xs_unfinished = list(self.data_interp.keys())
209
+
210
+        if extra_points is not None:
211
+            xs_unfinished += extra_points
212
+
213
+        if len(ys) == 0:
214
+            interp_ys = (0,) * len(xs_unfinished)
215
+        else:
216
+            interp_ys = np.interp(xs_unfinished, xs, ys)
217
+
218
+        data_interp = {x: y for x, y in zip(xs_unfinished, interp_ys)}
219
+
220
+        return data_interp
221
+
222
+    def plot(self):
223
+            if self.data:
224
+                return hv.Scatter(self.data)
225
+            else:
226
+                return hv.Scatter([])
227
+
228
+    def remove_unfinished(self):
229
+        self.data_interp = {}
230
+        self.losses_combined = deepcopy(self.losses)
231
+        self.neighbors_combined = deepcopy(self.neighbors)
0 232
new file mode 100644
... ...
@@ -0,0 +1,303 @@
1
+# -*- coding: utf-8 -*-
2
+import itertools
3
+
4
+import holoviews as hv
5
+import numpy as np
6
+from scipy import interpolate, special
7
+
8
+from .base_learner import BaseLearner
9
+from .utils import restore
10
+
11
+
12
+# Learner2D and helper functions.
13
+
14
+def _losses_per_triangle(ip):
15
+    tri = ip.tri
16
+    vs = ip.values.ravel()
17
+
18
+    gradients = interpolate.interpnd.estimate_gradients_2d_global(
19
+        tri, vs, tol=1e-6)
20
+    p = tri.points[tri.vertices]
21
+    g = gradients[tri.vertices]
22
+    v = vs[tri.vertices]
23
+    n_points_per_triangle = p.shape[1]
24
+
25
+    dev = 0
26
+    for j in range(n_points_per_triangle):
27
+        vest = v[:, j, None] + ((p[:, :, :] - p[:, j, None, :]) *
28
+                                g[:, j, None, :]).sum(axis=-1)
29
+        dev += abs(vest - v).max(axis=1)
30
+
31
+    q = p[:, :-1, :] - p[:, -1, None, :]
32
+    areas = abs(q[:, 0, 0] * q[:, 1, 1] - q[:, 0, 1] * q[:, 1, 0])
33
+    areas /= special.gamma(n_points_per_triangle)
34
+    areas = np.sqrt(areas)
35
+
36
+    vs_scale = vs[tri.vertices].ptp()
37
+    if vs_scale != 0:
38
+        dev /= vs_scale
39
+
40
+    return dev * areas
41
+
42
+class Learner2D(BaseLearner):
43
+    """Learns and predicts a function 'f: ℝ^2 → ℝ'.
44
+
45
+    Parameters
46
+    ----------
47
+    function : callable
48
+        The function to learn. Must take a tuple of two real
49
+        parameters and return a real number.
50
+    bounds : list of 2-tuples
51
+        A list ``[(a1, b1), (a2, b2)]`` containing bounds,
52
+        one per dimension.
53
+
54
+    Attributes
55
+    ----------
56
+    points_combined
57
+        Sample points so far including the unknown interpolated ones.
58
+    values_combined
59
+        Sampled values so far including the unknown interpolated ones.
60
+    points
61
+        Sample points so far with real results.
62
+    values
63
+        Sampled values so far with real results.
64
+
65
+    Notes
66
+    -----
67
+    Adapted from an initial implementation by Pauli Virtanen.
68
+
69
+    The sample points are chosen by estimating the point where the
70
+    linear and cubic interpolants based on the existing points have
71
+    maximal disagreement. This point is then taken as the next point
72
+    to be sampled.
73
+
74
+    In practice, this sampling protocol results to sparser sampling of
75
+    smooth regions, and denser sampling of regions where the function
76
+    changes rapidly, which is useful if the function is expensive to
77
+    compute.
78
+
79
+    This sampling procedure is not extremely fast, so to benefit from
80
+    it, your function needs to be slow enough to compute.
81
+    """
82
+
83
+    def __init__(self, function, bounds):
84
+        self.ndim = len(bounds)
85
+        if self.ndim != 2:
86
+            raise ValueError("Only 2-D sampling supported.")
87
+        self.bounds = tuple((float(a), float(b)) for a, b in bounds)
88
+        self._points = np.zeros([100, self.ndim])
89
+        self._values = np.zeros([100], dtype=float)
90
+        self._stack = []
91
+        self._interp = {}
92
+
93
+        xy_mean = np.mean(self.bounds, axis=1)
94
+        xy_scale = np.ptp(self.bounds, axis=1)
95
+
96
+        def scale(points):
97
+            return (points - xy_mean) / xy_scale
98
+
99
+        def unscale(points):
100
+            return points * xy_scale + xy_mean
101
+
102
+        self.scale = scale
103
+        self.unscale = unscale
104
+
105
+        # Keeps track till which index _points and _values are filled
106
+        self.n = 0
107
+
108
+        self._bounds_points = list(itertools.product(*bounds))
109
+
110
+        # Add the loss improvement to the bounds in the stack
111
+        self._stack = [(*p, np.inf) for p in self._bounds_points]
112
+
113
+        self.function = function
114
+
115
+    @property
116
+    def points_combined(self):
117
+        return self._points[:self.n]
118
+
119
+    @property
120
+    def values_combined(self):
121
+        return self._values[:self.n]
122
+
123
+    @property
124
+    def points(self):
125
+        return np.delete(self.points_combined,
126
+                         list(self._interp.values()), axis=0)
127
+
128
+    @property
129
+    def values(self):
130
+        return np.delete(self.values_combined,
131
+                         list(self._interp.values()), axis=0)
132
+
133
+    def ip(self):
134
+        points = self.scale(self.points)
135
+        return interpolate.LinearNDInterpolator(points, self.values)
136
+
137
+    @property
138
+    def n_real(self):
139
+        return self.n - len(self._interp)
140
+
141
+    def ip_combined(self):
142
+        points = self.scale(self.points_combined)
143
+        values = self.values_combined
144
+
145
+        # Interpolate the unfinished points
146
+        if self._interp:
147
+            n_interp = list(self._interp.values())
148
+            bounds_are_done = not any(p in self._interp
149
+                                      for p in self._bounds_points)
150
+            if bounds_are_done:
151
+                values[n_interp] = self.ip()(points[n_interp])
152
+            else:
153
+                # It is important not to return exact zeros because
154
+                # otherwise the algo will try to add the same point
155
+                # to the stack each time.
156
+                values[n_interp] = np.random.rand(len(n_interp)) * 1e-15
157
+
158
+        return interpolate.LinearNDInterpolator(points, values)
159
+
160
+    def add_point(self, point, value):
161
+        nmax = self.values_combined.shape[0]
162
+        if self.n >= nmax:
163
+            self._values = np.resize(self._values, [2*nmax + 10])
164
+            self._points = np.resize(self._points, [2*nmax + 10, self.ndim])
165
+
166
+        point = tuple(point)
167
+
168
+        # When the point is not evaluated yet, add an entry to self._interp
169
+        # that saves the point and index.
170
+        if value is None:
171
+            self._interp[point] = self.n
172
+            old_point = False
173
+        else:
174
+            old_point = point in self._interp
175
+
176
+        # If the point is new add it a new value to _points and _values,
177
+        # otherwise get the index of the value that is being replaced.
178
+        if old_point:
179
+            n = self._interp.pop(point)
180
+        else:
181
+            n = self.n
182
+            self.n += 1
183
+
184
+        self._points[n] = point
185
+        self._values[n] = value
186
+
187
+        # Remove the point if in the stack.
188
+        for i, (*_point, _) in enumerate(self._stack):
189
+            if point == tuple(_point):
190
+                self._stack.pop(i)
191
+                break
192
+
193
+    def _fill_stack(self, stack_till=None):
194
+        if stack_till is None:
195
+            stack_till = 1
196
+
197
+        if self.values_combined.shape[0] < self.ndim + 1:
198
+            raise ValueError("too few points...")
199
+
200
+        # Interpolate
201
+        ip = self.ip_combined()
202
+        tri = ip.tri
203
+
204
+        losses = _losses_per_triangle(ip)
205
+
206
+        def point_exists(p):
207
+            eps = np.finfo(float).eps * self.points_combined.ptp() * 100
208
+            if abs(p - self.points_combined).sum(axis=1).min() < eps:
209
+                return True
210
+            if self._stack:
211
+                _stack_points, _ = self._split_stack()
212
+                if abs(p - np.asarray(_stack_points)).sum(axis=1).min() < eps:
213
+                    return True
214
+            return False
215
+
216
+        for j, _ in enumerate(losses):
217
+            # Estimate point of maximum curvature inside the simplex
218
+            jsimplex = np.argmax(losses)
219
+            p = tri.points[tri.vertices[jsimplex]]
220
+            point_new = self.unscale(p.mean(axis=-2))
221
+
222
+            # XXX: not sure whether this is necessary it was there
223
+            # originally.
224
+            point_new = np.clip(point_new, *zip(*self.bounds))
225
+
226
+            # Check if it is really new
227
+            if point_exists(point_new):
228
+                losses[jsimplex] = 0
229
+                continue
230
+
231
+            # Add to stack
232
+            self._stack.append((*point_new, losses[jsimplex]))
233
+
234
+            if len(self._stack) >= stack_till:
235
+                break
236
+            else:
237
+                losses[jsimplex] = 0
238
+
239
+    def _split_stack(self, n=None):
240
+        points = []
241
+        loss_improvements = []
242
+        for *point, loss_improvement in self._stack[:n]:
243
+            points.append(point)
244
+            loss_improvements.append(loss_improvement)
245
+        return points, loss_improvements
246
+
247
+    def _choose_and_add_points(self, n):
248
+        if n <= len(self._stack):
249
+            points, loss_improvements = self._split_stack(n)
250
+            self.add_data(points, itertools.repeat(None))
251
+        else:
252
+            points = []
253
+            loss_improvements = []
254
+            n_left = n
255
+            while n_left > 0:
256
+                # The while loop is needed because `stack_till` could be larger
257
+                # than the number of triangles between the points. Therefore
258
+                # it could fill up till a length smaller than `stack_till`.
259
+                if self.n >= 2**self.ndim:
260
+                    # Only fill the stack if no more bounds left in _stack
261
+                    self._fill_stack(stack_till=n_left)
262
+                new_points, new_loss_improvements = self._split_stack(n_left)
263
+                points += new_points
264
+                loss_improvements += new_loss_improvements
265
+                self.add_data(new_points, itertools.repeat(None))
266
+                n_left -= len(new_points)
267
+
268
+        return points, loss_improvements
269
+
270
+    def choose_points(self, n, add_data=True):
271
+        if not add_data:
272
+            with restore(self):
273
+                return self._choose_and_add_points(n)
274
+        else:
275
+            return self._choose_and_add_points(n)
276
+
277
+    def loss(self, real=True):
278
+        n = self.n_real if real else self.n
279
+        bounds_are_not_done = any(p in self._interp
280
+                                  for p in self._bounds_points)
281
+        if n <= 4 or bounds_are_not_done:
282
+            return np.inf
283
+        ip = self.ip() if real else self.ip_combined()
284
+        losses = _losses_per_triangle(ip)
285
+        return losses.max()
286
+
287
+    def remove_unfinished(self):
288
+        self._points = self.points.copy()
289
+        self._values = self.values.copy()
290
+        self.n -= len(self._interp)
291
+        self._interp = {}
292
+
293
+    def plot(self, n_x=201, n_y=201):
294
+        x, y = self.bounds
295
+        lbrt = x[0], y[0], x[1], y[1]
296
+        if self.n_real >= 4:
297
+            x = np.linspace(-0.5, 0.5, n_x)
298
+            y = np.linspace(-0.5, 0.5, n_y)
299
+            ip = self.ip()
300
+            z = ip(x[:, None], y[None, :])
301
+            return hv.Image(np.rot90(z), bounds=lbrt)
302
+        else:
303
+            return hv.Image(np.zeros((2, 2)), bounds=lbrt)
0 304
new file mode 100644
... ...
@@ -0,0 +1,11 @@
1
+# -*- coding: utf-8 -*-
2
+from contextlib import contextmanager
3
+
4
+@contextmanager
5
+def restore(*learners):
6
+    states = [learner.__getstate__() for learner in learners]
7
+    try:
8
+        yield
9
+    finally:
10
+        for state, learner in zip(states, learners):
11
+            learner.__setstate__(state)
0 12
new file mode 100644
... ...
@@ -0,0 +1,541 @@
1
+# Copyright 2010 Pedro Gonnet
2
+# Copyright 2017 Christoph Groth
3
+
4
+import warnings
5
+from fractions import Fraction as Frac
6
+from collections import defaultdict
7
+import numpy as np
8
+from numpy.testing import assert_allclose
9
+from numpy.linalg import cond
10
+from scipy.linalg import norm, inv
11
+
12
+
13
+eps = np.spacing(1)
14
+
15
+def legendre(n):
16
+    """Return the first n Legendre polynomials.
17
+
18
+    The polynomials have *standard* normalization, i.e.
19
+    int_{-1}^1 dx L_n(x) L_m(x) = delta(m, n) * 2 / (2 * n + 1).
20
+
21
+    The return value is a list of list of fraction.Fraction instances.
22
+    """
23
+    result = [[Frac(1)], [Frac(0), Frac(1)]]
24
+    if n <= 2:
25
+        return result[:n]
26
+    for i in range(2, n):
27
+        # Use Bonnet's recursion formula.
28
+        new = (i + 1) * [Frac(0)]
29
+        new[1:] = (r * (2*i - 1) for r in result[-1])
30
+        new[:-2] = (n - r * (i - 1) for n, r in zip(new[:-2], result[-2]))
31
+        new[:] = (n / i for n in new)
32
+        result.append(new)
33
+    return result
34
+
35
+
36
+def newton(n):
37
+    """Compute the monomial coefficients of the Newton polynomial over the
38
+    nodes of the n-point Clenshaw-Curtis quadrature rule.
39
+    """
40
+    # The nodes of the Clenshaw-Curtis rule are x_i = -cos(i * Pi / (n-1)).
41
+    # Here, we calculate the coefficients c_i such that sum_i c_i * x^i
42
+    # = prod_i (x - x_i).  The coefficients are thus sums of products of
43
+    # cosines.
44
+    #
45
+    # This routine uses the relation
46
+    #   cos(a) cos(b) = (cos(a + b) + cos(a - b)) / 2
47
+    # to efficiently calculate the coefficients.
48
+    #
49
+    # The dictionary 'terms' descibes the terms that make up the
50
+    # monomial coefficients.  Each item ((d, a), m) corresponds to a
51
+    # term m * cos(a * Pi / n) to be added to prefactor of the
52
+    # monomial x^(n-d).
53
+
54
+    mod = 2 * (n-1)
55
+    terms = defaultdict(int)
56
+    terms[0, 0] += 1
57
+
58
+    for i in range(n):
59
+        newterms = []
60
+        for (d, a), m in terms.items():
61
+            for b in [i, -i]:
62
+                # In order to reduce the number of terms, cosine
63
+                # arguments are mapped back to the inteval [0, pi/2).
64
+                arg = (a + b) % mod
65
+                if arg > n-1:
66
+                    arg = mod - arg
67
+                if arg >= n // 2:
68
+                    if n % 2 and arg == n // 2:
69
+                        # Zero term: ignore
70
+                        continue
71
+                    newterms.append((d + 1, n - 1 - arg, -m))
72
+                else:
73
+                    newterms.append((d + 1, arg, m))
74
+        for d, s, m in newterms:
75
+            terms[d, s] += m
76
+
77
+    c = (n + 1) * [0]
78
+    for (d, a), m in terms.items():
79
+        if m and a != 0:
80
+            raise ValueError("Newton polynomial cannot be represented exactly.")
81
+        c[n - d] += m
82
+        # The check could be removed and the above line replaced by
83
+        # the following, but then the result would be no longer exact.
84
+        # c[n - d] += m * np.cos(a * np.pi / (n - 1))
85
+
86
+    cf = np.array(c, float)
87
+    assert all(int(cfe) == ce for cfe, ce in zip(cf, c)), 'Precision loss'
88
+
89
+    cf /= 2.**np.arange(n, -1, -1)
90
+    return cf
91
+
92
+
93
+def scalar_product(a, b):
94
+    """Compute the polynomial scalar product int_-1^1 dx a(x) b(x).
95
+
96
+    The args must be sequences of polynomial coefficients.  This
97
+    function is careful to use the input data type for calculations.
98
+    """
99
+    la = len(a)
100
+    lc = len(b) + la + 1
101
+
102
+    # Compute the even coefficients of the product of a and b.
103
+    c = lc * [a[0].__class__()]
104
+    for i, bi in enumerate(b):
105
+        if bi == 0:
106
+            continue
107
+        for j in range(i % 2, la, 2):
108
+            c[i + j] += a[j] * bi
109
+
110
+    # Calculate the definite integral from -1 to 1.
111
+    return 2 * sum(c[i] / (i + 1) for i in range(0, lc, 2))
112
+
113
+
114
+def calc_bdef(ns):
115
+    """Calculate the decompositions of Newton polynomials (over the nodes
116
+    of the n-point Clenshaw-Curtis quadrature rule) in terms of
117
+    Legandre polynomials.
118
+
119
+    The parameter 'ns' is a sequence of numers of points of the
120
+    quadrature rule.  The return value is a corresponding sequence of
121
+    normalized Legendre polynomial coefficients.
122
+    """
123
+    legs = legendre(max(ns) + 1)
124
+    result = []
125
+    for n in ns:
126
+        poly = []
127
+        a = list(map(Frac, newton(n)))
128
+        for b in legs[:n + 1]:
129
+            igral = scalar_product(a, b)
130
+
131
+            # Normalize & store.  (The polynomials returned by
132
+            # legendre() have standard normalization that is not
133
+            # orthonormal.)
134
+            poly.append(np.sqrt((2*len(b) - 1) / 2) * igral)
135
+
136
+        result.append(np.array(poly))
137
+    return result
138
+
139
+
140
+# Nodes and Newton polynomials.
141
+n = (5, 9, 17, 33)
142
+xi = [-np.cos(np.arange(n[j])/(n[j]-1) * np.pi) for j in range(4)]
143
+# Make `xi` perfectly anti-symmetric, important for splitting the intervals
144
+xi = [(row - row[::-1]) / 2 for row in xi]
145
+b_def = calc_bdef(n)
146
+
147
+
148
+def calc_V(xi, n):
149
+    V = [np.ones(xi.shape), xi.copy()]
150
+    for i in range(2, n):
151
+        V.append((2*i-1) / i * xi * V[-1] - (i-1) / i * V[-2])
152
+    for i in range(n):
153
+        V[i] *= np.sqrt(i + 0.5)
154
+    return np.array(V).T
155
+
156
+# compute the coefficients
157
+V = [calc_V(*args) for args in zip(xi, n)]
158
+V_inv = list(map(inv, V))
159
+Vcond = list(map(cond, V))
160
+
161
+# shift matrix
162
+T_lr = [V_inv[3] @ calc_V((xi[3] + a) / 2, n[3]) for a in [-1, 1]]
163
+
164
+# compute the integral
165
+w = np.sqrt(0.5)                # legendre
166
+
167
+
168
+k = np.arange(n[3])
169
+alpha = np.sqrt((k+1)**2 / (2*k+1) / (2*k+3))
170
+gamma = np.concatenate([[0, 0], np.sqrt(k[2:]**2 / (4*k[2:]**2-1))])
171
+
172
+def _downdate(c, nans, depth):
173
+    b = b_def[depth].copy()
174
+    m = n[depth] - 1
175
+    for i in nans:
176
+        b[m + 1] /= alpha[m]
177
+        xii = xi[depth][i]
178
+        b[m] = (b[m] + xii * b[m + 1]) / alpha[m - 1]
179
+        for j in range(m - 1, 0, -1):
180
+            b[j] = ((b[j] + xii * b[j + 1] - gamma[j + 1] * b[j + 2])
181
+                    / alpha[j - 1])
182
+        b = b[1:]
183
+
184
+        c[:m] -= c[m] / b[m] * b[:m]
185
+        c[m] = 0
186
+        m -= 1
187
+
188
+
189
+def _zero_nans(fx):
190
+    nans = []
191
+    for i in range(len(fx)):
192
+        if not np.isfinite(fx[i]):
193
+            nans.append(i)
194
+            fx[i] = 0.0
195
+    return nans
196
+
197
+
198
+def _calc_coeffs(fx, depth):
199
+    """Caution: this function modifies fx."""
200
+    nans = _zero_nans(fx)
201
+    c_new = V_inv[depth] @ fx
202
+    if nans:
203
+        fx[nans] = np.nan
204
+        _downdate(c_new, nans, depth)
205
+    return c_new
206
+
207
+
208
+class DivergentIntegralError(ValueError):
209
+    def __init__(self, msg, igral, err, nr_points):
210
+        self.igral = igral
211
+        self.err = err
212
+        self.nr_points = nr_points
213
+        super().__init__(msg)
214
+
215
+
216
+class _Interval:
217
+    __slots__ = ['a', 'b', 'c', 'c_old', 'fx', 'igral', 'err', 'tol',
218
+                 'depth', 'rdepth', 'ndiv']
219
+
220
+    @classmethod
221
+    def make_first(cls, f, a, b, tol):
222
+        points = (a+b)/2 + (b-a) * xi[3] / 2
223
+        fx = f(points)
224
+        nans = _zero_nans(fx)
225
+        ival = _Interval()
226
+        ival.c = np.zeros((4, n[3]))
227
+        ival.c[3] = V_inv[3] @ fx
228
+        ival.c[2, :n[2]] = V_inv[2] @ fx[:n[3]:2]
229
+        fx[nans] = np.nan
230
+        ival.fx = fx
231
+        ival.c_old = np.zeros(fx.shape)
232
+        ival.a = a
233
+        ival.b = b
234
+        ival.igral = (b-a) * w * ival.c[3, 0]
235
+        c_diff = norm(ival.c[3] - ival.c[2])
236
+        ival.err = (b-a) * c_diff
237
+        if c_diff / norm(ival.c[3]) > 0.1:
238
+            ival.err = max( ival.err , (b-a) * norm(ival.c[3]) )
239
+        ival.tol = tol
240
+        ival.depth = 3
241
+        ival.ndiv = 0
242
+        ival.rdepth = 1
243
+        return ival, points
244
+
245
+    def split(self, f, ndiv_max=20):
246
+        a = self.a
247
+        b = self.b
248
+        m = (a + b) / 2
249
+        f_center = self.fx[(len(self.fx)-1)//2]
250
+
251
+        ivals = []
252
+        nr_points = 0
253
+        for aa, bb, f_left, f_right, T in [
254
+                (a, m, self.fx[0], f_center, T_lr[0]),
255
+                (m, b, f_center, self.fx[-1], T_lr[1])]:
256
+            ival = _Interval()
257
+            ivals.append(ival)
258
+            ival.a = aa
259
+            ival.b = bb
260
+            ival.tol = self.tol / np.sqrt(2)
261
+            ival.depth = 0
262
+            ival.rdepth = self.rdepth + 1
263
+            ival.c = np.zeros((4, n[3]))
264
+            fx = np.concatenate(
265
+                ([f_left],
266
+                 f((aa + bb) / 2 + (bb - aa) * xi[0][1:-1] / 2),
267
+                 [f_right]))
268
+            nr_points += n[0] - 2
269
+
270
+            ival.c[0, :n[0]] = c_new = _calc_coeffs(fx, 0)
271
+            ival.fx = fx
272
+
273
+            ival.c_old = T @ self.c[self.depth]
274
+            c_diff = norm(ival.c[0] - ival.c_old)
275
+            ival.err = (bb - aa) * c_diff
276
+            ival.igral = (bb - aa) * ival.c[0, 0] * w
277
+            ival.ndiv = (self.ndiv
278
+                         + (abs(self.c[0, 0]) > 0
279
+                            and ival.c[0, 0] / self.c[0, 0] > 2))
280
+            if ival.ndiv > ndiv_max and 2*ival.ndiv > ival.rdepth:
281
+                return (aa, bb, bb-aa), nr_points
282
+
283
+        return ivals, nr_points
284
+
285
+    def refine(self, f):
286
+        """Increase degree of interval."""
287
+        self.depth = depth = self.depth + 1
288
+        a = self.a
289
+        b = self.b
290
+        points = (a+b)/2 + (b-a)*xi[depth]/2
291
+        fx = np.empty(n[depth])
292
+        fx[0:n[depth]:2] = self.fx
293
+        fx[1:n[depth]-1:2] = f(points[1:n[depth]-1:2])
294
+        fx = fx[:n[depth]]
295
+        self.c[depth, :n[depth]] = c_new = _calc_coeffs(fx, depth)
296
+        self.fx = fx
297
+        c_diff = norm(self.c[depth - 1] - self.c[depth])
298
+        self.err = (b-a) * c_diff
299
+        self.igral = (b-a) * w * c_new[0]
300
+        nc = norm(c_new)
301
+        split = nc > 0 and c_diff / nc > 0.1
302
+
303
+        return points, split, n[depth] - n[depth - 1]
304
+
305
+
306
+def algorithm_4 (f, a, b, tol):
307
+    """ALGORITHM_4 evaluates an integral using adaptive quadrature. The
308
+    algorithm uses Clenshaw-Curtis quadrature rules of increasing
309
+    degree in each interval and bisects the interval if either the
310
+    function does not appear to be smooth or a rule of maximum degree
311
+    has been reached. The error estimate is computed from the L2-norm
312
+    of the difference between two successive interpolations of the
313
+    integrand over the nodes of the respective quadrature rules.
314
+
315
+    INT = ALGORITHM_4 ( F , A , B , TOL ) approximates the integral of
316
+    F in the interval [A,B] up to the relative tolerance TOL. The
317
+    integrand F should accept a vector argument and return a vector
318
+    result containing the integrand evaluated at each element of the
319
+    argument.
320
+
321
+    [INT,ERR,NR_POINTS] = ALGORITHM_4 ( F , A , B , TOL ) returns ERR,
322
+    an estimate of the absolute integration error as well as
323
+    NR_POINTS, the number of function values for which the integrand
324
+    was evaluated. The value of ERR may be larger than the requested
325
+    tolerance, indicating that the integration may have failed.
326
+
327
+    ALGORITHM_4 halts with a warning if the integral is or appears to
328
+    be divergent.
329
+
330
+    Reference: "Increasing the Reliability of Adaptive Quadrature
331
+        Using Explicit Interpolants", P. Gonnet, ACM Transactions on
332
+        Mathematical Software, 37 (3), art. no. 26, 2008.
333
+    """
334
+
335
+    # compute the first interval
336
+    ival, points = _Interval.make_first(f, a, b, tol)
337
+    ivals = [ival]
338
+
339
+    # init some globals
340
+    igral = ival.igral
341
+    err = ival.err
342
+    igral_final = 0
343
+    err_final = 0
344
+    i_max = 0
345
+    nr_points = n[3]
346
+
347
+    # do we even need to go this way?
348
+    if err < igral * tol:
349
+        return igral, err, nr_points
350
+
351
+    # main loop
352
+    for _ in range(int(1e9)):
353
+        if ivals[i_max].depth == 3:
354
+            split = True
355
+        else:
356
+            points, split, nr_points_inc = ivals[i_max].refine(f)
357
+            nr_points += nr_points_inc
358
+
359
+        # can we safely ignore this interval?
360
+        if (points[1] <= points[0]
361
+            or points[-1] <= points[-2]
362
+            or ivals[i_max].err < (abs(ivals[i_max].igral) * eps
363
+                                   * Vcond[ivals[i_max].depth])):
364
+            err_final += ivals[i_max].err
365
+            igral_final += ivals[i_max].igral
366
+            ivals[i_max] = ivals[-1]
367
+            ivals.pop()
368
+        elif split:
369
+            result, nr_points_inc = ivals[i_max].split(f)
370
+            nr_points += nr_points_inc
371
+            if isinstance(result, tuple):
372
+                igral = np.sign(igral) * np.inf
373
+                raise DivergentIntegralError(
374
+                    'Possibly divergent integral in the interval'
375
+                    ' [{}, {}]! (h={})'.format(*result),
376
+                    igral, err, nr_points)
377
+            ivals.extend(result)
378
+            ivals[i_max] = ivals.pop()
379
+
380
+        # compute the running err and new max
381
+        i_max = 0
382
+        i_min = 0
383
+        err = err_final
384
+        igral = igral_final
385
+        for i in range(len(ivals)):
386
+            if ivals[i].err > ivals[i_max].err:
387
+                i_max = i
388
+            elif ivals[i].err < ivals[i_min].err:
389
+                i_min = i
390
+            err += ivals[i].err
391
+            igral += ivals[i].igral
392
+
393
+        # nuke smallest element if stack is larger than 200
394
+        if len(ivals) > 200:
395
+            err_final += ivals[i_min].err
396
+            igral_final += ivals[i_min].igral
397
+            ivals[i_min] = ivals[-1]
398
+            ivals.pop()
399
+            if i_max == len(ivals):
400
+                i_max = i_min
401
+
402
+        # get up and leave?
403
+        if (err == 0
404
+            or err < abs(igral) * tol
405
+            or (err_final > abs(igral) * tol
406
+                and err - err_final < abs(igral) * tol)
407
+            or not ivals):
408
+            break
409
+
410
+    return igral, err, nr_points, ivals
411
+
412
+
413
+################ Tests ################
414
+
415
+def f0(x):
416
+    return x * np.sin(1/x) * np.sqrt(abs(1 - x))
417
+
418
+
419
+def f7(x):
420
+    return x**-0.5
421
+
422
+
423
+def f24(x):
424
+    return np.floor(np.exp(x))
425
+
426
+
427
+def f21(x):
428
+    y = 0
429
+    for i in range(1, 4):
430
+        y += 1 / np.cosh(20**i * (x - 2 * i / 10))
431
+    return y
432
+
433
+
434
+def f63(x):
435
+    return abs(x - 0.987654321)**-0.45
436
+
437
+
438
+def fdiv(x):
439
+    return abs(x - 0.987654321)**-1.1
440
+
441
+
442
+def test_legendre():
443
+    legs = legendre(11)
444
+    comparisons = [(legs[0], [1], 1),
445
+                    (legs[1], [0, 1], 1),
446
+                    (legs[10], [-63, 0, 3465, 0, -30030, 0,
447
+                                90090, 0, -109395, 0, 46189], 256)]
448
+    for a, b, div in comparisons:
449
+        for c, d in zip(a, b):
450
+            assert c * div == d
451
+
452
+
453
+def test_scalar_product(n=33):
454
+    legs = legendre(n)
455
+    selection = [0, 5, 7, n-1]
456
+    for i in selection:
457
+        for j in selection:
458
+            assert (scalar_product(legs[i], legs[j])
459
+                    == ((i == j) and Frac(2, 2*i + 1)))
460
+
461
+
462
+def simple_newton(n):
463
+    """Slower than 'newton()' and prone to numerical error."""
464
+    from itertools import combinations
465
+
466
+    nodes = -np.cos(np.arange(n) / (n-1) * np.pi)
467
+    return [sum(np.prod(-np.asarray(sel))
468
+                for sel in combinations(nodes, n - d))
469
+            for d in range(n + 1)]
470
+
471
+
472
+def test_newton():
473
+    assert_allclose(newton(9), simple_newton(9), atol=1e-15)
474
+
475
+
476
+def test_b_def(depth=1):
477
+    legs = [np.array(leg, float) for leg in legendre(n[depth] + 1)]
478
+    result = np.zeros(len(legs[-1]))
479
+    for factor, leg in zip(b_def[depth], legs):
480
+        factor *= np.sqrt((2*len(leg) - 1) / 2)
481
+        result[:len(leg)] += factor * leg
482
+    assert_allclose(result, newton(n[depth]), rtol=1e-15)
483
+
484
+
485
+def test_downdate(depth=3):
486
+    fx = np.abs(xi[depth])
487
+    fx[1::2] = np.nan
488
+    c_downdated = _calc_coeffs(fx, depth)
489
+
490
+    depth -= 1
491
+    fx = np.abs(xi[depth])
492
+    c = _calc_coeffs(fx, depth)
493
+
494
+    assert_allclose(c_downdated[:len(c)], c, rtol=0, atol=1e-9)
495
+
496
+
497
+def test_integration():
498
+    old_settings = np.seterr(all='ignore')
499
+
500
+    igral, err, nr_points = algorithm_4(f0, 0, 3, 1e-5)
501
+    assert_allclose(igral, 1.98194117954329, 1e-15)
502
+    assert_allclose(err, 1.9563545589988155e-05, 1e-10)
503
+    assert nr_points == 1129
504
+
505
+    igral, err, nr_points = algorithm_4(f7, 0, 1, 1e-6)
506
+    assert_allclose(igral, 1.9999998579359648, 1e-15)
507
+    assert_allclose(err, 1.8561437334964041e-06, 1e-10)
508
+    assert nr_points == 693
509
+
510
+    igral, err, nr_points = algorithm_4(f24, 0, 3, 1e-3)
511
+    assert_allclose(igral, 17.664696186312934, 1e-15)
512
+    assert_allclose(err, 0.017602618074957457, 1e-10)
513
+    assert nr_points == 4519
514
+
515
+    igral, err, nr_points = algorithm_4(f21, 0, 1, 1e-3)
516
+    assert_allclose(igral, 0.16310022131213361, 1e-15)
517
+    assert_allclose(err, 0.00011848806384952786, 1e-10)
518
+    assert nr_points == 191
519
+
520
+    igral, err, nr_points = algorithm_4(f63, 0, 1, 1e-10)
521
+    assert_allclose(igral, 1.967971650560763, 1e-15)
522
+    assert_allclose(err, 2.9049859499240667e-09, 1e-7)
523
+    assert nr_points == 2715
524
+
525
+    try:
526
+        igral, err, nr_points = algorithm_4(fdiv, 0, 1, 1e-6)
527
+    except DivergentIntegralError as e:
528
+        assert e.igral == np.inf
529
+        assert_allclose(e.err, 284.56192231467958, 1e-10)
530
+        assert e.nr_points == 431
531
+
532
+    np.seterr(**old_settings)
533
+
534
+
535
+if __name__ == '__main__':
536
+    test_legendre()
537
+    test_scalar_product()
538
+    test_newton()
539
+    test_b_def()
540
+    test_downdate()
541
+    test_integration()
... ...
@@ -175,11 +175,12 @@
175 175
    "metadata": {},
176 176
    "outputs": [],
177 177
    "source": [
178
-    "def func(xy):\n",
178
+    "def func(xy, wait=True):\n",
179 179
     "    import numpy as np\n",
180 180
     "    from time import sleep\n",
181 181
     "    from random import random\n",
182
-    "    sleep(random())\n",
182
+    "    if wait:\n",
183
+    "        sleep(random())\n",
183 184
     "    x, y = xy\n",
184 185
     "    a = 0.2\n",
185 186
     "    return x + np.exp(-(x**2 + y**2 - 0.75**2)**2/a**4)\n",
... ...
@@ -204,7 +205,7 @@
204 205
    "source": [
205 206
     "%%output size=100\n",
206 207
     "%%opts Contours (alpha=0.3)\n",
207
-    "from adaptive.learner import *\n",
208
+    "import holoviews as hv\n",
208 209
     "\n",
209 210
     "def plot(learner):\n",
210 211
     "    tri = learner.ip().tri\n",
... ...
@@ -230,7 +231,7 @@
230 231
     "learner2 = adaptive.learner.Learner2D(func, bounds=[(-1, 1), (-1, 1)])\n",
231 232
     "lin = np.linspace(-1, 1, len(learner.points)**0.5)\n",
232 233
     "xy = [(x, y) for x in lin for y in lin]\n",
233
-    "learner2.add_data(xy, map(func, xy))\n",
234
+    "learner2.add_data(xy, map(partial(func, wait=False), xy))\n",
234 235
     "learner2.plot().relabel('Homogeneous grid') + learner.plot().relabel('With adaptive')"
235 236
    ]
236 237
   },
... ...
@@ -275,11 +276,104 @@
275 276
    "metadata": {},
276 277
    "outputs": [],
277 278
    "source": [
278
-    "learner = adaptive.AverageLearner(g, None, 0.01)\n",
279
+    "learner = adaptive.learner.AverageLearner(g, None, 0.01)\n",
279 280
     "runner = adaptive.Runner(learner, goal=lambda l: l.loss() < 1)\n",
280 281
     "adaptive.live_plot(runner)"
281 282
    ]
282 283
   },
284
+  {
285
+   "cell_type": "markdown",
286
+   "metadata": {},
287
+   "source": [
288
+    "# 1D integration learner with `cquad`"
289
+   ]
290
+  },
291
+  {
292
+   "cell_type": "markdown",
293
+   "metadata": {},
294
+   "source": [
295
+    "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).\n",
296
+    "\n",
297
+    "Let's try the following function with cusps (that is difficult to integrate):"
298
+   ]
299
+  },
300
+  {
301
+   "cell_type": "code",
302
+   "execution_count": null,
303
+   "metadata": {},
304
+   "outputs": [],
305
+   "source": [
306
+    "import numpy as np\n",
307
+    "import holoviews as hv\n",
308
+    "\n",
309
+    "def f24(x):\n",
310
+    "    return np.floor(np.exp(x))\n",
311
+    "\n",
312
+    "xs = np.linspace(0, 3, 200)\n",
313
+    "hv.Scatter((xs, [f24(x) for x in xs]))"
314
+   ]
315
+  },
316
+  {
317
+   "cell_type": "markdown",
318
+   "metadata": {},
319
+   "source": [
320
+    "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."
321
+   ]
322
+  },
323
+  {
324
+   "cell_type": "code",
325
+   "execution_count": null,
326
+   "metadata": {},
327
+   "outputs": [],
328
+   "source": [
329
+    "import scipy.integrate\n",
330
+    "scipy.integrate.quad(f24, 0, 3)"
331
+   ]
332
+  },
333
+  {
334
+   "cell_type": "markdown",
335
+   "metadata": {},
336
+   "source": [
337
+    "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."
338
+   ]
339
+  },
340
+  {
341
+   "cell_type": "code",
342
+   "execution_count": null,
343
+   "metadata": {},
344
+   "outputs": [],
345
+   "source": [
346
+    "learner = adaptive.learner.IntegratorLearner(f24, bounds=(0, 3), tol=1e-10)\n",
347
+    "runner = adaptive.Runner(learner, executor=adaptive.runner.SequentialExecutor(), goal=lambda l: l.done())"
348
+   ]
349
+  },
350
+  {
351
+   "cell_type": "markdown",
352
+   "metadata": {},
353
+   "source": [
354
+    "Now we could do the live plotting again, but lets just wait untill the runner is done."
355
+   ]
356
+  },
357
+  {
358
+   "cell_type": "code",
359
+   "execution_count": null,
360
+   "metadata": {},
361
+   "outputs": [],
362
+   "source": [
363
+    "if not runner.task.done():\n",
364
+    "    raise RuntimeError('Wait for the runner to finish before executing the cells below!')"
365
+   ]
366
+  },
367
+  {
368
+   "cell_type": "code",
369
+   "execution_count": null,
370
+   "metadata": {},
371
+   "outputs": [],
372
+   "source": [
373
+    "print('The integral value is {} with the corresponding error of {}'.format(learner.igral, learner.err))\n",
374
+    "learner.plot()"
375
+   ]
376
+  },
283 377
   {
284 378
    "cell_type": "markdown",
285 379
    "metadata": {},
... ...
@@ -562,6 +656,61 @@
562 656
     "learner.plot().opts(style=dict(size=6)) * reconstructed_learner.plot()"
563 657
    ]
564 658
   },
659
+  {
660
+   "cell_type": "markdown",
661
+   "metadata": {},
662
+   "source": [
663
+    "### Timing functions"
664
+   ]
665
+  },
666
+  {
667
+   "cell_type": "markdown",
668
+   "metadata": {},
669
+   "source": [
670
+    "To time the runner you **cannot** simply use \n",
671
+    "```python\n",
672
+    "now = datetime.now()\n",
673
+    "runner = adaptive.Runner(...)\n",
674
+    "print(datetime.now() - now)\n",
675
+    "```\n",
676
+    "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.\n",
677
+    "\n",
678
+    "Therefore you need to create an `async` function and hook it into the `ioloop` like so:"
679
+   ]
680
+  },
681
+  {
682
+   "cell_type": "code",
683
+   "execution_count": null,
684
+   "metadata": {},
685
+   "outputs": [],
686
+   "source": [
687
+    "import asyncio\n",
688
+    "\n",
689
+    "async def time(runner):\n",
690
+    "    from datetime import datetime\n",
691
+    "    now = datetime.now()\n",
692
+    "    await runner.task\n",
693
+    "    return datetime.now() - now\n",
694
+    "\n",
695
+    "ioloop = asyncio.get_event_loop()\n",
696
+    "\n",
697
+    "learner = adaptive.learner.IntegratorLearner(f24, bounds=(0, 3), tol=1e-3)\n",
698
+    "runner = adaptive.Runner(learner, executor=adaptive.runner.SequentialExecutor(),\n",
699
+    "                         goal=lambda l: l.done())\n",
700
+    "\n",
701
+    "timer = ioloop.create_task(time(runner))"
702
+   ]
703
+  },
704
+  {
705
+   "cell_type": "code",
706
+   "execution_count": null,
707
+   "metadata": {},
708
+   "outputs": [],
709
+   "source": [
710
+    "# The result will only be set when the runner is done.\n",
711
+    "timer.result()"
712
+   ]
713
+  },
565 714
   {
566 715
    "cell_type": "markdown",
567 716
    "metadata": {},
... ...
@@ -1,3 +1,4 @@
1
+#!/usr/bin/env python3
1 2
 # -*- coding: utf-8 -*-
2 3
 
3 4
 from setuptools import setup
... ...
@@ -21,6 +22,7 @@ setup(
21 22
         'Topic :: Software Development :: Build Tools',
22 23
         'Programming Language :: Python :: 3.6',
23 24
     ],
24
-    packages=['adaptive'],
25
+    packages=['adaptive',
26
+              'adaptive.learner'],
25 27
     install_requires=requirements,
26 28
 )
27 29
new file mode 100644
... ...
@@ -0,0 +1,66 @@
1
+import numpy as np
2
+from adaptive.learner import IntegratorLearner
3
+from algorithm_4 import algorithm_4
4
+
5
+def same_ivals(f, a, b, tol, verbose):
6
+        igral, err, nr_points, ivals = algorithm_4(f, a, b, tol)
7
+
8
+        learner = IntegratorLearner(f, bounds=(a, b), tol=tol)
9
+        for i in range(nr_points):
10
+            points, loss_improvement = learner.choose_points(1)
11
+            learner.add_data(points, map(learner.function, points))
12
+        if verbose:
13
+            print('igral diff, ', learner.igral-igral, 'err diff', learner.err - err)
14
+        return learner.equal(ivals, verbose=verbose)
15
+
16
+
17
+def same_ivals_up_to(f, a, b, tol):
18
+        igral, err, nr_points, ivals = algorithm_4(f, a, b, tol)
19
+
20
+        learner = IntegratorLearner(f, bounds=(a, b), tol=tol)
21
+        j = 0
22
+        equal_till = 0
23
+        for i in range(nr_points):
24
+            points, loss_improvement = learner.choose_points(1)
25
+            learner.add_data(points, map(learner.function, points))
26
+            if not learner._stack:
27
+                try:
28
+                    j += 1
29
+                    if learner.equal(ivals):
30
+                        equal_till = i + 1
31
+                except:
32
+                    all_equal = False
33
+
34
+        return 'equal_till nr_points={} of {}'.format(equal_till, nr_points)
35
+
36
+if __name__ == '__main__':
37
+    old_settings = np.seterr(all='ignore')
38
+    from algorithm_4 import f0, f7, f24, f21, f63, fdiv
39
+    for i, args in enumerate([[f0, 0, 3, 1e-5],
40
+                              [f7, 0, 1, 1e-6],
41
+                              [f21, 0, 1, 1e-3],
42
+                              [f24, 0, 3, 1e-3],
43
+                              [f63, 0, 1, 1e-10]]):
44
+        print('\nFunction {}'.format(i))
45
+        if same_ivals(*args, verbose=True):
46
+            print(True)
47
+        else:
48
+            print(same_ivals_up_to(*args))
49
+    
50
+    # This function should raise a DivergentIntegralError.
51
+    print('Function ', i+1)
52
+    f, a, b, tol = [fdiv, 0, 1, 1e-6]
53
+    try:
54
+        igral, err, nr_points, ivals = algorithm_4(f, a, b, tol)
55
+    except Exception:
56
+        print('The integral is diverging.')
57
+
58
+    try:
59
+        learner = IntegratorLearner(f, bounds=(a, b), tol=tol)
60
+        for i in range(nr_points):
61
+            points, loss_improvement = learner.choose_points(1)
62
+            learner.add_data(points, map(learner.function, points))
63
+    except Exception:
64
+        print('The integral is diverging.')
65
+
66
+    np.seterr(**old_settings)