Source code for scitools.convergencerate

#!/usr/bin/env python
"""
Module for estimating convergence rate of numerical algorithms,
based on data from a set of experiments.

Recommended usage:
Vary only discretization parameter (h in spatial problems, h/dt**q for
some q so that h/dt**q=const, in time-dependent problems with time step dt).
Use class OneDiscretizationPrm, or the function convergence_rate,
or the analyze_filedata convenience function. Start with reading
the convergence_rate function (too see an easily adapted example).

(There is support for fitting more general error models, like
C1*h**r1 + C2*h*dt**r2, with nonlinear least squares, but sound
fits are more challenging to obtain.)
"""

from Scientific.Functions.LeastSquares import leastSquaresFit
from numpy import zeros, array, asarray, log10, transpose, linalg, linspace
import sys
from scitools.std import plot

__all__ = ['ManyDiscretizationPrm',
           'OneDiscretizationPrm',
           ]

log = log10
inv_log = lambda x: 10**x

# The classes in this module have only static methods, i.e.,
# classes are merely name spaces for related functions.


[docs]class OneDiscretizationPrm(object): """ Tools for fitting an error model with only one discretization parameter: e = C*h^2. """
[docs] def error_model(p, d): """Return e = C*h**a, where p=[C, a] and h=d[0].""" C, a = p h = d[0] return C*h**a
[docs] def loglog_error_model(p, d): """As error_model if log-log data was used in the estimation.""" C, a = p h = d[0] return log(C) + a*log(h)
[docs] def linear_loglog_fit(d, e): """ Linear least squares algorithm. Suitable for problems with only one distinct discretization parameter. *d* is the sequence of discretization parameter values, and *e* is the sequence of corresponding error values. The error model the data is supposed to fit reads ``log(e[i]) = log(C[i]) + a*log(d[i])``. """ A = transpose(array([d, zeros(len(d))+1])) sol = linalg.lstsq(A, e) a, logC = sol[0] C = inv_log(logC) return a, C
[docs] def nonlinear_fit(d, e, p0): """ ======== =========================================================== Name Description ======== =========================================================== d list of values of the (single) discretization parameter in each experiment: d[i] provides the values of the discretization, parameter in experiement no. i. e list of error values; e = (e_1, e_2, ...), e[i] is the error associated with the parameters d[i] p0 starting values for the unknown parameters vector return a, C; a is the exponent, C is the factor in front. ======== =========================================================== """ if len(d) != len(e): raise ValueError('d and e must have the same length') if not isinstance(d[0], (float,int)): raise TypeError('d must be an array of numbers, not %s' % \ str(type(d[0]))) # transform d and e to the data format required by # the Scientific package: data = [] for d_i, e_i in zip(d, e): data.append(((d_i,) , e_i)) # recall (a,) conversion to tuple sol = leastSquaresFit(OneDiscretizationPrm.error_model, p0, data) C = sol[0][0] a = sol[0][1] return a, C
[docs] def pairwise_rates(d, e): """ Compare convergence rates, where each rate is based on a formula for two successive experiments. """ if len(d) != len(e): raise ValueError('d and e must have the same length') rates = [] for i in range(1, len(d)): try: rates.append( log(e[i-1]/e[i])/log(float(d[i-1])/d[i]) ) except (ZeroDivisionError, OverflowError): rates.append(0) # estimate C from the last data point: try: C = e[-1]/d[-1]**rates[-1] except: C = 0 return rates, C
[docs] def analyze(d, e, initial_guess=None, plot_title='', filename='tmp.ps'): """ Run linear, nonlinear and successive rates models. d: list/array of discretization parameter. e: errors corresponding to d. Plot results for comparison of the three approaches. """ # convert to NumPy arrays: d = asarray(d, float); e = asarray(e, float) # linear least squares fit: a1, C1 = OneDiscretizationPrm.linear_loglog_fit(log(d), log(e)) print "linear LS fit: const=%e rate=%.1f" % (C1, a1) # nonlinear least squares fit (no log-log): a2, C2 = OneDiscretizationPrm.nonlinear_fit(d, e, initial_guess) print "nonlinear LS fit: const=%e rate=%.1f" % (C2, a2) # pairwise estimate of the rate: rates, C3 = OneDiscretizationPrm.pairwise_rates(d, e) a3 = rates[-1] print "pairwise fit: const=%e rate=%.1f" % (C3, a3) print "all rates:", rates if C1 < 0 or C2 < 0 or C3 < 0: raise ValueError('Some fits give negative const value! Cannot plot.') return # else log plot: log_d1 = linspace(log(d[0]), log(d[-1]), 2) log_e1 = log(C1) + a1*log_d1 log_e2 = log(C2) + a2*log_d1 log_e3 = log(C3) + a3*log_d1 plot(log(d), log(e), 'yo', log_d1, log_e1, 'r-', log_d1, log_e2, 'b-', log_d1, log_e3, 'g-', legend=('data', 'linear LS log-log fit: %.1f*h^%.1f' % (log(C1), a1), 'nonlinear LS fit: %.1f*h^%.1f' % (log(C2), a2), 'successive rates, last two experiments: %.1f*h^%.1f' % (log(C3), a3)), #axis=[log_d1[-1], log_d1[0], 1.3*log(e[-1]), 1.5*log(e[0])], hardcopy=filename)
analyze = staticmethod(analyze) error_model = staticmethod(error_model) loglog_error_model = staticmethod(loglog_error_model) linear_loglog_fit = staticmethod(linear_loglog_fit) nonlinear_fit = staticmethod(nonlinear_fit) pairwise_rates = staticmethod(pairwise_rates) # convenience function:
def convergence_rate(discretization_prm, error): """ Given two lists/arrays with discretization parameters and corresponding errors in a numerical method (element no. i in the two lists must correspond to each other), this function assumes an error formula of the form E=C*d^r, where E is the error and d is the discretization parameter. The function returns C and r. Method used: OneDiscretizationPrm.pairwise_rates is called and the final r value is used for return. A check that the rates are converging (the last three) is done. """ rates, C = OneDiscretizationPrm.pairwise_rates(discretization_prm, error) # check that there is no divergence at the end of # the series of experiments differences = [rates[i] - rates[i-1] for i in range(len(rates)-1)] # the differences between the rates should decrease, at least # toward the end try: if differences[-1] > differences[-2]: raise ValueError('The pairwise convergence rates do not '\ 'decrease toward the end:\n%s' % \ str(rates)) except IndexError: pass # not enough data to check the differences list return C, rates[-1]
[docs]class ManyDiscretizationPrm(object): """ General tool for fitting an error model containing an arbitrary number of discretization parameters. The error is a weighted sum of each discretization parameter raised to a real expoenent. The weights and exponents are the unknown parameters to be fitted by a nonlinear least squares procedure. """
[docs] def error_model(p, d): """ Evaluate the theoretical error model (sum of C*h^r terms): sum_i p[2*i]*d[i]**p[2*i+1] ====== ======================================================= Name Description ====== ======================================================= p sequence ofvalues of parameters (estimated) d sequence of values of (known) discretization parameters return error evaluated ====== ======================================================= Note that ``len(p)`` must be ``2*len(d)`` in this model since there are two parameters (constant and exponent) for each discretization parameter. """ if len(p) != 2*len(d): raise ValueError('len(p)=%d != 2*len(d)=%d' % (len(p),2*len(d))) sum = 0 for i in range(len(d)): sum += p[2*i] * d[i]**p[2*i+1] return sum
[docs] def nonlinear_fit(d, e, initial_guess): """ @param d: list of values of the set of discretization parameters in each experiment: d = ((d_1,d_2,d_3),(d_1,d_2,d_3,),...); d[i] provides the values of the discretization parameters in experiement no. i. @param e: list of error values; e = (e_1, e_2, ...): e[i] is the error associated with the parameters d[i]. @param initial_guess: the starting value for the unknown parameters vector. @return: list of fitted paramters. """ if len(d) != len(e): raise ValueError('len(d) != len(e)') # transform d and e to the data format required by # the Scientific package: data = [] for d_i, e_i in zip(d, e): if isinstance(d_i, (float, int)): data.append(((d_i,), e_i)) else: # d_i is tuple, list, array, NumArray, ... data.append((d_i, e_i)) sol = leastSquaresFit(ManyDiscretizationPrm.error_model, initial_guess, data) # return list of fitted parameters (p in error_model) # (sol[1] is a measure of the quality of the fit) return sol[0]
error_model = staticmethod(error_model) nonlinear_fit = staticmethod(nonlinear_fit)
def _test1(): """Single discretization parameter test.""" import random random.seed(1234) n = 7 h = 1 e = []; d = [] for i in range(7): h /= 2.0 error = OneDiscretizationPrm.error_model((4,2), (h,)) error += random.gauss(0, 0.1*error) # perturb data d.append(h) e.append(error) OneDiscretizationPrm.analyze(d, e, initial_guess=(3,2)) def analyze_filedata(): # read filename and initial guess of C and r in error formula E=C*h^r f = open(sys.argv[1], 'r') C = float(sys.argv[2]) r = float(sys.argv[3]) try: plot_title = sys.argv[4] except: plot_title = '' from scitools.filetable import read data = read(f) print data OneDiscretizationPrm.analyze(data[:,0], data[:,1], initial_guess=(C,r), plot_title=plot_title) # extensions: # example with dx, dy and dt # same example, but with factors to get a common rate # dx, dt tables and experiments with whole table, one # column and one row, and the diagonal if __name__ == '__main__': if len(sys.argv) == 1: print 'Usage: %s filename C-guess r-guess [plot-title]' % sys.argv[0] elif sys.argv[1] == 'example': _test1() else: analyze_filedata()