Deleted vectorized computation feature. Deleted neldermead_contour. Fixed the demos.
[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 // Implementation note:
251 //   We implement the following "rules" of the Spendley et al.
252 //   method.
253 //   Rule 1 is strictly applied, but the reflection is done
254 //   by reflection the high point, since we minimize a function
255 //   instead of maximizing it, like Spendley.
256 //   The Rule 2 is NOT implemented, as we expect that the
257 //   function evaluation is not subject to errors.
258 //   The Rule 3 is applied, ie reflection with respect
259 //   to next to high point.
260 //   A shrink step is included, with shrinkage factor sigma.
261 //
262 //   "Rule 1. Ascertain the lowest reading y, of yi ... Yk+1
263 //   Complete a new simplex Sp by excluding the point Vp corresponding to
264 //   y, and replacing it by V* defined as above."
265 //   "Rule 2. If a result has occurred in (k + 1) successive simplexes, and is not
266 //   then eliminated by application of Rule 1, do not move in the direction
267 //   indicated by Rule 1, or at all, but discard the result and replace it by a new
268 //   observation at the same point."
269 //   "Rule 3. If y is the lowest reading in So , and if the next observation made,
270 //   y* , is the lowest reading in the new simplex S , do not apply Rule 1 and
271 //   return to So from Sp . Move out of S, by rejecting the second lowest reading
272 //   (which is also the second lowest reading in So)."
273 //
274 function this = neldermead_fixed (this)
275   //
276   // Order the vertices for the first time
277   //
278   simplex = this.simplex0;
279   n = optimbase_cget ( this.optbase , "-numberofvariables" );
280   fvinitial = optimbase_get ( this.optbase , "-fx0" );
281   // Sort function values and x points by increasing function value order
282   this = neldermead_log (this,sprintf("Sort"));
283   simplex = optimsimplex_sort ( simplex );
284   neldermead_outputcmd ( this, "init" , simplex )
285   //
286   // Compute center of simplex
287   //
288   currentcenter = optimsimplex_center ( simplex );
289   newfvmean = optimsimplex_fvmean ( simplex );
290   currentxopt = optimbase_cget ( this.optbase , "-x0" );
291   //
292   // Set indices for "clarity"
293   //
294   ilow = 1
295   ihigh = n + 1
296   inext = n
297   //
298   // Initialize
299   //
300   terminate = 0;
301   iter = 0;
302   //
303   // main N-M loop
304   //
305   while (terminate == 0)
306     this.optbase = optimbase_incriter ( this.optbase );
307     iter = iter + 1;
308     xlow = optimsimplex_getx ( simplex , ilow )
309     flow = optimsimplex_getfv ( simplex , ilow )
310     xhigh = optimsimplex_getx ( simplex , ihigh )
311     fhigh = optimsimplex_getfv ( simplex , ihigh )
312     //
313     // Store history
314     //
315     xcoords = optimsimplex_getallx ( simplex )
316     this = neldermead_storehistory ( this , n , flow , xlow , xcoords );
317     deltafv = abs(optimsimplex_deltafvmax ( simplex ));
318     currentfopt = flow;
319     previousxopt = currentxopt;
320     currentxopt = xlow;
321     previouscenter = currentcenter;
322     currentcenter = optimsimplex_center ( simplex );
323     oldfvmean = newfvmean;
324     newfvmean = optimsimplex_fvmean ( simplex );
325     totaliter = optimbase_get ( this.optbase , "-iterations" );
326     funevals = optimbase_get ( this.optbase , "-funevals" );
327     ssize = optimsimplex_size ( simplex )
328     this = neldermead_log (this,sprintf("================================================================="));
329     this = neldermead_log (this,sprintf("Iteration #%d (total = %d)",iter,totaliter));
330     this = neldermead_log (this,sprintf("Function Eval #%d",funevals));
331     this = neldermead_log (this,sprintf("Xopt : %s",strcat(string(xlow)," ")));
332     this = neldermead_log (this,sprintf("Fopt : %e",flow));
333     this = neldermead_log (this,sprintf("DeltaFv : %e",deltafv));
334     this = neldermead_log (this,sprintf("Center : %s",strcat(string(currentcenter)," ")));
335     this = neldermead_log (this,sprintf("Size : %e",ssize));
336     str = optimsimplex_tostring ( simplex )
337     for i = 1:n+1
338       this = neldermead_log (this,str(i));
339     end
340     this.optbase = optimbase_set ( this.optbase , "-xopt" , xlow );
341     this.optbase = optimbase_set ( this.optbase , "-fopt" , flow );
342     neldermead_outputcmd ( this, "iter" , simplex )
343     //
344     // Update termination flag
345     //
346     if ( iter > 1 ) then
347       [this , terminate , status] = neldermead_termination (this , ...
348         fvinitial , oldfvmean , newfvmean , previouscenter , currentcenter , simplex );
349       if (terminate==1) then
350         this = neldermead_log (this,sprintf("Terminate with status : %s",status));
351         break;
352       end
353     end
354     //
355     // Compute xbar, center of better vertices
356     //
357     this = neldermead_log (this,sprintf("Reflect"));
358     xbar = optimsimplex_xbar ( simplex );
359     this = neldermead_log (this,sprintf("xbar="+strcat(string(xbar)," ")+""));
360     //
361     // Reflect the worst point with respect to center
362     //
363     xr = neldermead_interpolate ( xbar , xhigh , this.rho );
364     [ this , fr ] = neldermead_function ( this , xr );
365     this = neldermead_log (this,sprintf("xr="+strcat(string(xr)," ")+", f(xr)=%f",fr));
366     //
367     // Replace worst point by xr if it is better
368     //
369     if ( fr < fhigh ) then
370       this = neldermead_log (this,sprintf("  > Perform reflect"));
371       simplex = optimsimplex_setve ( simplex , ihigh , fr , xr )
372     else
373       // Reflect / xnext
374       xnext = optimsimplex_getx ( simplex , inext );
375       fnext = optimsimplex_getfv ( simplex , inext );
376       xbar2 = optimsimplex_xbar ( simplex , inext );
377       this = neldermead_log (this,sprintf("xbar2="+strcat(string(xbar2)," ")+""));
378       xr2 = neldermead_interpolate ( xbar2 , xnext , this.rho );
379       [ this , fr2 ] = neldermead_function ( this ,xr2 );
380       this = neldermead_log (this,sprintf("xr2="+strcat(string(xr2)," ")+", f(xr2)=%f",fr2));
381       if ( fr2 < fnext ) then
382         this = neldermead_log (this,sprintf("  > Perform reflect / next"));
383         simplex = optimsimplex_setve ( simplex , inext , fr2 , xr2 )
384       else
385         //  Shrink
386         this = neldermead_log (this,sprintf("  > Perform Shrink"));
387         [ simplex , this ] = optimsimplex_shrink ( simplex , neldermead_costf , this.sigma , this )
388       end
389     end
390     //
391     // Sort simplex
392     //
393     simplex = optimsimplex_sort ( simplex );
394   end
395   this.optbase = optimbase_set ( this.optbase , "-xopt" , xlow.' );
396   this.optbase = optimbase_set ( this.optbase , "-fopt" , flow );
397   this.optbase = optimbase_set ( this.optbase , "-status" , status );
398   this.simplexopt = simplex;
399 endfunction
400 //
401 // neldermead_interpolate --
402 //   Computes xi as an interpolation between x1 and x2, with factor as :
403 //     xi = (1+fac)x1 - fac x2
404 //
405 function xi = neldermead_interpolate ( x1 , x2 , fac )
406   xi = (1 + fac)*x1 - fac*x2;
407 endfunction
408
409 //
410 // neldermead_termination --
411 //   Returns 1 if the algorithm terminates.
412 //   Returns 0 if the algorithm must continue.
413 // Arguments
414 //   this : the current object
415 //   fvinitial : initial function value
416 //   newfvmean, oldfvmean : the old and new function value average on the simplex
417 //   previousxopt : the previous value of x
418 //   currentxopt : the current value of x
419 //   simplex : the simplex
420 //     The best point in the simplex is expected to be stored at 1
421 //     The worst point in the simplex is expected to be stored at n+1
422 //   terminate : 1 if the algorithm terminates, 0 if the algorithm must continue.
423 //   this.status : termination status
424 //     status = "continue"
425 //     status = "maxiter"
426 //     status = "maxfuneval"
427 //     status = "tolf"
428 //     status = "tolx"
429 //     status = "tolfstdev"
430 //     status = "tolsize"
431 //     status = "tolsizedeltafv"
432 // Notes
433 //   Use the function average on the simplex instead of the best function value.
434 //   This is because the function average changes at each iteration.
435 //   Instead, the best function value as a step-by-step evolution and may not
436 //   change in 2 iterations, leading to astop of the algorithm.
437 // TODO : set the fvinitial, oldfvmean, newfvmean.
438 //
439 function [this , terminate , status ] = neldermead_termination (this , ...
440   fvinitial , oldfvmean , newfvmean , previousxopt , currentxopt , ...
441   simplex )
442   terminate = 0;
443   status = "continue";
444   //
445   // Termination Criteria from parent optimization class
446   //
447   [ this.optbase , terminate , status] = optimbase_terminate ( this.optbase , ...
448     fvinitial , newfvmean , previousxopt , currentxopt );
449   //
450   // Criteria #5 : standard deviation of function values
451   //
452   if (terminate == 0) then
453     if this.tolfstdeviationmethod == "enabled" then
454       fv = optimsimplex_getallfv ( simplex )
455       sd = st_deviation(fv);
456       this.optbase = optimbase_stoplog  ( this.optbase,sprintf("  > st_deviation(fv)=%e < tolfstdeviation=%e",...
457         sd, this.tolfstdeviation));
458       if sd < this.tolfstdeviation then
459         terminate = 1;
460         status = "tolfstdev";
461       end
462     end
463   end
464   //
465   // Criteria #6 : simplex absolute + relative size
466   //
467   if (terminate == 0) then
468     if this.tolsimplexizemethod == "enabled" then
469       ssize = optimsimplex_size ( simplex , "sigmaplus" );
470       this.optbase = optimbase_stoplog  ( this.optbase,sprintf("  > simplex size=%e < %e + %e * %e",...
471         ssize, this.tolsimplexizeabsolute , this.tolsimplexizerelative , this.simplexsize0 ));
472       if ssize < this.tolsimplexizeabsolute + this.tolsimplexizerelative * this.simplexsize0 then
473         terminate = 1;
474         status = "tolsize";
475       end
476     end
477   end
478   //
479   // Criteria #7 : simplex absolute size + difference in function values (Matlab-like)
480   //
481   if (terminate == 0) then
482     if this.tolssizedeltafvmethod == "enabled" then
483       ssize = optimsimplex_size ( simplex , "sigmaplus" );
484       this.optbase = optimbase_stoplog  ( this.optbase,sprintf("  > simplex size=%e < %e",...
485         ssize, this.tolsimplexizeabsolute));
486       shiftfv = abs(optimsimplex_deltafvmax( simplex ))
487       this.optbase = optimbase_stoplog  ( this.optbase,sprintf("  > abs(fv(n+1) - fv(1))=%e < toldeltafv=%e",...
488         shiftfv, this.toldeltafv));
489       if ssize < this.tolsimplexizeabsolute & shiftfv < this.toldeltafv then
490         terminate = 1;
491         status = "tolsizedeltafv";
492       end
493     end
494   end
495   //
496   // Criteria #8 : Kelley stagnation, based on
497   // a sufficient decrease condition
498   //
499   if ( terminate == 0 ) then
500     if ( this.kelleystagnationflag==1 ) then
501       [ sg , this ] = optimsimplex_gradientfv ( simplex , neldermead_costf , "forward" , this );
502       nsg = sg.' * sg;
503       sgstr = strcat(string(sg)," ");
504       this.optbase = optimbase_stoplog ( this.optbase , sprintf ( "Test Stagnation : nsg = %e, sg = "+sgstr, nsg) );
505       this.optbase = optimbase_stoplog ( this.optbase , ...
506         sprintf ( "Test Stagnation : newfvmean=%e >= oldfvmean=%e - %e * %e" , newfvmean, oldfvmean , this.kelleyalpha , nsg ) );
507       if ( newfvmean >= oldfvmean - this.kelleyalpha * nsg ) then
508         terminate = 1;
509         status = "kelleystagnation";
510       end
511     end
512   end
513 endfunction
514   
515
516 //
517 // neldermead_outputcmd --
518 //   Calls back user's output command
519 // Arguments
520 //   this : the current object
521 //   state : the state of the algorithm,
522 //     "init", "done", "iter"
523 //   simplex : the current simplex
524 //
525 function  neldermead_outputcmd ( this, ...
526    state , simplex )
527   outputcmd = optimbase_cget ( this.optbase , "-outputcommand" );
528   if typeof(outputcmd) <> "string" then
529     brutedata = optimbase_outstruct ( this.optbase );
530     data = tlist(["T_NMDATA",...
531       "x","fval","iteration","funccount",...
532       "simplex"]);
533     data.x = brutedata.x;
534     data.fval = brutedata.fval;
535     data.iteration = brutedata.iteration;
536     data.funccount = brutedata.funccount;
537     data.simplex = simplex;
538     optimbase_outputcmd ( this.optbase , state , data );
539   end
540 endfunction
541 //
542 // neldermead_storehistory --
543 //   Stores the history into the data structure.
544 // Arguments, input
545 //   this : current object
546 //   n : the number of unknown parameters
547 //   fopt : the current value of the function at optimum
548 //   xopt : arrary with size n, current optimum
549 //   xcoords : array with size n x n+1, coordinates of the n+1 vertices
550 //
551 function this = neldermead_storehistory ( this , n , fopt , xopt , xcoords )
552   storehistory = optimbase_cget ( this.optbase , "-storehistory" );
553   iterations = optimbase_get ( this.optbase , "-iterations" );
554   if storehistory == 1 then
555     this.optbase = optimbase_histset ( this.optbase , iterations , "-fopt" , fopt );
556     this.optbase = optimbase_histset ( this.optbase , iterations , "-xopt" , xopt(1:n).' );
557     this.historysimplex ( iterations , 1:n+1,1:n) = xcoords(1:n+1,1:n);
558   end
559 endfunction
560
561 //
562 // neldermead_istorestart --
563 //   Returns 1 if the optimization is to restart.
564 // Arguments
565 //   istorestart : 1 of the the optimization is to restart.
566 //
567 function [ this , istorestart ] = neldermead_istorestart ( this )
568   select this.restartdetection
569   case "oneill"
570     [ this , istorestart ] = neldermead_isroneill ( this )
571   case "kelley"
572     [ this , istorestart ] = neldermead_isrkelley ( this )
573   else
574     errmsg = msprintf(gettext("%s: Unknown restart detection %s"),"neldermead_istorestart", this.restartdetection)
575     error(errmsg)
576   end
577 endfunction
578 //
579 // neldermead_isrkelley--
580 //   Returns 1 if the optimization is to restart.
581 //   Use kelleystagnation as a criteria for restart.
582 // Arguments
583 //   istorestart : 1 of the the optimization is to restart.
584 //
585 function [ this , istorestart ] = neldermead_isrkelley ( this )
586   istorestart = 0
587   if ( this.kelleystagnationflag==1 ) then
588     status = optimbase_get ( this.optbase , "-status" );
589     if ( status =="kelleystagnation" ) then
590        istorestart = 1
591     end
592   end
593 endfunction
594 //
595 // neldermead_isroneill --
596 //   Returns 1 if the optimization is to restart.
597 //   Use O'Neill method as a criteria for restart.
598 //   It is a axis by axis search for optimality.
599 // Arguments
600 //   xopt : the optimum, as a st of n values
601 //   fopt : function value at optimum
602 //   eps : a small value
603 //   step : a list of n values, representing
604 //     the "size" of each parameter
605 //   istorestart : 1 of the the optimization is to restart.
606 //
607 function [ this , istorestart ] = neldermead_isroneill ( this )
608   n = optimbase_cget ( this.optbase , "-numberofvariables" );
609   //
610   // If required, make a vector step from the scalar step
611   //
612   defaultstep = this.restartstep;
613   stepn = length ( defaultstep );
614   if ( stepn <> n ) then
615     step = defaultstep * ones(n,1);
616   else
617     step = defaultstep;
618   end
619
620   xopt = optimbase_get ( this.optbase , "-xopt" );
621   fopt = optimbase_get ( this.optbase , "-fopt" );
622
623     istorestart = 0
624     for ix = 1:n
625       stepix = step ( ix )
626       del = stepix * this.restarteps
627       if ( del==0.0 ) then
628          del = eps
629       end
630       xoptix =  xopt ( ix )
631       xopt ( ix ) = xoptix + del
632       [ this , fv ] = neldermead_function ( this , xopt )
633       if ( fv < fopt ) then
634         istorestart = 1
635         break
636       end
637       xopt ( ix ) = xoptix - del
638       [ this , fv ] = neldermead_function ( this , xopt )
639       if ( fv < fopt ) then
640         istorestart = 1
641         break
642       end
643       xopt ( ix ) = xoptix
644     end
645 endfunction
646
647 //
648 // neldermead_startup --
649 //   Startup the algorithm.
650 //   Computes the initial simplex, depending on the -simplex0method.
651 //
652 function this = neldermead_startup (this)
653   [ this.optbase , hasbounds ] = optimbase_hasbounds ( this.optbase );
654   if ( hasbounds ) then
655     [ this.optbase , isok , errmsg ] = optimbase_checkbounds ( this.optbase );
656     if ( ~isok ) then
657       error ( msprintf(gettext("%s: %s"), "neldermead_startup" , errmsg ))
658     end
659   end
660   x0 = optimbase_cget ( this.optbase , "-x0" );
661   select this.simplex0method
662   case "given" then
663     [ simplex0 , this ] = optimsimplex_new ( this.coords0 , ...
664       neldermead_costf , this );
665   case "axes" then
666     simplex0 = optimsimplex_new ( );
667     [ simplex0 , this ] = optimsimplex_axes ( simplex0 , ...
668       x0.' , neldermead_costf , this.simplex0length , this );
669   case "spendley" then
670     simplex0 = optimsimplex_new ( );
671     [ simplex0 , this ] = optimsimplex_spendley ( simplex0 , ...
672       x0.' , neldermead_costf , this.simplex0length , this );
673   case "pfeffer" then
674     simplex0 = optimsimplex_new ( );
675     [ simplex0 , this ] = optimsimplex_pfeffer ( simplex0 , ...
676       x0.' , neldermead_costf , this.simplex0deltausual , ...
677       this.simplex0deltazero , this );
678   case "randbounds" then
679     simplex0 = optimsimplex_new ( );
680     if ( this.boxnbpoints == "2n" ) then
681       this.boxnbpointseff = 2 * this.optbase.numberofvariables;
682     else
683       this.boxnbpointseff = this.boxnbpoints;
684     end
685     if ( ~hasbounds ) then
686       error ( msprintf(gettext("%s: Randomized bounds initial simplex is not available without bounds." ), "neldermead_startup"))
687     end
688     [ simplex0 , this ] = optimsimplex_randbounds ( simplex0 , x0.' , ...
689       neldermead_costf , this.optbase.boundsmin , this.optbase.boundsmax , ...
690       this.boxnbpointseff  , this );
691   else
692     errmsg = msprintf(gettext("%s: Unknown -simplex0method : %s"), "neldermead_startup", this.simplex0method);
693     error(errmsg);
694   end
695   //
696   // Scale the simplex into the bounds and the nonlinear inequality constraints, if any
697   //
698   if ( hasbounds | this.optbase.nbineqconst > 0 ) then
699     this = neldermead_log (this,sprintf("Scaling initial simplex into nonlinear inequality constraints..."));
700     nbve = optimsimplex_getnbve ( simplex0 );
701     for ive = 1 : nbve
702       // optimsimplex returns a row vector
703       x = optimsimplex_getx ( simplex0 , ive );
704       this = neldermead_log (this,sprintf("Scaling vertex #%d/%d at ["+...
705         strcat(string(x)," ")+"]... " , ...
706         ive , nbve ));
707       // Transpose x, because x0 is a column vector
708       [ this , status , xp ] = _scaleinconstraints ( this , x.' , x0 );
709       if ( ~status ) then
710         errmsg = msprintf(gettext("%s: Impossible to scale the vertex #%d/%d at [%s] into inequality constraints"), ...
711           "neldermead_startup", ive , nbve , strcat(string(x)," "));
712         error(errmsg);
713       end
714       if ( or ( x <> xp ) ) then
715         // Set the index, so that, if an additionnal cost function argument is provided,
716         // it can be appended at the end.
717         index = 1;
718         [ this , fv ] = neldermead_function ( this , xp , index );
719         // Transpose xp, which is a column vector
720         simplex0 = optimsimplex_setve ( simplex0 , ive , fv , xp.' );
721       end
722     end
723   end
724   //
725   // Store the simplex
726   //
727   this.simplex0 = optimsimplex_destroy ( this.simplex0 );
728   this.simplex0 = simplex0;
729   this.simplexsize0 = optimsimplex_size ( simplex0 );
730   fx0 = optimsimplex_getfv ( this.simplex0 , 1 );
731   this.optbase = optimbase_set ( this.optbase , "-fx0" , fx0 );
732   this.optbase = optimbase_set ( this.optbase , "-xopt" , x0.' );
733   this.optbase = optimbase_set ( this.optbase , "-fopt" , fx0 );
734   this.optbase = optimbase_set ( this.optbase , "-iterations" , 0 );
735   if ( this.kelleystagnationflag == 1 ) then
736     this = neldermead_kelleystag ( this );
737   end
738 endfunction
739
740 //
741 // neldermead_kelleystag --
742 //   Initialize Kelley's stagnation detection system when normalization is required,
743 //   by computing kelleyalpha.
744 //   If the simplex gradient is zero, then
745 //   use alpha0 as alpha.
746 // Arguments
747 //   status : the status after the failing
748 //     optimization process
749 //   simplex : the simplex computed at the end of the failing
750 //     optimization process
751 //
752 function this = neldermead_kelleystag ( this )
753     if this.kelleystagnationflag == 1 then
754       if this.kelleynormalizationflag == 0 then
755         this.kelleyalpha = this.kelleystagnationalpha0
756       else
757         [sg,this] = optimsimplex_gradientfv ( this.simplex0 , neldermead_costf , "forward" , this )
758         nsg = sg.' * sg
759         sigma0 = optimsimplex_size ( this.simplex0 , "sigmaplus" )
760         if nsg==0.0 then 
761           this.kelleyalpha = this.kelleystagnationalpha0
762         else
763           this.kelleyalpha = this.kelleystagnationalpha0 * sigma0 / nsg
764         end
765       end
766       this = neldermead_log (this,sprintf("Test Stagnation Kelley : alpha0 = %e", this.kelleyalpha));
767     end
768 endfunction
769   //
770   // _scaleinconstraints --
771   //   Given a point to scale and a reference point which satisfies the constraints, 
772   //   scale the point towards the reference point until it satisfies all the constraints.
773   //   Returns a list of key values with the following
774   //   keys : -status 0/1 -x x, where status
775   //   is 0 if the procedure has failed after -boxnbnlloops
776   //   iterations.
777   // Arguments
778   //   x : the point to scale
779   //   xref : the reference point
780   //   status : %T or %F
781   //   p : scaled point
782   //
783 function [ this , status , p ] = _scaleinconstraints ( this , x , xref )
784   p = x
785   //
786   // Project the point into the bounds
787   //
788   [ this.optbase , hasbounds ] = optimbase_hasbounds ( this.optbase );
789   if ( hasbounds ) then
790     [ this.optbase , p ] = optimbase_proj2bnds ( this.optbase ,  p );
791     this = neldermead_log (this,sprintf(" > After projection into bounds p = [%s]" , ...
792       strcat(string(p)," ")));
793   end
794   //
795   // Adjust point to satisfy nonlinear inequality constraints
796   //
797   nbnlc = optimbase_cget ( this.optbase , "-nbineqconst" )
798   if ( nbnlc == 0 ) then
799     scaled = %T
800   else
801     scaled = %F
802     //
803     // Try the current point and see if the constraints are satisfied.
804     // If not, move the point "halfway" to the centroid,
805     // which should satisfy the constraints, if
806     // the constraints are convex.
807     // Perform this loop until the constraints are satisfied.
808     // If all loops have been performed without success, the scaling
809     // has failed.
810     //
811     for i = 1 : this.nbineqloops
812       [ this , constlist ] = neldermead_function ( this , p , index=2 );
813       feasible = %T
814       for ic = 1 : this.optbase.nbineqconst;
815         const = constlist ( ic );
816         if ( const < 0.0 ) then
817           this = neldermead_log (this,sprintf("Constraint #%d/%d is not satisfied", ...
818             ic , this.optbase.nbineqconst ));
819           feasible = %F;
820           break;
821         end
822       end
823       if ( feasible ) then
824         scaled = %T;
825         break;
826       else
827         this = neldermead_log (this,sprintf("Scaling inequality constraint at loop #%d/%d", ...
828           i , this.nbineqloops));
829         p = ( xref + p ) * this.ineqscaling;
830       end
831     end
832     this = neldermead_log (this,sprintf(" > After scaling into inequality constraints p = [%s]" , ...
833       strcat(string(p)," ") ) );
834   end
835   if ( scaled ) then
836     status = %T
837   else
838     status = %F
839     this = neldermead_log (this,sprintf(" > Impossible to scale into constraints after %d loops" , ...
840       this.optbase.nbineqconst ));
841   end
842 endfunction
843 //
844 // neldermead_box --
845 //   The Nelder-Mead algorithm, with variable-size simplex
846 //   and modifications by Box for bounds and
847 //   inequality constraints.
848 //
849 function this = neldermead_box ( this )
850   //
851   // Order the vertices for the first time
852   //
853   simplex = this.simplex0;
854   n = optimbase_cget ( this.optbase , "-numberofvariables" );
855   fvinitial = optimbase_get ( this.optbase , "-fx0" );
856   // Sort function values and x points by increasing function value order
857   this = neldermead_log (this,"Step #1 : order");
858   simplex = optimsimplex_sort ( simplex );
859   currentcenter = optimsimplex_center ( simplex );
860   currentxopt = optimbase_cget ( this.optbase , "-x0" );
861   newfvmean = optimsimplex_fvmean ( simplex );
862   nbve = optimsimplex_getnbve ( simplex );
863   ihigh = nbve;
864   inext = ihigh - 1
865   ilow = 1
866   [ this.optbase , hasbounds ] = optimbase_hasbounds ( this.optbase );
867   nbnlc = optimbase_cget ( this.optbase , "-nbineqconst" )
868   //
869   // Initialize
870   //
871   terminate = 0;
872   iter = 0;
873   //
874   // Nelder-Mead Loop
875   //
876   while ( terminate == 0 )
877     this.optbase = optimbase_incriter ( this.optbase );
878     iter = iter + 1;
879     xlow = optimsimplex_getx ( simplex , ilow )
880     flow = optimsimplex_getfv ( simplex , ilow )
881     xhigh = optimsimplex_getx ( simplex , ihigh )
882     fhigh = optimsimplex_getfv ( simplex , ihigh )
883     xn = optimsimplex_getx ( simplex , inext )
884     fn = optimsimplex_getfv ( simplex , inext )
885     //
886     // Store history
887     //
888     xcoords = optimsimplex_getallx ( simplex )
889     this = neldermead_storehistory ( this , n , flow , xlow , xcoords );
890     deltafv = abs(optimsimplex_deltafvmax ( simplex ));
891     currentfopt = flow;
892     previousxopt = currentxopt;
893     currentxopt = xlow;
894     previouscenter = currentcenter;
895     currentcenter = optimsimplex_center ( simplex );
896     oldfvmean = newfvmean;
897     newfvmean = optimsimplex_fvmean ( simplex );
898     totaliter = optimbase_get ( this.optbase , "-iterations" );
899     funevals = optimbase_get ( this.optbase , "-funevals" );
900     ssize = optimsimplex_size ( simplex )
901     this = neldermead_log (this,sprintf("================================================================="));
902     this = neldermead_log (this,sprintf("Iteration #%d (total = %d)",iter,totaliter));
903     this = neldermead_log (this,sprintf("Function Eval #%d",funevals));
904     this = neldermead_log (this,sprintf("Xopt : [%s]",strcat(string(xlow)," ")));
905     this = neldermead_log (this,sprintf("Fopt : %e",flow));
906     this = neldermead_log (this,sprintf("DeltaFv : %e",deltafv));
907     this = neldermead_log (this,sprintf("Center : [%s]",strcat(string(currentcenter)," ")));
908     this = neldermead_log (this,sprintf("Size : %e",ssize));
909     str = optimsimplex_tostring ( simplex )
910     for i = 1:nbve
911       this = neldermead_log (this,str(i));
912     end
913     neldermead_outputcmd ( this, "iter" , simplex )
914
915     //
916     // Update termination flag
917     //
918     if ( iter > 1 ) then
919       [this , terminate , status] = neldermead_termination (this , ...
920         fvinitial , oldfvmean , newfvmean , previouscenter , currentcenter , simplex );
921       if (terminate==1) then
922         this = neldermead_log (this,sprintf("Terminate with status : %s",status));
923         break
924       end
925     end
926     //
927     // Compute xbar, center of better vertices
928     //
929     this = neldermead_log (this,sprintf("Reflect"));
930     xbar = optimsimplex_xbar ( simplex );
931     this = neldermead_log (this,sprintf("xbar=[%s]" , strcat(string(xbar)," ")));
932     //
933     // Reflect the worst point with respect to center
934     //
935     xr = neldermead_interpolate ( xbar , xhigh , this.rho );
936     // Adjust point to satisfy bounds and nonlinear inequality constraints
937     if ( hasbounds | nbnlc > 0 ) then
938       [ this , status , xr ] = _scaleinconstraints ( this , xr , xbar )
939       if ( ~status ) then
940         status = "impossibleconstr"
941         break
942       end
943     end
944     if ( nbnlc > 0 ) then
945       // Set the index, so that, if an additionnal cost function argument is provided,
946       // it can be appended at the end.
947       index = 1;
948       [ this , fr ] = neldermead_function ( this , xr , index );
949     else
950       [ this , fr ] = neldermead_function ( this , xr );
951     end
952     this = neldermead_log (this,sprintf("xr=[%s], f(xr)=%f", strcat(string(xr)," ") , fr));
953     if ( fr >= flow & fr < fn ) then
954       this = neldermead_log (this,sprintf("  > Perform reflection"));
955       simplex = optimsimplex_setve ( simplex , ihigh , fr , xr )
956     elseif ( fr < flow ) then
957       // Expand
958       this = neldermead_log (this,sprintf("Expand"));
959       xe = neldermead_interpolate ( xbar , xhigh , this.rho*this.chi );
960       // Adjust point to satisfy bounds and nonlinear inequality constraints
961       if ( hasbounds | nbnlc > 0 ) then
962         [ this , status , xe ] = _scaleinconstraints ( this , xe , xbar )
963         if ( ~status ) then
964           status = "impossibleconstr"
965           break
966         end
967       end
968       if ( nbnlc > 0 ) then
969         // Set the index, so that, if an additionnal cost function argument is provided,
970         // it can be appended at the end.
971         index = 1;
972         [ this , fe ] = neldermead_function ( this , xe , index );
973       else
974         [ this , fe ] = neldermead_function ( this , xe );
975       end
976       this = neldermead_log (this,sprintf("xe=[%s], f(xe)=%f", strcat(string(xe)," ") , fe ));
977       if (fe < fr) then
978         this = neldermead_log (this,sprintf("  > Perform Expansion"));
979         simplex = optimsimplex_setve ( simplex , ihigh , fe , xe )
980       else
981         this = neldermead_log (this,sprintf("  > Perform reflection"));
982         simplex = optimsimplex_setve ( simplex , ihigh , fr , xr )
983       end
984     elseif ( fr >= fn & fr < fhigh ) then
985       // Outside contraction
986       this = neldermead_log (this,sprintf("Contract - outside"));
987       xc = neldermead_interpolate ( xbar , xhigh , this.rho*this.gamma );
988       // Adjust point to satisfy bounds and nonlinear inequality constraints
989       if ( hasbounds | nbnlc > 0 ) then
990         [ this , status , xc ] = _scaleinconstraints ( this , xc , xbar )
991         if ( ~status ) then
992           status = "impossibleconstr"
993           break
994         end
995       end
996       if ( nbnlc > 0 ) then
997         // Set the index, so that, if an additionnal cost function argument is provided,
998         // it can be appended at the end.
999         index = 1;
1000         [ this , fc ] = neldermead_function ( this , xc , index );
1001       else
1002         [ this , fc ] = neldermead_function ( this , xc );
1003       end
1004       this = neldermead_log (this,sprintf("xc=[%s], f(xc)=%f", strcat(string(xc)," ") , fc));
1005       if ( fc <= fr ) then
1006         this = neldermead_log (this,sprintf("  > Perform Outside Contraction"));
1007         simplex = optimsimplex_setve ( simplex , ihigh , fc , xc )
1008       else
1009         //  Shrink
1010         this = neldermead_log (this,sprintf("  > Perform Shrink"));
1011         [ simplex , this ] = optimsimplex_shrink ( simplex , neldermead_costf , this.sigma , this );
1012       end
1013     else
1014       // ( fr >= fn & fr >= fhigh )  
1015       // Inside contraction
1016       this = neldermead_log (this,sprintf("Contract - inside"));
1017       xc = neldermead_interpolate ( xbar , xhigh , -this.gamma );
1018       // Adjust point to satisfy bounds and nonlinear inequality constraints
1019       if ( hasbounds | nbnlc > 0 ) then
1020         [ this , status , xc ] = _scaleinconstraints ( this , xc , xbar )
1021         if ( ~status ) then
1022           status = "impossibleconstr"
1023           break
1024         end
1025       end
1026       if ( nbnlc > 0 ) then
1027         // Set the index, so that, if an additionnal cost function argument is provided,
1028         // it can be appended at the end.
1029         index = 1;
1030         [ this , fc ] = neldermead_function ( this , xc , index );
1031       else
1032         [ this , fc ] = neldermead_function ( this , xc );
1033       end
1034       this = neldermead_log (this,sprintf("xc=[%s], f(xc)=%f", strcat(string(xc)," ") , fc));
1035       if ( fc < fhigh ) then
1036         this = neldermead_log (this,sprintf("  > Perform Inside Contraction"));
1037         simplex = optimsimplex_setve ( simplex , ihigh , fc , xc )
1038       else
1039         //  Shrink
1040         this = neldermead_log (this,sprintf("  > Perform Shrink"));
1041         [ simplex , this ] = optimsimplex_shrink ( simplex , neldermead_costf , this.sigma , this )
1042       end
1043     end
1044     //
1045     // Sort simplex
1046     //
1047     this = neldermead_log (this,sprintf("Sort"));
1048     simplex  = optimsimplex_sort ( simplex );
1049   end
1050   this.optbase = optimbase_set ( this.optbase , "-xopt" , xlow.' );
1051   this.optbase = optimbase_set ( this.optbase , "-fopt" , flow );
1052   this.optbase = optimbase_set ( this.optbase , "-status" , status );
1053   this.simplexopt = simplex;
1054 endfunction
1055