diff options
| author | Christoph Groth <christoph.groth@cea.fr> | 2018-01-12 11:29:54 +0100 |
|---|---|---|
| committer | Christoph Groth <christoph.groth@cea.fr> | 2018-02-23 12:55:00 +0100 |
| commit | 781c2a57f4803cbc0cbc0368f9e353fc4d90ca9a (patch) | |
| tree | 58aa61a59c86b86347d4d8c6caac03e0a0442d0e /vquad | |
| parent | b09d4b8143d7eda458cd3cccf4f0af54e4d66272 (diff) | |
sort ivals instead of using i_max
Diffstat (limited to 'vquad')
| -rw-r--r-- | vquad/core.py | 40 |
1 files changed, 17 insertions, 23 deletions
diff --git a/vquad/core.py b/vquad/core.py index d7fcec4..7f586f9 100644 --- a/vquad/core.py +++ b/vquad/core.py @@ -98,6 +98,9 @@ class _Interval: self.err = w * coeffs_diff self.unreliable_err = coeffs_diff > hint * norm(coeffs) + def __lt__(self, other): + return self.err < other.err + class Vquad: """Evaluate an integral using adaptive quadrature. @@ -129,7 +132,6 @@ class Vquad: self.nr_points = len(vals) self.igral_excess = 0 self.err_excess = 0 - self.i_max = 0 def split(self, ival): m = (ival.a + ival.b) / 2 @@ -169,45 +171,37 @@ class Vquad: return points def improve(self): - i_max = self.i_max + ival = self.ivals[-1] - if self.ivals[i_max].level == max_level: + if ival.level == max_level: split = True else: - points = self.refine(self.ivals[i_max]) - split = self.ivals[i_max].unreliable_err + 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 (self.ivals[i_max].err - < (abs(self.ivals[i_max].igral) * eps - * tbls.V_cond_nums[self.ivals[i_max].level]))): + 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. - self.err_excess += self.ivals[i_max].err - self.igral_excess += self.ivals[i_max].igral - self.ivals[i_max] = self.ivals[-1] + self.err_excess += ival.err + self.igral_excess += ival.igral self.attic.append(self.ivals.pop()) return if split: - self.ivals.extend(self.split(self.ivals[i_max])) - self.ivals[i_max] = self.ivals.pop() + self.ivals.extend(self.split(self.ivals.pop())) + self.ivals.sort() def totals(self): - # Compute the total error and new max. - i_max = 0 - err = self.err_excess igral = self.igral_excess - for i in range(len(self.ivals)): - if self.ivals[i].err > self.ivals[i_max].err: - i_max = i - err += self.ivals[i].err - igral += self.ivals[i].igral - - self.i_max = i_max + err = self.err_excess + for ival in self.ivals: + igral += ival.igral + err += ival.err return igral, err def improve_until(self, tol): |