#!/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()