diff options
| author | Christoph Groth <christoph.groth@cea.fr> | 2018-02-21 16:51:41 +0100 |
|---|---|---|
| committer | Christoph Groth <christoph.groth@cea.fr> | 2018-02-23 12:55:00 +0100 |
| commit | 581bad7324c4c1164ba4adec69c7b11b0e25fb7e (patch) | |
| tree | eece21a8d282601ad08b2d3880d90be996e79afc /vquad | |
| parent | ee6289434905c0b3674e0f071ab0e27f77dc1417 (diff) | |
move interval quality check into refine()
Diffstat (limited to 'vquad')
| -rw-r--r-- | vquad/core.py | 31 |
1 files changed, 18 insertions, 13 deletions
diff --git a/vquad/core.py b/vquad/core.py index 7435334..f6baedb 100644 --- a/vquad/core.py +++ b/vquad/core.py @@ -159,14 +159,26 @@ class Vquad: return children def refine(self, ival): - """Increase degree of interval.""" + """Increase degree of interval. + + Returns True if the refined interval is OK, and False if it is + borderline and should not be refined/split further. This + happens when neigboring points can be only barely resolved by + floating point numbers, or when the estimated relative error is + already at the limit of numerical accuracy and cannot be reduced + further. + """ ival.level += 1 points = ival.points() vals = np.empty(points.shape) vals[0::2] = ival.vals vals[1::2] = self.f(points[1::2]) ival.interpolate(vals, ival.coeffs) - return points + + return (points[1] - points[0] > points[0] * min_sep + and points[-1] - points[-2] > points[-2] * min_sep + and ival.err > (abs(ival.igral) + * eps * tbls.V_cond_nums[ival.level])) def improve(self): ival = self.ivals[-1] @@ -174,21 +186,14 @@ class Vquad: if ival.level == max_level: split = True else: - points = self.refine(ival) - split = ival.unreliable_err - - if (points[1] - points[0] < points[0] * min_sep - or points[-1] - points[-2] < points[-2] * min_sep - or ival.err < (abs(ival.igral) - * eps * tbls.V_cond_nums[ival.level])): - # Remove the interval (while remembering the excess integral - # and error), since it is either too narrow, or the estimated - # relative error is already at the limit of numerical accuracy - # and cannot be reduced further. + if not self.refine(ival): + # Remove the interval but remember the excess integral and + # error. self.err_excess += ival.err self.igral_excess += ival.igral self.attic.append(self.ivals.pop()) return + split = ival.unreliable_err if split: # Replace current interval by its children. |
