/usr/lib/python3/dist-packages/pytools/convergence.py is in python3-pytools 2015.1.6-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 | from __future__ import absolute_import
import numpy as np
from six.moves import range
from six.moves import zip
# {{{ eoc estimation --------------------------------------------------------------
def estimate_order_of_convergence(abscissae, errors):
"""Assuming that abscissae and errors are connected by a law of the form
error = constant * abscissa ^ (order),
this function finds, in a least-squares sense, the best approximation of
constant and order for the given data set. It returns a tuple (constant, order).
"""
assert len(abscissae) == len(errors)
if len(abscissae) <= 1:
raise RuntimeError("Need more than one value to guess order of convergence.")
coefficients = np.polyfit(np.log10(abscissae), np.log10(errors), 1)
return 10**coefficients[-1], coefficients[-2]
class EOCRecorder(object):
def __init__(self):
self.history = []
def add_data_point(self, abscissa, error):
self.history.append((abscissa, error))
def estimate_order_of_convergence(self, gliding_mean=None):
abscissae = np.array([a for a, e in self.history])
errors = np.array([e for a, e in self.history])
size = len(abscissae)
if gliding_mean is None:
gliding_mean = size
data_points = size - gliding_mean + 1
result = np.zeros((data_points, 2), float)
for i in range(data_points):
result[i, 0], result[i, 1] = estimate_order_of_convergence(
abscissae[i:i+gliding_mean], errors[i:i+gliding_mean])
return result
def order_estimate(self):
return self.estimate_order_of_convergence()[0, 1]
def max_error(self):
return max(err for absc, err in self.history)
def pretty_print(self, abscissa_label="h", error_label="Error", gliding_mean=2):
from pytools import Table
tbl = Table()
tbl.add_row((abscissa_label, error_label, "Running EOC"))
gm_eoc = self.estimate_order_of_convergence(gliding_mean)
for i, (absc, err) in enumerate(self.history):
if i < gliding_mean-1:
tbl.add_row((str(absc), str(err), ""))
else:
tbl.add_row((str(absc), str(err), str(gm_eoc[i-gliding_mean+1, 1])))
if len(self.history) > 1:
return str(tbl) + "\n\nOverall EOC: %s" \
% self.estimate_order_of_convergence()[0, 1]
else:
return str(tbl)
def __str__(self):
return self.pretty_print()
def write_gnuplot_file(self, filename):
outfile = open(filename, "w")
for absc, err in self.history:
outfile.write("%f %f\n" % (absc, err))
result = self.estimate_order_of_convergence()
const = result[0, 0]
order = result[0, 1]
outfile.write("\n")
for absc, err in self.history:
outfile.write("%f %f\n" % (absc, const * absc**(-order)))
# }}}
# {{{ p convergence verifier
class PConvergenceVerifier(object):
def __init__(self):
self.orders = []
self.errors = []
def add_data_point(self, order, error):
self.orders.append(order)
self.errors.append(error)
def __str__(self):
from pytools import Table
tbl = Table()
tbl.add_row(("p", "error"))
for p, err in zip(self.orders, self.errors):
tbl.add_row((str(p), str(err)))
return str(tbl)
def __call__(self):
orders = np.array(self.orders, np.float64)
errors = np.abs(np.array(self.errors, np.float64))
rel_change = np.diff(1e-20 + np.log10(errors)) / np.diff(orders)
assert (rel_change < -0.2).all()
# }}}
# vim: foldmethod=marker
|