/usr/share/pyshared/pytools/convergence.py is in python-pytools 2011.5-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 | import numpy as np
# {{{ 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 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 write_gnuplot_file(self, filename):
outfile = file(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)))
# }}}
# vim: foldmethod=marker
|