source:orange/orange/mathutil.py@6538:a5f65d7f0b2c

Revision 6538:a5f65d7f0b2c, 6.6 KB checked in by Mitar <Mitar@…>, 4 years ago (diff)

Made XPM version of the icon 32x32.

Line
1#   Mathematical utility routines
2#   Copyright (C) 1999, Wesley Phoa
3#
4#   Reference: Numerical Recipes in C
5
6# Although lacking the proper copyright notice, the author of this module
7# states on his site (http://www.margaretmorgan.com/wesley/python/)
8# that the code is released under GNU GPL
9
10from operator import *
11from Numeric import *
12
13class BracketingException(Exception):
14    pass
15
16class RootFindingException(Exception):
17    pass
18
19class MinimizationException(Exception):
20    pass
21
22GOLDEN = (1+5**.5)/2
23LITTLE = 1e-10
24SQPREC = 1e-4
25
26#
27# MISCELLANEOUS
28#
29
30def sgn(x):
31    if x==0:
32        return 0
33    else:
34        return x/abs(x)
35
36def along(f, x, v):
37    """\
38Given a multivariate function f, a point x and a vector v,
39return the univariate function t |-> f(x+tv).
40    """
41    return lambda t,f=f,x=x,v=v: apply(f, add(x, multiply(t, v)))
42
43#
44# UNIVARIATE ROOT FINDING
45#
46
47def bracket_root(f, interval, max_iterations=50):
48    """\
49Given a univariate function f and a tuple interval=(x1,x2),
50return a new tuple (bracket, fnvals) where bracket=(x1,x2)
51brackets a root of f and fnvals=(f(x1),f(x2)).
52    """
53    x1, x2 = interval
54    if x1==x2:
55        raise BracketingException("initial interval has zero width")
56    elif x2<x1:
57        x1, x2 = x2, x1
58    f1, f2 = f(x1), f(x2)
59    for j in range(max_iterations):
60        while f1*f2 >= 0:  # not currently bracketed
61            if abs(f1)<abs(f2):
62                x1 = x1 + GOLDEN*(x1-x2)
63            else:
64                x2 = x2 + GOLDEN*(x2-x1)
65            f1, f2 = f(x1), f(x2)
66        return (x1, x2), (f1, f2)
67    raise BracketingException("too many iterations")
68
69def ridder_root(f, bracket, fnvals=None, accuracy=1e-6, max_iterations=50):
70    """\
71Given a univariate function f and a tuple bracket=(x1,x2) bracketing a root,
72find a root x of f using Ridder's method. Parameter fnvals=(f(x1),f(x2)) is optional.
73    """
74    x1, x2 = bracket
75    if fnvals==None:
76        f1, f2 = f(x1), f(x2)
77    else:
78        f1, f2 = fnvals
79    if f1==0:
80        return x1
81    elif f2==0:
82        return x2
83    elif f1*f2>=0:
84        raise BracketingException("initial interval does not bracket a root")
85    x4 = 123456789.
86    for j in range(max_iterations):
87        x3 = (x1+x2)/2
88        f3 = f(x3)
89        temp = f3*f3 - f1*f2
90        x4, x4old = x3 + (x3-x1)*sgn(f1-f2)*f3/temp**.5, x4
91        f4 = f(x4)
92        if f1*f4<0:  # x1 and x4 bracket root
93            x2, f2 = x4, f4
94        else:  # x4 and x2 bracket root
95            x1, f1 = x4, f4
96        if min(abs(x1-x2),abs(x4-x4old))<accuracy or temp==0:
97            return x4
98    raise RootFindingException("too many iterations")
99
100def root(f, interval=(0.,1.), accuracy=1e-4, max_iterations=50):
101    """\
102Given a univariate function f and an optional interval (x1,x2),
103find a root of f using bracket_root and ridder_root.
104    """
105    bracket, fnvals = bracket_root(f, interval, max_iterations)
106    return ridder_root(f, bracket, fnvals, accuracy, max_iterations)
107
108#
109# UNIVARIATE MINIMIZATION
110#
111
112def bracket_min(f, interval, max_iterations=50):
113    """\
114Given a univariate function f and a tuple interval=(x1,x2),
115return a new tuple (bracket, fnval) where bracket=(x1,x2,x3)
116brackets a minimum of f and fnvals=(f(x1),f(x2),f(x3)).
117    """
118    x1, x2 = interval
119    f1, f2 = f(x1), f(x2)
120    if f2>f1:  # ensure x1 --> x2 is downhill direction
121        x1, x2, f1, f2 = x2, x1, f2, f1
122    x3 = x2 + GOLDEN*(x2-x1)
123    f3 = f(x3)
124    for j in range(max_iterations):
125        if f2<f3:
126            if x1>x3:  # ensure x1<x2<x3
127                x1, x3, f1, f3 = x3, x1, f3, f1
128            return (x1, x2, x3), (f1, f2, f3)
129        else:
130            x1, x2, x3 = x1, x3, x3 + GOLDEN*(x3-x1)
131            f1, f2, f3 = f(x1), f(x2), f(x3)
132    raise BracketingException("too many iterations")
133
134def brent_min(f, bracket, fnvals=None, tolerance=1e-6, max_iterations=50):
135    """\
136Given a univariate function f and a tuple bracket=(x1,x2,x3) bracketing a minimum,
137find a local minimum of f (with fn value) using Brent's method.
138Optionally pass in the tuple fnvals=(f(x1),f(x2),f(x3)) as a parameter.
139    """
140    x1, x2, x3 = bracket
141    if fnvals==None:
142        f1, f2, f3 = f(x1), f(xx), f(x3)
143    else:
144        f1, f2, f3 = fnvals
145    if not f1>f2<f3:
146        raise MinimizationException("initial triple does not bracket a minimum")
147    if not x1<x3:  # ensure x1, x2, x3 in ascending order
148        x1, f1, x3, f3 = x3, f3, x1, f1
149    a, b = x1, x3
150
151    e = 0.
152    x = w = v = x2
153    fw = fv = fx = f(x)
154
155    for j in range(max_iterations):
156        xm = (a+b)/2
157        accuracy = tolerance*abs(x) + LITTLE
158        if abs(x-xm) < (2*accuracy - (b-a)/2):
159            return x, fx
160
161        if abs(e)>accuracy:
162            r = (x-w)*(fx-fv)
163            q = (x-v)*(fx-fw)
164            p = (x-v)*q - (x-w)*r
165            q = 2*(q-r)
166            if q>0:
167                p = -p
168            q = abs(q)
169            etemp = e
170            e = d
171            if abs(p)>=abs(q*etemp)/2 or p<=q*(a-x) or p>=q*(b-x):
172                if x>=xm:
173                    e = a-x
174                else:
175                    e = b-x
176                d = (2-GOLDEN)*e
177            else:  # accept parabolic fit
178                d = p/q
179                u = x+d
180                if u-a<2*accuracy or b-u<2*accuracy:
181                    d = accuracy*sgn(xm-x)
182        else:
183            if x>=xm:
184                e = a-x
185            else:
186                e = b-x
187            d = (2-GOLDEN)*e
188
189        if abs(d)>=accuracy:
190            u = x+d
191        else:
192            u = x+accuracy*sgn(d)
193        fu = f(u)
194
195        if fu<=fx:
196            if u>=x:
197                a = x
198            else:
199                b = x
200            v, w, x = w, x, u
201            fv, fw, fx = fw, fx, fu
202        else:
203            if u<x:
204                a = u
205            else:
206                b = u
207            if fu<-fw or w==x:
208                v, w, fv, fw = w, u, fw, fu
209            elif fu<=fw or v==x or v==w:
210                v, fv = u, fu
211
212    raise MinimizationException("too many iterations")
213
214def minimum(f, interval=(0.,1.), tolerance=1e-4, max_iterations=50, return_fnval=0):
215    """\
216Given a univariate function f and an optional interval (x1,x2),
217find a local minimum of f using bracket_min and brent_min.
218    """
219    bracket, fnvals = bracket_min(f, interval, max_iterations)
220    min, fnval = brent_min(f, bracket, fnvals, tolerance, max_iterations)
221    if return_fnval:
222        return min, fnval
223    else:
224        return min
225
226#
227# MULTIVARIATE MINIMIZATION
228#
229
230def powell_min(f, p0, tolerance=1e-4, max_iterations=200, return_fnval=0):
231    """\
232Given a multivariate function f and a starting point p0,
233find a local minimum of f using Powell's direction set method.
234    """
235    p = p0
236    fp = apply(f, p)
237    n = len(p)
238    directions = identity(n).tolist()  # the n coordinate vectors
239
240    for j in range(max_iterations):
241        maxdrop_i = 0
242        maxdrop = 0.
243        pold=p
244        fpold=fp
245
246        for i in range(n):
247            fptemp = fp
248            v = directions[i]
249            t, fp = minimum(along(f, p, v), return_fnval=1)
250            p = add(p, multiply(t, v))
251            if fptemp-fp>maxdrop:
252                maxdrop_i = i
253                maxdrop = fpold-fp
254
255        totaldrop = fpold-fp
256        if 2*totaldrop<=tolerance*(abs(fp)+abs(fpold)):
257            if return_fnval:
258                return tuple(p), fp
259            else:
260                return tuple(p)
261
262        vnew = subtract(p, pold)
263        pex = add(pold, multiply(2, vnew))
264        fpex = apply(f, pex)
265        if fpex<fp:
266            temp = fpold-fp-maxdrop
267            if 2*(fpold-2*fp+fpex)*temp*temp<totaldrop*totaldrop*maxdrop:
268                t, fp = minimum(along(f, p, vnew), return_fnval=1)
269                p = add(p, multiply(t, vnew))
270                directions[maxdrop_i] = vnew
271
272    raise MinimizationException("too many iterations")
273
274if __name__=='__main__':
275
276    from math import *
277
278    print root(cos)
279    print minimum(sin)
280
281    def f(x, y):
282        return x*x + x*y + y*y + 2*x - y
283
284    print powell_min(f, (0., 0.))
Note: See TracBrowser for help on using the repository browser.