source: sasview/park-1.2.1/park/simplex.py @ e9bd127

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since e9bd127 was 3570545, checked in by Mathieu Doucet <doucetm@…>, 13 years ago

Adding park Part 2

  • Property mode set to 100755
File size: 10.3 KB
Line 
1#__docformat__ = "restructuredtext en"
2# ******NOTICE***************
3# from optimize.py module by Travis E. Oliphant
4#
5# You may copy and use this module as you see fit with no
6# guarantee implied provided you keep this notice in all copies.
7# *****END NOTICE************
8#
9# Modified by Paul Kienzle to support bounded minimization
10"""
11Downhill simplex optimizer.
12"""
13
14__all__ = ['simplex']
15
16__docformat__ = "restructuredtext en"
17
18import numpy
19__version__="0.7"
20
21def wrap_function(function, bounds):
22    ncalls = [0]
23    if bounds is not None:
24        lo, hi = [numpy.asarray(v) for v in bounds]
25        def function_wrapper(x):
26            ncalls[0] += 1
27            if numpy.any((x<lo)|(x>hi)):
28                return numpy.inf
29            else:
30                return function(x)
31    else:
32        def function_wrapper(x):
33            ncalls[0] += 1
34            return function(x)
35    return ncalls, function_wrapper
36
37class Result:
38    """
39    Results from the fit.
40   
41    x : ndarray
42        Best parameter set
43    fx : float
44        Best value
45    iters : int
46        Number of iterations
47    calls : int
48        Number of function calls
49    success : boolean
50        True if the fit completed successful, false if terminated early
51        because of too many iterations.
52    """
53    def __init__(self, x, fx, iters, calls, status):
54        self.x,self.fx,self.iters,self.calls=x,fx,iters,calls
55        self.status = status
56    def __str__(self):
57        return "Minimum %g at %s after %d calls"%(self.fx,self.x,self.calls)
58
59def dont_abort(): return False
60
61def simplex(f, x0=None, bounds=None, radius=0.05,
62            xtol=1e-4, ftol=1e-4, maxiter=None,
63            update_handler=None, abort_test=dont_abort):
64    """
65    Minimize a function using Nelder-Mead downhill simplex algorithm.
66
67    This optimizer is also known as Amoeba (from Numerical Recipes) and
68    the Nealder-Mead simplex algorithm.  This is not the simplex algorithm
69    for solving constrained linear systems.
70
71    Downhill simplex is a robust derivative free algorithm for finding
72    minima.  It proceeds by choosing a set of points (the simplex) forming
73    an n-dimensional triangle, and transforming that triangle so that the
74    worst vertex is improved, either by stretching, shrinking or reflecting
75    it about the center of the triangle.  This algorithm is not known for
76    its speed, but for its simplicity and robustness, and is a good algorithm
77    to start your problem with.
78
79    *Parameters*:
80
81        f : callable f(x,*args)
82            The objective function to be minimized.
83        x0 : ndarray
84            Initial guess.
85        bounds : (ndarray,ndarray) or None
86            Bounds on the parameter values for the function.
87        radius: float
88            Size of the initial simplex.  For bounded parameters (those
89            which have finite lower and upper bounds), radius is clipped
90            to a value in (0,0.5] representing the portion of the
91            range to use as the size of the initial simplex.
92
93    *Returns*: Result (`park.simplex.Result`)
94
95        x : ndarray
96            Parameter that minimizes function.
97        fx : float
98            Value of function at minimum: ``fopt = func(xopt)``.
99        iters : int
100            Number of iterations performed.
101        calls : int
102            Number of function calls made.
103        success : boolean
104            True if fit completed successfully.
105
106    *Other Parameters*:
107
108        xtol : float
109            Relative error in xopt acceptable for convergence.
110        ftol : number
111            Relative error in func(xopt) acceptable for convergence.
112        maxiter : int=200*N
113            Maximum number of iterations to perform.  Defaults
114        update_handler : callable
115            Called after each iteration, as callback(xk,fxk), where xk
116            is the current parameter vector and fxk is the function value.
117            Returns True if the fit should continue.
118
119
120    *Notes*
121
122        Uses a Nelder-Mead simplex algorithm to find the minimum of
123        function of one or more variables.
124
125    """
126    fcalls, func = wrap_function(f, bounds)
127    x0 = numpy.asfarray(x0).flatten()
128    #print "x0",x0
129    N = len(x0)
130    rank = len(x0.shape)
131    if not -1 < rank < 2:
132        raise ValueError, "Initial guess must be a scalar or rank-1 sequence."
133
134    if maxiter is None:
135        maxiter = N * 200
136
137    rho = 1; chi = 2; psi = 0.5; sigma = 0.5;
138
139    if rank == 0:
140        sim = numpy.zeros((N+1,), dtype=x0.dtype)
141    else:
142        sim = numpy.zeros((N+1,N), dtype=x0.dtype)
143    fsim = numpy.zeros((N+1,), float)
144    sim[0] = x0
145    fsim[0] = func(x0)
146
147    # Metropolitan simplex: simplex has vertices at x0 and at
148    # x0 + j*radius for each unit vector j.  Radius is a percentage
149    # change from the initial value, or just the radius if the initial
150    # value is 0.  For bounded problems, the radius is a percentage of
151    # the bounded range in dimension j.
152    val = x0*(1+radius)
153    val[val == 0] = radius
154    if bounds is not None:
155        radius = numpy.clip(radius,0,0.5)
156        lo,hi = [numpy.asarray(v) for v in bounds]
157
158        # Keep the initial simplex inside the bounds
159        x0[x0<lo] = lo
160        x0[x0>hi] = hi
161        bounded = ~numpy.isinf(lo) & ~numpy.isinf(hi)
162        val[bounded] = x0[bounded] + (hi[bounded]-lo[bounded])*radius
163        val[val<lo] = lo
164        val[val>hi] = hi
165       
166        # If the initial point was at or beyond an upper bound, then bounds
167        # projection will put x0 and x0+j*radius at the same point.  We
168        # need to detect these collisions and reverse the radius step
169        # direction when such collisions occur.  The only time the collision
170        # can occur at the lower bound is when upper and lower bound are
171        # identical.  In that case, we are already done.
172        collision = val==x0
173        if numpy.any(collision):
174            reverse = x0*(1-radius)
175            reverse[reverse==0] = -radius
176            reverse[bounded] = x0[bounded] - (hi[bounded]-lo[bounded])*radius
177            val[collision] = reverse[collision]
178       
179        # Make tolerance relative for bounded parameters
180        tol = numpy.ones(x0.shape)*xtol
181        tol[bounded] = (hi[bounded]-lo[bounded])*xtol
182        xtol = tol
183
184    # Compute values at the simplex vertices
185    for k in range(0,N):
186        y = x0+0
187        y[k] = val[k]
188        sim[k+1] = y
189        fsim[k+1] = func(y)
190
191    #print sim
192    ind = numpy.argsort(fsim)
193    fsim = numpy.take(fsim,ind,0)
194    # sort so sim[0,:] has the lowest function value
195    sim = numpy.take(sim,ind,0)
196    #print sim
197
198    iterations = 1
199    while iterations < maxiter:
200        if numpy.all(abs(sim[1:]-sim[0]) <= xtol) \
201            and max(abs(fsim[0]-fsim[1:])) <= ftol:
202            #print abs(sim[1:]-sim[0])
203            break
204
205        xbar = numpy.sum(sim[:-1],0) / N
206        xr = (1+rho)*xbar - rho*sim[-1]
207        #print "xbar" ,xbar,rho,sim[-1],N
208        #break
209        fxr = func(xr)
210        doshrink = 0
211
212        if fxr < fsim[0]:
213            xe = (1+rho*chi)*xbar - rho*chi*sim[-1]
214            fxe = func(xe)
215
216            if fxe < fxr:
217                sim[-1] = xe
218                fsim[-1] = fxe
219            else:
220                sim[-1] = xr
221                fsim[-1] = fxr
222        else: # fsim[0] <= fxr
223            if fxr < fsim[-2]:
224                sim[-1] = xr
225                fsim[-1] = fxr
226            else: # fxr >= fsim[-2]
227                # Perform contraction
228                if fxr < fsim[-1]:
229                    xc = (1+psi*rho)*xbar - psi*rho*sim[-1]
230                    fxc = func(xc)
231
232                    if fxc <= fxr:
233                        sim[-1] = xc
234                        fsim[-1] = fxc
235                    else:
236                        doshrink=1
237                else:
238                    # Perform an inside contraction
239                    xcc = (1-psi)*xbar + psi*sim[-1]
240                    fxcc = func(xcc)
241
242                    if fxcc < fsim[-1]:
243                        sim[-1] = xcc
244                        fsim[-1] = fxcc
245                    else:
246                        doshrink = 1
247
248                if doshrink:
249                    for j in xrange(1,N+1):
250                        sim[j] = sim[0] + sigma*(sim[j] - sim[0])
251                        fsim[j] = func(sim[j])
252
253        ind = numpy.argsort(fsim)
254        sim = numpy.take(sim,ind,0)
255        fsim = numpy.take(fsim,ind,0)
256        if update_handler is not None:
257            update_handler(sim[0],fsim[0])
258        iterations += 1
259        if abort_test(): break
260
261    status = 0 if iterations < maxiter else 1
262    res = Result(sim[0],fsim[0],iterations,fcalls[0], status)
263    return res
264
265def main():
266    import time
267    def rosen(x):  # The Rosenbrock function
268        return numpy.sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0,axis=0)
269
270
271    x0 = [0.8,1.2,0.7]
272    print "Nelder-Mead Simplex"
273    print "==================="
274    start = time.time()
275    x = simplex(rosen,x0)
276    print x
277    print "Time:",time.time() - start
278
279    x0 = [0]*3
280    print "Nelder-Mead Simplex"
281    print "==================="
282    print "starting at zero"
283    start = time.time()
284    x = simplex(rosen,x0)
285    print x
286    print "Time:",time.time() - start
287
288    x0 = [0.8,1.2,0.7]
289    lo,hi = [0]*3, [1]*3
290    print "Bounded Nelder-Mead Simplex"
291    print "==========================="
292    start = time.time()
293    x = simplex(rosen,x0,bounds=(lo,hi))
294    print x
295    print "Time:",time.time() - start
296
297
298    x0 = [0.8,1.2,0.7]
299    lo,hi = [0.999]*3, [1.001]*3
300    print "Bounded Nelder-Mead Simplex"
301    print "==========================="
302    print "tight bounds"
303    print "simplex is smaller than 1e-7 in every dimension, but you can't"
304    print "see this without uncommenting the print statement simplex function"
305    start = time.time()
306    x = simplex(rosen,x0,bounds=(lo,hi),xtol=1e-4)
307    print x
308    print "Time:",time.time() - start
309
310
311    x0 = [0]*3
312    hi,lo = [-0.999]*3, [-1.001]*3
313    print "Bounded Nelder-Mead Simplex"
314    print "==========================="
315    print "tight bounds, x0=0 outside bounds from above"
316    start = time.time()
317    x = simplex(lambda x:rosen(-x),x0,bounds=(lo,hi),xtol=1e-4)
318    print x
319    print "Time:",time.time() - start
320
321
322    x0 = [0.8,1.2,0.7]
323    lo,hi = [-numpy.inf]*3, [numpy.inf]*3
324    print "Bounded Nelder-Mead Simplex"
325    print "==========================="
326    print "infinite bounds"
327    start = time.time()
328    x = simplex(rosen,x0,bounds=(lo,hi),xtol=1e-4)
329    print x
330    print "Time:",time.time() - start
331
332if __name__ == "__main__":
333    main()
Note: See TracBrowser for help on using the repository browser.