# Mathematical utility routines # Copyright (C) 1999, Wesley Phoa # # Reference: Numerical Recipes in C # Although lacking the proper copyright notice, the author of this module # states on his site (http://www.margaretmorgan.com/wesley/python/) # that the code is released under GNU GPL from operator import * from Numeric import * class BracketingException(Exception): pass class RootFindingException(Exception): pass class MinimizationException(Exception): pass GOLDEN = (1+5**.5)/2 LITTLE = 1e-10 SQPREC = 1e-4 # # MISCELLANEOUS # def sgn(x): if x==0: return 0 else: return x/abs(x) def along(f, x, v): """\ Given a multivariate function f, a point x and a vector v, return the univariate function t |-> f(x+tv). """ return lambda t,f=f,x=x,v=v: apply(f, add(x, multiply(t, v))) # # UNIVARIATE ROOT FINDING # def bracket_root(f, interval, max_iterations=50): """\ Given a univariate function f and a tuple interval=(x1,x2), return a new tuple (bracket, fnvals) where bracket=(x1,x2) brackets a root of f and fnvals=(f(x1),f(x2)). """ x1, x2 = interval if x1==x2: raise BracketingException("initial interval has zero width") elif x2= 0: # not currently bracketed if abs(f1)=0: raise BracketingException("initial interval does not bracket a root") x4 = 123456789. for j in range(max_iterations): x3 = (x1+x2)/2 f3 = f(x3) temp = f3*f3 - f1*f2 x4, x4old = x3 + (x3-x1)*sgn(f1-f2)*f3/temp**.5, x4 f4 = f(x4) if f1*f4<0: # x1 and x4 bracket root x2, f2 = x4, f4 else: # x4 and x2 bracket root x1, f1 = x4, f4 if min(abs(x1-x2),abs(x4-x4old))f1: # ensure x1 --> x2 is downhill direction x1, x2, f1, f2 = x2, x1, f2, f1 x3 = x2 + GOLDEN*(x2-x1) f3 = f(x3) for j in range(max_iterations): if f2x3: # ensure x1f2accuracy: r = (x-w)*(fx-fv) q = (x-v)*(fx-fw) p = (x-v)*q - (x-w)*r q = 2*(q-r) if q>0: p = -p q = abs(q) etemp = e e = d if abs(p)>=abs(q*etemp)/2 or p<=q*(a-x) or p>=q*(b-x): if x>=xm: e = a-x else: e = b-x d = (2-GOLDEN)*e else: # accept parabolic fit d = p/q u = x+d if u-a<2*accuracy or b-u<2*accuracy: d = accuracy*sgn(xm-x) else: if x>=xm: e = a-x else: e = b-x d = (2-GOLDEN)*e if abs(d)>=accuracy: u = x+d else: u = x+accuracy*sgn(d) fu = f(u) if fu<=fx: if u>=x: a = x else: b = x v, w, x = w, x, u fv, fw, fx = fw, fx, fu else: if umaxdrop: maxdrop_i = i maxdrop = fpold-fp totaldrop = fpold-fp if 2*totaldrop<=tolerance*(abs(fp)+abs(fpold)): if return_fnval: return tuple(p), fp else: return tuple(p) vnew = subtract(p, pold) pex = add(pold, multiply(2, vnew)) fpex = apply(f, pex) if fpex