63b39f368a5b3b9b2199384e5c016a2809e3cf3d
[scilab.git] / scilab / modules / optimization / macros / neldermead / neldermead_search.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 2009 - INRIA - Michael Baudin
3 //
4 // This file must be used under the terms of the CeCILL.
5 // This source file is licensed as described in the file COPYING, which
6 // you should have received as part of this distribution.  The terms
7 // are also available at
8 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
9
10 //
11 // References
12 //   Sequential Application of Simplex Designs in Optimisation and Evolutionary Operation
13 //   Spendley, W. and Hext, G. R. and Himsworth, F. R.
14 //   1962
15 //
16 //   Convergence Properties of the Nelder-Mead Simplex Method in Low Dimensions
17 //   Jeffrey C. Lagarias and James A. Reeds nd Margaret H. Wright and Paul E. Wright
18 //   SIAM Journal of Optimization, 1998, volume 9, pp. 112--147
19 //
20 //   A Simplex Method for Function Minimization
21 //   Nelder, J. A.  and Mead, R.
22 //   1965
23 //
24 //   Iterative Methods for Optimization
25 //   C. T. Kelley,
26 //   SIAM Frontiers in Applied Mathematics
27 //   1999
28 //
29 //   Detection and Remediation of Stagnation in the Nelder--Mead Algorithm Using a Sufficient Decrease Condition
30 //   Kelley C. T.
31 //   SIAM J. on Optimization},
32 //   1999
33 //
34
35 //
36 // neldermead_search --
37 //   Search the minimum with Nelder-Mead algorithm.
38 //
39 function this = neldermead_search ( this )
40   if ( this.startupflag == 0) then
41     this = neldermead_startup ( this );
42     this.startupflag = 1;
43   end
44   neldermead_outputcmd ( this, "init" , this.simplex0 )
45   if this.restartflag == 1 then
46     this = neldermead_autorestart ( this )
47   else
48     this = neldermead_algo ( this );
49   end
50   neldermead_outputcmd ( this, "done" , this.simplexopt )
51 endfunction
52 //
53 // neldermead_algo --
54 //   Performs an optimization without restart
55 //
56 function this = neldermead_algo ( this )
57     select this.method
58     case "fixed" then
59       this = neldermead_fixed (this);
60     case "variable" then
61       this = neldermead_variable (this);
62     case "box" then
63       this = neldermead_box (this);
64     else
65       errmsg = msprintf(gettext("%s: Unknown -method %s"), "neldermead_algo", this.method)
66       error(errmsg)
67     end
68 endfunction
69 //
70 // neldermead_autorestart --
71 //   Performs an optimization with automatic restart
72 //
73 function this = neldermead_autorestart ( this )
74   for irestart = 1: this.restartmax
75     this = neldermead_log (this,sprintf("Restart #%d/%d", irestart,this.restartmax));
76     this = neldermead_algo ( this );
77     [ this , istorestart ] = neldermead_istorestart ( this );
78     if istorestart==0 then
79       this = neldermead_log (this,"Must not restart");
80       this.restartnb  = irestart
81       break
82     else
83       this = neldermead_log (this,"Must restart");
84     end
85     if ( irestart == this.restartmax ) then
86       this = neldermead_log (this,"Stopping after all restarts performed");
87       this.restartnb  = irestart
88       this.optbase = optimbase_set ( this.optbase , "-status" , "maxrestart" );
89     else
90       this = neldermead_updatesimp ( this );
91     end
92   end
93 endfunction
94
95
96 //
97 // neldermead_variable --
98 //   The original Nelder-Mead algorithm, with variable-size simplex.
99 //
100 function this = neldermead_variable ( this )
101   //
102   // Order the vertices for the first time
103   //
104   simplex = this.simplex0;
105   n = optimbase_cget ( this.optbase , "-numberofvariables" );
106   if (n==0) then
107     errmsg = msprintf(gettext("%s: The number of variable is zero."), "neldermead_variable")
108     error(errmsg)
109   end
110   fvinitial = optimbase_get ( this.optbase , "-fx0" );
111   // Sort function values and x points by increasing function value order
112   this = neldermead_log (this,"Step #1 : order");
113   simplex = optimsimplex_sort ( simplex );
114   currentcenter = optimsimplex_center ( simplex );
115   currentxopt = optimbase_cget ( this.optbase , "-x0" );
116   newfvmean = optimsimplex_fvmean ( simplex );
117   //
118   // Initialize
119   //
120   terminate = 0;
121   iter = 0;
122   //
123   // Nelder-Mead Loop
124   //
125   while ( terminate == 0 )
126     this.optbase = optimbase_incriter ( this.optbase );
127     iter = iter + 1;
128     xlow = optimsimplex_getx ( simplex , 1 )
129     flow = optimsimplex_getfv ( simplex , 1 )
130     xhigh = optimsimplex_getx ( simplex , n+1 )
131     fhigh = optimsimplex_getfv ( simplex , n+1 )
132     xn = optimsimplex_getx ( simplex , n )
133     fn = optimsimplex_getfv ( simplex , n )
134     //
135     // Store history
136     //
137     xcoords = optimsimplex_getallx ( simplex )
138     this = neldermead_storehistory ( this , n , flow , xlow , xcoords );
139     deltafv = abs(optimsimplex_deltafvmax ( simplex ));
140     currentfopt = flow;
141     previousxopt = currentxopt;
142     currentxopt = xlow;
143     previouscenter = currentcenter;
144     currentcenter = optimsimplex_center ( simplex );
145     oldfvmean = newfvmean;
146     newfvmean = optimsimplex_fvmean ( simplex );
147     totaliter = optimbase_get ( this.optbase , "-iterations" );
148     funevals = optimbase_get ( this.optbase , "-funevals" );
149     ssize = optimsimplex_size ( simplex )
150     this = neldermead_log (this,sprintf("================================================================="));
151     this = neldermead_log (this,sprintf("Iteration #%d (total = %d)",iter,totaliter));
152     this = neldermead_log (this,sprintf("Function Eval #%d",funevals));
153     this = neldermead_log (this,sprintf("Xopt : %s",strcat(string(xlow)," ")));
154     this = neldermead_log (this,sprintf("Fopt : %e",flow));
155     this = neldermead_log (this,sprintf("DeltaFv : %e",deltafv));
156     this = neldermead_log (this,sprintf("Center : %s",strcat(string(currentcenter)," ")));
157     this = neldermead_log (this,sprintf("Size : %e",ssize));
158     str = optimsimplex_tostring ( simplex )
159     for i = 1:n+1
160       this = neldermead_log (this,str(i));
161     end
162     this.optbase = optimbase_set ( this.optbase , "-xopt" , xlow );
163     this.optbase = optimbase_set ( this.optbase , "-fopt" , flow );
164     neldermead_outputcmd ( this, "iter" , simplex )
165
166     //
167     // Update termination flag
168     //
169     if ( iter > 1 ) then
170       [this , terminate , status] = neldermead_termination (this , ...
171         fvinitial , oldfvmean , newfvmean , previouscenter , currentcenter , simplex );
172       if (terminate==1) then
173         this = neldermead_log (this,sprintf("Terminate with status : %s",status));
174         break
175       end
176     end
177     //
178     // Compute xbar, center of better vertices
179     //
180     this = neldermead_log (this,sprintf("Reflect"));
181     xbar   = optimsimplex_xbar ( simplex );
182     this = neldermead_log (this,sprintf("xbar="+strcat(string(xbar)," ")+""));
183     //
184     // Reflect the worst point with respect to center
185     //
186     xr = neldermead_interpolate ( xbar , xhigh , this.rho );
187     [ this , fr ] = neldermead_function ( this , xr );
188     this = neldermead_log (this,sprintf("xr=["+strcat(string(xr)," ")+"], f(xr)=%f",fr));
189     if ( fr >= flow & fr < fn ) then
190       this = neldermead_log (this,sprintf("  > Perform reflection"));
191       simplex = optimsimplex_setve ( simplex , n+1 , fr , xr )
192     elseif ( fr < flow ) then
193       // Expand
194       this = neldermead_log (this,sprintf("Expand"));
195       xe = neldermead_interpolate ( xbar , xhigh , this.rho*this.chi );
196       [ this ,fe] = neldermead_function ( this ,xe);
197       this = neldermead_log (this,sprintf("xe="+strcat(string(xe)," ")+", f(xe)=%f",fe));
198       if (fe < fr) then
199         this = neldermead_log (this,sprintf("  > Perform Expansion"));
200         simplex = optimsimplex_setve ( simplex , n+1 , fe , xe )
201       else
202         this = neldermead_log (this,sprintf("  > Perform reflection"));
203         simplex = optimsimplex_setve ( simplex , n+1 , fr , xr )
204       end
205     elseif ( fr >= fn & fr < fhigh ) then
206       // Outside contraction
207       this = neldermead_log (this,sprintf("Contract - outside"));
208       xc = neldermead_interpolate ( xbar , xhigh , this.rho*this.gamma );
209       [ this ,fc] = neldermead_function ( this ,xc);
210       this = neldermead_log (this,sprintf("xc="+strcat(string(xc)," ")+", f(xc)=%f",fc));
211       if ( fc <= fr ) then
212         this = neldermead_log (this,sprintf("  > Perform Outside Contraction"));
213         simplex = optimsimplex_setve ( simplex , n+1 , fc , xc )
214       else
215         //  Shrink
216         this = neldermead_log (this,sprintf("  > Perform Shrink"));
217         [ simplex , this ] = optimsimplex_shrink ( simplex , neldermead_costf , this.sigma , this );
218       end
219     else
220       // ( fr >= fn & fr >= fhigh )  
221       // Inside contraction
222       this = neldermead_log (this,sprintf("Contract - inside"));
223       xc = neldermead_interpolate ( xbar , xhigh , -this.gamma );
224       [ this ,fc] = neldermead_function ( this ,xc);
225       this = neldermead_log (this,sprintf("xc="+strcat(string(xc)," ")+", f(xc)=%f",fc));
226       if ( fc < fhigh ) then
227         this = neldermead_log (this,sprintf("  > Perform Inside Contraction"));
228         simplex = optimsimplex_setve ( simplex , n+1 , fc , xc )
229       else
230         //  Shrink
231         this = neldermead_log (this,sprintf("  > Perform Shrink"));
232         [ simplex , this ] = optimsimplex_shrink ( simplex , neldermead_costf , this.sigma , this )
233       end
234     end
235     //
236     // Sort simplex
237     //
238     this = neldermead_log (this,sprintf("Sort"));
239     simplex  = optimsimplex_sort ( simplex );
240   end
241   this.optbase = optimbase_set ( this.optbase , "-xopt" , xlow.' );
242   this.optbase = optimbase_set ( this.optbase , "-fopt" , flow );
243   this.optbase = optimbase_set ( this.optbase , "-status" , status );
244   this.simplexopt = simplex;
245 endfunction
246
247 //
248 // neldermead_fixed --
249 //   The simplex algorithm with fixed size simplex.
250 //
251 function this = neldermead_fixed (this)
252   //
253   // Order the vertices for the first time
254   //
255   simplex = this.simplex0;
256   n = optimbase_cget ( this.optbase , "-numberofvariables" );
257   fvinitial = optimbase_get ( this.optbase , "-fx0" );
258   // Sort function values and x points by increasing function value order
259   this = neldermead_log (this,sprintf("Sort"));
260   simplex = optimsimplex_sort ( simplex );
261   neldermead_outputcmd ( this, "init" , simplex )
262   //
263   // Compute center of simplex
264   //
265   currentcenter = optimsimplex_center ( simplex );
266   newfvmean = optimsimplex_fvmean ( simplex );
267   currentxopt = optimbase_cget ( this.optbase , "-x0" );
268   //
269   // Initialize
270   //
271   terminate = 0;
272   iter = 0;
273   //
274   // main N-M loop
275   //
276   while (terminate == 0)
277     this.optbase = optimbase_incriter ( this.optbase );
278     iter = iter + 1;
279     xlow = optimsimplex_getx ( simplex , 1 )
280     flow = optimsimplex_getfv ( simplex , 1 )
281     xhigh = optimsimplex_getx ( simplex , n+1 )
282     fhigh = optimsimplex_getfv ( simplex , n+1 )
283     //
284     // Store history
285     //
286     xcoords = optimsimplex_getallx ( simplex )
287     this = neldermead_storehistory ( this , n , flow , xlow , xcoords );
288     deltafv = abs(optimsimplex_deltafvmax ( simplex ));
289     currentfopt = flow;
290     previousxopt = currentxopt;
291     currentxopt = xlow;
292     previouscenter = currentcenter;
293     currentcenter = optimsimplex_center ( simplex );
294     oldfvmean = newfvmean;
295     newfvmean = optimsimplex_fvmean ( simplex );
296     totaliter = optimbase_get ( this.optbase , "-iterations" );
297     funevals = optimbase_get ( this.optbase , "-funevals" );
298     ssize = optimsimplex_size ( simplex )
299     this = neldermead_log (this,sprintf("================================================================="));
300     this = neldermead_log (this,sprintf("Iteration #%d (total = %d)",iter,totaliter));
301     this = neldermead_log (this,sprintf("Function Eval #%d",funevals));
302     this = neldermead_log (this,sprintf("Xopt : %s",strcat(string(xlow)," ")));
303     this = neldermead_log (this,sprintf("Fopt : %e",flow));
304     this = neldermead_log (this,sprintf("DeltaFv : %e",deltafv));
305     this = neldermead_log (this,sprintf("Center : %s",strcat(string(currentcenter)," ")));
306     this = neldermead_log (this,sprintf("Size : %e",ssize));
307     str = optimsimplex_tostring ( simplex )
308     for i = 1:n+1
309       this = neldermead_log (this,str(i));
310     end
311     this.optbase = optimbase_set ( this.optbase , "-xopt" , xlow );
312     this.optbase = optimbase_set ( this.optbase , "-fopt" , flow );
313     neldermead_outputcmd ( this, "iter" , simplex )
314     //
315     // Update termination flag
316     //
317     if ( iter > 1 ) then
318       [this , terminate , status] = neldermead_termination (this , ...
319         fvinitial , oldfvmean , newfvmean , previouscenter , currentcenter , simplex );
320       if (terminate==1) then
321         this = neldermead_log (this,sprintf("Terminate with status : %s",status));
322         break;
323       end
324     end
325     //
326     // Compute xbar, center of better vertices
327     //
328     this = neldermead_log (this,sprintf("Reflect"));
329     xbar = optimsimplex_xbar ( simplex );
330     this = neldermead_log (this,sprintf("xbar="+strcat(string(xbar)," ")+""));
331     //
332     // Reflect the worst point with respect to center
333     //
334     xr = neldermead_interpolate ( xbar , xhigh , this.rho );
335     [ this ,fr] = neldermead_function ( this ,xr);
336     this = neldermead_log (this,sprintf("xr="+strcat(string(xr)," ")+", f(xr)=%f",fr));
337     //
338     // Replace worst point by xr if it is better
339     //
340     if ( fr < fhigh ) then
341       this = neldermead_log (this,sprintf("  > Perform reflect"));
342       simplex = optimsimplex_setve ( simplex , n+1 , fr , xr )
343     else
344       //  Shrink
345       this = neldermead_log (this,sprintf("  > Perform Shrink"));
346       [ simplex , this ] = optimsimplex_shrink ( simplex , neldermead_costf , this.sigma , this )
347     end
348     //
349     // Sort simplex
350     //
351     simplex = optimsimplex_sort ( simplex );
352   end
353   this.optbase = optimbase_set ( this.optbase , "-xopt" , xlow.' );
354   this.optbase = optimbase_set ( this.optbase , "-fopt" , flow );
355   this.optbase = optimbase_set ( this.optbase , "-status" , status );
356   this.simplexopt = simplex;
357 endfunction
358 //
359 // neldermead_interpolate --
360 //   Computes xi as an interpolation between x1 and x2, with factor as :
361 //     xi = (1+fac)x1 - fac x2
362 //
363 function xi = neldermead_interpolate ( x1 , x2 , fac )
364   xi = (1 + fac)*x1 - fac*x2;
365 endfunction
366
367 //
368 // neldermead_termination --
369 //   Returns 1 if the algorithm terminates.
370 //   Returns 0 if the algorithm must continue.
371 // Arguments
372 //   this : the current object
373 //   fvinitial : initial function value
374 //   newfvmean, oldfvmean : the old and new function value average on the simplex
375 //   previousxopt : the previous value of x
376 //   currentxopt : the current value of x
377 //   simplex : the simplex
378 //     The best point in the simplex is expected to be stored at 1
379 //     The worst point in the simplex is expected to be stored at n+1
380 //   terminate : 1 if the algorithm terminates, 0 if the algorithm must continue.
381 //   this.status : termination status
382 //     status = "continue"
383 //     status = "maxiter"
384 //     status = "maxfuneval"
385 //     status = "tolf"
386 //     status = "tolx"
387 //     status = "tolfstdev"
388 //     status = "tolsize"
389 //     status = "tolsizedeltafv"
390 // Notes
391 //   Use the function average on the simplex instead of the best function value.
392 //   This is because the function average changes at each iteration.
393 //   Instead, the best function value as a step-by-step evolution and may not
394 //   change in 2 iterations, leading to astop of the algorithm.
395 // TODO : set the fvinitial, oldfvmean, newfvmean.
396 //
397 function [this , terminate , status ] = neldermead_termination (this , ...
398   fvinitial , oldfvmean , newfvmean , previousxopt , currentxopt , ...
399   simplex )
400   terminate = 0;
401   status = "continue";
402   //
403   // Termination Criteria from parent optimization class
404   //
405   [ this.optbase , terminate , status] = optimbase_terminate ( this.optbase , ...
406     fvinitial , newfvmean , previousxopt , currentxopt );
407   //
408   // Criteria #5 : standard deviation of function values
409   //
410   if (terminate == 0) then
411     if this.tolfstdeviationmethod == "enabled" then
412       fv = optimsimplex_getallfv ( simplex )
413       sd = st_deviation(fv);
414       this.optbase = optimbase_stoplog  ( this.optbase,sprintf("  > st_deviation(fv)=%e < tolfstdeviation=%e",...
415         sd, this.tolfstdeviation));
416       if sd < this.tolfstdeviation then
417         terminate = 1;
418         status = "tolfstdev";
419       end
420     end
421   end
422   //
423   // Criteria #6 : simplex absolute + relative size
424   //
425   if (terminate == 0) then
426     if this.tolsimplexizemethod == "enabled" then
427       ssize = optimsimplex_size ( simplex , "sigmaplus" );
428       this.optbase = optimbase_stoplog  ( this.optbase,sprintf("  > simplex size=%e < %e + %e * %e",...
429         ssize, this.tolsimplexizeabsolute , this.tolsimplexizerelative , this.simplexsize0 ));
430       if ssize < this.tolsimplexizeabsolute + this.tolsimplexizerelative * this.simplexsize0 then
431         terminate = 1;
432         status = "tolsize";
433       end
434     end
435   end
436   //
437   // Criteria #7 : simplex absolute size + difference in function values (Matlab-like)
438   //
439   if (terminate == 0) then
440     if this.tolssizedeltafvmethod == "enabled" then
441       ssize = optimsimplex_size ( simplex , "sigmaplus" );
442       this.optbase = optimbase_stoplog  ( this.optbase,sprintf("  > simplex size=%e < %e",...
443         ssize, this.tolsimplexizeabsolute));
444       shiftfv = abs(optimsimplex_deltafvmax( simplex ))
445       this.optbase = optimbase_stoplog  ( this.optbase,sprintf("  > abs(fv(n+1) - fv(1))=%e < toldeltafv=%e",...
446         shiftfv, this.toldeltafv));
447       if ssize < this.tolsimplexizeabsolute & shiftfv < this.toldeltafv then
448         terminate = 1;
449         status = "tolsizedeltafv";
450       end
451     end
452   end
453   //
454   // Criteria #8 : Kelley stagnation, based on
455   // a sufficient decrease condition
456   //
457   if ( terminate == 0 ) then
458     if ( this.kelleystagnationflag==1 ) then
459       [ sg , this ] = optimsimplex_gradientfv ( simplex , neldermead_costf , "forward" , this );
460       nsg = sg.' * sg;
461       sgstr = strcat(string(sg)," ");
462       this.optbase = optimbase_stoplog ( this.optbase , sprintf ( "Test Stagnation : nsg = %e, sg = "+sgstr, nsg) );
463       this.optbase = optimbase_stoplog ( this.optbase , ...
464         sprintf ( "Test Stagnation : newfvmean=%e >= oldfvmean=%e - %e * %e" , newfvmean, oldfvmean , this.kelleyalpha , nsg ) );
465       if ( newfvmean >= oldfvmean - this.kelleyalpha * nsg ) then
466         terminate = 1;
467         status = "kelleystagnation";
468       end
469     end
470   end
471 endfunction
472   
473
474 //
475 // neldermead_outputcmd --
476 //   Calls back user's output command
477 // Arguments
478 //   this : the current object
479 //   state : the state of the algorithm,
480 //     "init", "done", "iter"
481 //   simplex : the current simplex
482 //
483 function  neldermead_outputcmd ( this, ...
484    state , simplex )
485   outputcmd = optimbase_cget ( this.optbase , "-outputcommand" );
486   if typeof(outputcmd) <> "string" then
487     brutedata = optimbase_outstruct ( this.optbase );
488     data = tlist(["T_NMDATA",...
489       "x","fval","iteration","funccount",...
490       "simplex"]);
491     data.x = brutedata.x;
492     data.fval = brutedata.fval;
493     data.iteration = brutedata.iteration;
494     data.funccount = brutedata.funccount;
495     data.simplex = simplex;
496     optimbase_outputcmd ( this.optbase , state , data );
497   end
498 endfunction
499 //
500 // neldermead_storehistory --
501 //   Stores the history into the data structure.
502 // Arguments, input
503 //   this : current object
504 //   n : the number of unknown parameters
505 //   fopt : the current value of the function at optimum
506 //   xopt : arrary with size n, current optimum
507 //   xcoords : array with size n x n+1, coordinates of the n+1 vertices
508 //
509 function this = neldermead_storehistory ( this , n , fopt , xopt , xcoords )
510   storehistory = optimbase_cget ( this.optbase , "-storehistory" );
511   iterations = optimbase_get ( this.optbase , "-iterations" );
512   if storehistory == 1 then
513     this.optbase = optimbase_histset ( this.optbase , iterations , "-fopt" , fopt );
514     this.optbase = optimbase_histset ( this.optbase , iterations , "-xopt" , xopt(1:n).' );
515     this.historysimplex ( iterations , 1:n+1,1:n) = xcoords(1:n+1,1:n);
516   end
517 endfunction
518
519 //
520 // neldermead_istorestart --
521 //   Returns 1 if the optimization is to restart.
522 // Arguments
523 //   istorestart : 1 of the the optimization is to restart.
524 //
525 function [ this , istorestart ] = neldermead_istorestart ( this )
526   select this.restartdetection
527   case "oneill"
528     [ this , istorestart ] = neldermead_isroneill ( this )
529   case "kelley"
530     [ this , istorestart ] = neldermead_isrkelley ( this )
531   else
532     errmsg = msprintf(gettext("%s: Unknown restart detection %s"),"neldermead_istorestart", this.restartdetection)
533     error(errmsg)
534   end
535 endfunction
536 //
537 // neldermead_isrkelley--
538 //   Returns 1 if the optimization is to restart.
539 //   Use kelleystagnation as a criteria for restart.
540 // Arguments
541 //   istorestart : 1 of the the optimization is to restart.
542 //
543 function [ this , istorestart ] = neldermead_isrkelley ( this )
544   istorestart = 0
545   if ( this.kelleystagnationflag==1 ) then
546     status = optimbase_get ( this.optbase , "-status" );
547     if ( status =="kelleystagnation" ) then
548        istorestart = 1
549     end
550   end
551 endfunction
552 //
553 // neldermead_isroneill --
554 //   Returns 1 if the optimization is to restart.
555 //   Use O'Neill method as a criteria for restart.
556 //   It is a axis by axis search for optimality.
557 // Arguments
558 //   xopt : the optimum, as a st of n values
559 //   fopt : function value at optimum
560 //   eps : a small value
561 //   step : a list of n values, representing
562 //     the "size" of each parameter
563 //   istorestart : 1 of the the optimization is to restart.
564 //
565 function [ this , istorestart ] = neldermead_isroneill ( this )
566   n = optimbase_cget ( this.optbase , "-numberofvariables" );
567   //
568   // If required, make a vector step from the scalar step
569   //
570   defaultstep = this.restartstep;
571   stepn = length ( defaultstep );
572   if ( stepn <> n ) then
573     step = defaultstep * ones(n,1);
574   else
575     step = defaultstep;
576   end
577
578   xopt = optimbase_get ( this.optbase , "-xopt" );
579   fopt = optimbase_get ( this.optbase , "-fopt" );
580
581     istorestart = 0
582     for ix = 1:n
583       stepix = step ( ix )
584       del = stepix * this.restarteps
585       if ( del==0.0 ) then
586          del = eps
587       end
588       xoptix =  xopt ( ix )
589       xopt ( ix ) = xoptix + del
590       [ this , fv ] = neldermead_function ( this , xopt )
591       if ( fv < fopt ) then
592         istorestart = 1
593         break
594       end
595       xopt ( ix ) = xoptix - del
596       [ this , fv ] = neldermead_function ( this , xopt )
597       if ( fv < fopt ) then
598         istorestart = 1
599         break
600       end
601       xopt ( ix ) = xoptix
602     end
603 endfunction
604
605 //
606 // neldermead_startup --
607 //   Startup the algorithm.
608 //   Computes the initial simplex, depending on the -simplex0method.
609 //
610 function this = neldermead_startup (this)
611   [ this.optbase , hasbounds ] = optimbase_hasbounds ( this.optbase );
612   if ( hasbounds ) then
613     [ this.optbase , isok , errmsg ] = optimbase_checkbounds ( this.optbase );
614     if ( ~isok ) then
615       error ( msprintf(gettext("%s: %s"), "neldermead_startup" , errmsg ))
616     end
617   end
618   x0 = optimbase_cget ( this.optbase , "-x0" );
619   select this.simplex0method
620   case "given" then
621     [ simplex0 , this ] = optimsimplex_new ( this.coords0 , ...
622       neldermead_costf , this );
623   case "axes" then
624     simplex0 = optimsimplex_new ( );
625     [ simplex0 , this ] = optimsimplex_axes ( simplex0 , ...
626       x0.' , neldermead_costf , this.simplex0length , this );
627   case "spendley" then
628     simplex0 = optimsimplex_new ( );
629     [ simplex0 , this ] = optimsimplex_spendley ( simplex0 , ...
630       x0.' , neldermead_costf , this.simplex0length , this );
631   case "pfeffer" then
632     simplex0 = optimsimplex_new ( );
633     [ simplex0 , this ] = optimsimplex_pfeffer ( simplex0 , ...
634       x0.' , neldermead_costf , this.simplex0deltausual , ...
635       this.simplex0deltazero , this );
636   case "randbounds" then
637     simplex0 = optimsimplex_new ( );
638     if ( this.boxnbpoints == "2n" ) then
639       this.boxnbpointseff = 2 * this.optbase.numberofvariables;
640     else
641       this.boxnbpointseff = this.boxnbpoints;
642     end
643     if ( ~hasbounds ) then
644       error ( msprintf(gettext("%s: Randomized bounds initial simplex is not available without bounds." ), "neldermead_startup"))
645     end
646     [ simplex0 , this ] = optimsimplex_randbounds ( simplex0 , x0.' , ...
647       neldermead_costf , this.optbase.boundsmin , this.optbase.boundsmax , ...
648       this.boxnbpointseff  , this );
649   else
650     errmsg = msprintf(gettext("%s: Unknown -simplex0method : %s"), "neldermead_startup", this.simplex0method);
651     error(errmsg);
652   end
653   //
654   // Scale the simplex into the bounds and the nonlinear inequality constraints, if any
655   //
656   if ( hasbounds | this.optbase.nbineqconst > 0 ) then
657     this = neldermead_log (this,sprintf("Scaling initial simplex into nonlinear inequality constraints..."));
658     nbve = optimsimplex_getnbve ( simplex0 );
659     for ive = 1 : nbve
660       // optimsimplex returns a row vector
661       x = optimsimplex_getx ( simplex0 , ive );
662       this = neldermead_log (this,sprintf("Scaling vertex #%d/%d at ["+...
663         strcat(string(x)," ")+"]... " , ...
664         ive , nbve ));
665       // Transpose x, because x0 is a column vector
666       [ this , status , xp ] = _scaleinconstraints ( this , x.' , x0 );
667       if ( ~status ) then
668         errmsg = msprintf(gettext("%s: Impossible to scale the vertex #%d/%d at [%s] into inequality constraints"), ...
669           "neldermead_startup", ive , nbve , strcat(string(x)," "));
670         error(errmsg);
671       end
672       if ( or ( x <> xp ) ) then
673         [ this , fv ] = neldermead_function ( this , xp );
674         // Transpose xp, which is a column vector
675         simplex0 = optimsimplex_setve ( simplex0 , ive , fv , xp.' );
676       end
677     end
678   end
679   //
680   // Store the simplex
681   //
682   this.simplex0 = optimsimplex_destroy ( this.simplex0 );
683   this.simplex0 = simplex0;
684   this.simplexsize0 = optimsimplex_size ( simplex0 );
685   fx0 = optimsimplex_getfv ( this.simplex0 , 1 );
686   this.optbase = optimbase_set ( this.optbase , "-fx0" , fx0 );
687   this.optbase = optimbase_set ( this.optbase , "-xopt" , x0.' );
688   this.optbase = optimbase_set ( this.optbase , "-fopt" , fx0 );
689   this.optbase = optimbase_set ( this.optbase , "-iterations" , 0 );
690   if ( this.kelleystagnationflag == 1 ) then
691     this = neldermead_kelleystag ( this );
692   end
693 endfunction
694
695 //
696 // neldermead_kelleystag --
697 //   Initialize Kelley's stagnation detection system when normalization is required,
698 //   by computing kelleyalpha.
699 //   If the simplex gradient is zero, then
700 //   use alpha0 as alpha.
701 // Arguments
702 //   status : the status after the failing
703 //     optimization process
704 //   simplex : the simplex computed at the end of the failing
705 //     optimization process
706 //
707 function this = neldermead_kelleystag ( this )
708     if this.kelleystagnationflag == 1 then
709       if this.kelleynormalizationflag == 0 then
710         this.kelleyalpha = this.kelleystagnationalpha0
711       else
712         [sg,this] = optimsimplex_gradientfv ( this.simplex0 , neldermead_costf , "forward" , this )
713         nsg = sg.' * sg
714         sigma0 = optimsimplex_size ( this.simplex0 , "sigmaplus" )
715         if nsg==0.0 then 
716           this.kelleyalpha = this.kelleystagnationalpha0
717         else
718           this.kelleyalpha = this.kelleystagnationalpha0 * sigma0 / nsg
719         end
720       end
721       this = neldermead_log (this,sprintf("Test Stagnation Kelley : alpha0 = %e", this.kelleyalpha));
722     end
723 endfunction
724   //
725   // _scaleinconstraints --
726   //   Given a point to scale and a reference point which satisfies the constraints, 
727   //   scale the point towards the reference point until it satisfies all the constraints.
728   //   Returns a list of key values with the following
729   //   keys : -status 0/1 -x x, where status
730   //   is 0 if the procedure has failed after -boxnbnlloops
731   //   iterations.
732   // Arguments
733   //   x : the point to scale
734   //   xref : the reference point
735   //   status : %T or %F
736   //   p : scaled point
737   //
738 function [ this , status , p ] = _scaleinconstraints ( this , x , xref )
739   p = x
740   //
741   // Project the point into the bounds
742   //
743   [ this.optbase , hasbounds ] = optimbase_hasbounds ( this.optbase );
744   if ( hasbounds ) then
745     [ this.optbase , p ] = optimbase_proj2bnds ( this.optbase ,  p );
746     this = neldermead_log (this,sprintf(" > After projection into bounds p = [%s]" , ...
747       strcat(string(p)," ")));
748   end
749   //
750   // Adjust point to satisfy nonlinear inequality constraints
751   //
752   nbnlc = optimbase_cget ( this.optbase , "-nbineqconst" )
753   if ( nbnlc == 0 ) then
754     scaled = %T
755   else
756     scaled = %F
757     //
758     // Try the current point and see if the constraints are satisfied.
759     // If not, move the point "halfway" to the centroid,
760     // which should satisfy the constraints, if
761     // the constraints are convex.
762     // Perform this loop until the constraints are satisfied.
763     // If all loops have been performed without success, the scaling
764     // has failed.
765     //
766     for i = 1 : this.nbineqloops
767       [ this , constlist ] = neldermead_function ( this , p , index=2 );
768       feasible = %T
769       for ic = 1 : this.optbase.nbineqconst;
770         const = constlist ( ic );
771         if ( const < 0.0 ) then
772           this = neldermead_log (this,sprintf("Constraint #%d/%d is not satisfied", ...
773             ic , this.optbase.nbineqconst ));
774           feasible = %F;
775           break;
776         end
777       end
778       if ( feasible ) then
779         scaled = %T;
780         break;
781       else
782         this = neldermead_log (this,sprintf("Scaling inequality constraint at loop #%d/%d", ...
783           i , this.nbineqloops));
784         p = ( xref + p ) * this.ineqscaling;
785       end
786     end
787     this = neldermead_log (this,sprintf(" > After scaling into inequality constraints p = [%s]" , ...
788       strcat(string(p)," ") ) );
789   end
790   if ( scaled ) then
791     status = %T
792   else
793     status = %F
794     this = neldermead_log (this,sprintf(" > Impossible to scale into constraints after %d loops" , ...
795       this.optbase.nbineqconst ));
796   end
797 endfunction
798 //
799 // neldermead_box --
800 //   The Nelder-Mead algorithm, with variable-size simplex
801 //   and modifications by Box for bounds and
802 //   inequality constraints.
803 //
804 function this = neldermead_box ( this )
805   //
806   // Order the vertices for the first time
807   //
808   simplex = this.simplex0;
809   n = optimbase_cget ( this.optbase , "-numberofvariables" );
810   fvinitial = optimbase_get ( this.optbase , "-fx0" );
811   // Sort function values and x points by increasing function value order
812   this = neldermead_log (this,"Step #1 : order");
813   simplex = optimsimplex_sort ( simplex );
814   currentcenter = optimsimplex_center ( simplex );
815   currentxopt = optimbase_cget ( this.optbase , "-x0" );
816   newfvmean = optimsimplex_fvmean ( simplex );
817   nbve = optimsimplex_getnbve ( simplex );
818   ihigh = nbve;
819   inext = ihigh - 1
820   ilow = 1
821   [ this.optbase , hasbounds ] = optimbase_hasbounds ( this.optbase );
822   nbnlc = optimbase_cget ( this.optbase , "-nbineqconst" )
823   //
824   // Initialize
825   //
826   terminate = 0;
827   iter = 0;
828   //
829   // Nelder-Mead Loop
830   //
831   while ( terminate == 0 )
832     this.optbase = optimbase_incriter ( this.optbase );
833     iter = iter + 1;
834     xlow = optimsimplex_getx ( simplex , ilow )
835     flow = optimsimplex_getfv ( simplex , ilow )
836     xhigh = optimsimplex_getx ( simplex , ihigh )
837     fhigh = optimsimplex_getfv ( simplex , ihigh )
838     xn = optimsimplex_getx ( simplex , inext )
839     fn = optimsimplex_getfv ( simplex , inext )
840     //
841     // Store history
842     //
843     xcoords = optimsimplex_getallx ( simplex )
844     this = neldermead_storehistory ( this , n , flow , xlow , xcoords );
845     deltafv = abs(optimsimplex_deltafvmax ( simplex ));
846     currentfopt = flow;
847     previousxopt = currentxopt;
848     currentxopt = xlow;
849     previouscenter = currentcenter;
850     currentcenter = optimsimplex_center ( simplex );
851     oldfvmean = newfvmean;
852     newfvmean = optimsimplex_fvmean ( simplex );
853     totaliter = optimbase_get ( this.optbase , "-iterations" );
854     funevals = optimbase_get ( this.optbase , "-funevals" );
855     ssize = optimsimplex_size ( simplex )
856     this = neldermead_log (this,sprintf("================================================================="));
857     this = neldermead_log (this,sprintf("Iteration #%d (total = %d)",iter,totaliter));
858     this = neldermead_log (this,sprintf("Function Eval #%d",funevals));
859     this = neldermead_log (this,sprintf("Xopt : [%s]",strcat(string(xlow)," ")));
860     this = neldermead_log (this,sprintf("Fopt : %e",flow));
861     this = neldermead_log (this,sprintf("DeltaFv : %e",deltafv));
862     this = neldermead_log (this,sprintf("Center : [%s]",strcat(string(currentcenter)," ")));
863     this = neldermead_log (this,sprintf("Size : %e",ssize));
864     str = optimsimplex_tostring ( simplex )
865     for i = 1:nbve
866       this = neldermead_log (this,str(i));
867     end
868     neldermead_outputcmd ( this, "iter" , simplex )
869
870     //
871     // Update termination flag
872     //
873     if ( iter > 1 ) then
874       [this , terminate , status] = neldermead_termination (this , ...
875         fvinitial , oldfvmean , newfvmean , previouscenter , currentcenter , simplex );
876       if (terminate==1) then
877         this = neldermead_log (this,sprintf("Terminate with status : %s",status));
878         break
879       end
880     end
881     //
882     // Compute xbar, center of better vertices
883     //
884     this = neldermead_log (this,sprintf("Reflect"));
885     xbar = optimsimplex_xbar ( simplex );
886     this = neldermead_log (this,sprintf("xbar=[%s]" , strcat(string(xbar)," ")));
887     //
888     // Reflect the worst point with respect to center
889     //
890     xr = neldermead_interpolate ( xbar , xhigh , this.rho );
891     // Adjust point to satisfy bounds and nonlinear inequality constraints
892     if ( hasbounds | nbnlc > 0 ) then
893       [ this , status , xr ] = _scaleinconstraints ( this , xr , xbar )
894       if ( ~status ) then
895         status = "impossibleconstr"
896         break
897       end
898     end
899     [ this , fr ] = neldermead_function ( this , xr );
900     this = neldermead_log (this,sprintf("xr=[%s], f(xr)=%f", strcat(string(xr)," ") , fr));
901     if ( fr >= flow & fr < fn ) then
902       this = neldermead_log (this,sprintf("  > Perform reflection"));
903       simplex = optimsimplex_setve ( simplex , ihigh , fr , xr )
904     elseif ( fr < flow ) then
905       // Expand
906       this = neldermead_log (this,sprintf("Expand"));
907       xe = neldermead_interpolate ( xbar , xhigh , this.rho*this.chi );
908       // Adjust point to satisfy bounds and nonlinear inequality constraints
909       if ( hasbounds | nbnlc > 0 ) then
910         [ this , status , xe ] = _scaleinconstraints ( this , xe , xbar )
911         if ( ~status ) then
912           status = "impossibleconstr"
913           break
914         end
915       end
916       [ this ,fe] = neldermead_function ( this ,xe);
917       this = neldermead_log (this,sprintf("xe=[%s], f(xe)=%f", strcat(string(xe)," ") , fe ));
918       if (fe < fr) then
919         this = neldermead_log (this,sprintf("  > Perform Expansion"));
920         simplex = optimsimplex_setve ( simplex , ihigh , fe , xe )
921       else
922         this = neldermead_log (this,sprintf("  > Perform reflection"));
923         simplex = optimsimplex_setve ( simplex , ihigh , fr , xr )
924       end
925     elseif ( fr >= fn & fr < fhigh ) then
926       // Outside contraction
927       this = neldermead_log (this,sprintf("Contract - outside"));
928       xc = neldermead_interpolate ( xbar , xhigh , this.rho*this.gamma );
929       // Adjust point to satisfy bounds and nonlinear inequality constraints
930       if ( hasbounds | nbnlc > 0 ) then
931         [ this , status , xc ] = _scaleinconstraints ( this , xc , xbar )
932         if ( ~status ) then
933           status = "impossibleconstr"
934           break
935         end
936       end
937       [ this ,fc] = neldermead_function ( this ,xc);
938       this = neldermead_log (this,sprintf("xc=[%s], f(xc)=%f", strcat(string(xc)," ") , fc));
939       if ( fc <= fr ) then
940         this = neldermead_log (this,sprintf("  > Perform Outside Contraction"));
941         simplex = optimsimplex_setve ( simplex , ihigh , fc , xc )
942       else
943         //  Shrink
944         this = neldermead_log (this,sprintf("  > Perform Shrink"));
945         [ simplex , this ] = optimsimplex_shrink ( simplex , neldermead_costf , this.sigma , this );
946       end
947     else
948       // ( fr >= fn & fr >= fhigh )  
949       // Inside contraction
950       this = neldermead_log (this,sprintf("Contract - inside"));
951       xc = neldermead_interpolate ( xbar , xhigh , -this.gamma );
952       // Adjust point to satisfy bounds and nonlinear inequality constraints
953       if ( hasbounds | nbnlc > 0 ) then
954         [ this , status , xc ] = _scaleinconstraints ( this , xc , xbar )
955         if ( ~status ) then
956           status = "impossibleconstr"
957           break
958         end
959       end
960       [ this ,fc] = neldermead_function ( this ,xc);
961       this = neldermead_log (this,sprintf("xc=[%s], f(xc)=%f", strcat(string(xc)," ") , fc));
962       if ( fc < fhigh ) then
963         this = neldermead_log (this,sprintf("  > Perform Inside Contraction"));
964         simplex = optimsimplex_setve ( simplex , ihigh , fc , xc )
965       else
966         //  Shrink
967         this = neldermead_log (this,sprintf("  > Perform Shrink"));
968         [ simplex , this ] = optimsimplex_shrink ( simplex , neldermead_costf , this.sigma , this )
969       end
970     end
971     //
972     // Sort simplex
973     //
974     this = neldermead_log (this,sprintf("Sort"));
975     simplex  = optimsimplex_sort ( simplex );
976   end
977   this.optbase = optimbase_set ( this.optbase , "-xopt" , xlow.' );
978   this.optbase = optimbase_set ( this.optbase , "-fopt" , flow );
979   this.optbase = optimbase_set ( this.optbase , "-status" , status );
980   this.simplexopt = simplex;
981 endfunction
982