Fixed profile of the cost function to match optim cost functions and features
[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   withderivatives = optimbase_cget ( this.optbase , "-withderivatives" );
41   if ( withderivatives ) then
42     errmsg = msprintf(gettext("%s: The -withderivatives option is true but all algorithms in neldermead are derivative-free."), "neldermead_search")
43     error(errmsg)
44   end
45   if ( ~this.startupflag ) then
46     this = neldermead_startup ( this );
47     this.startupflag = %t;
48   end
49   neldermead_outputcmd ( this, "init" , this.simplex0 , "init" )
50   if ( this.restartflag ) then
51     this = neldermead_autorestart ( this )
52   else
53     this = neldermead_algo ( this );
54   end
55   neldermead_outputcmd ( this, "done" , this.simplexopt , "done" )
56 endfunction
57 //
58 // neldermead_algo --
59 //   Performs an optimization without restart
60 //
61 function this = neldermead_algo ( this )
62     select this.method
63     case "fixed" then
64       this = neldermead_fixed (this);
65     case "variable" then
66       this = neldermead_variable (this);
67     case "box" then
68       this = neldermead_box (this);
69     case "mine" then
70       this = this.mymethod ( this );
71     else
72       errmsg = msprintf(gettext("%s: Unknown -method %s"), "neldermead_algo", this.method)
73       error(errmsg)
74     end
75 endfunction
76 //
77 // neldermead_autorestart --
78 //   Performs an optimization with automatic restart
79 //
80 function this = neldermead_autorestart ( this )
81   for irestart = 1: this.restartmax
82     this = neldermead_log (this,sprintf("Restart #%d/%d", irestart,this.restartmax));
83     this = neldermead_algo ( this );
84     [ this , istorestart ] = neldermead_istorestart ( this );
85     if istorestart==0 then
86       this = neldermead_log (this,"Must not restart");
87       this.restartnb  = irestart
88       break
89     else
90       this = neldermead_log (this,"Must restart");
91     end
92     if ( irestart == this.restartmax ) then
93       this = neldermead_log (this,"Stopping after all restarts performed");
94       this.restartnb  = irestart
95       this.optbase = optimbase_set ( this.optbase , "-status" , "maxrestart" );
96     else
97       this = neldermead_updatesimp ( this );
98     end
99   end
100 endfunction
101
102
103 //
104 // neldermead_variable --
105 //   The original Nelder-Mead algorithm, with variable-size simplex.
106 //
107 function this = neldermead_variable ( this )
108   // Check settings correspond to algo
109   [ this.optbase , hascons ] = optimbase_hasconstraints ( this.optbase );
110   if ( hascons ) then
111       errmsg = msprintf(gettext("%s: Problem has constraints, but variable algorithm ignores them."), "neldermead_variable")
112       error(errmsg)
113   end
114   //
115   // Order the vertices for the first time
116   //
117   simplex = this.simplex0;
118   n = optimbase_cget ( this.optbase , "-numberofvariables" );
119   if (n==0) then
120     errmsg = msprintf(gettext("%s: The number of variable is zero."), "neldermead_variable")
121     error(errmsg)
122   end
123   fvinitial = optimbase_get ( this.optbase , "-fx0" );
124   // Sort function values and x points by increasing function value order
125   this = neldermead_log (this,"Step #1 : order");
126   simplex = optimsimplex_sort ( simplex );
127   currentcenter = optimsimplex_center ( simplex );
128   currentxopt = optimbase_cget ( this.optbase , "-x0" );
129   newfvmean = optimsimplex_fvmean ( simplex );
130   //
131   // Initialize
132   //
133   terminate = %f;
134   iter = 0;
135   step = "init";
136   //
137   // Nelder-Mead Loop
138   //
139   while ( ~terminate )
140     this.optbase = optimbase_incriter ( this.optbase );
141     iter = iter + 1;
142     xlow = optimsimplex_getx ( simplex , 1 )
143     flow = optimsimplex_getfv ( simplex , 1 )
144     xhigh = optimsimplex_getx ( simplex , n+1 )
145     fhigh = optimsimplex_getfv ( simplex , n+1 )
146     xn = optimsimplex_getx ( simplex , n )
147     fn = optimsimplex_getfv ( simplex , n )
148     //
149     // Store history
150     //
151     xcoords = optimsimplex_getallx ( simplex )
152     this = neldermead_storehistory ( this , n , flow , xlow , xcoords );
153     deltafv = abs(optimsimplex_deltafvmax ( simplex ));
154     currentfopt = flow;
155     previousxopt = currentxopt;
156     currentxopt = xlow;
157     previouscenter = currentcenter;
158     currentcenter = optimsimplex_center ( simplex );
159     oldfvmean = newfvmean;
160     newfvmean = optimsimplex_fvmean ( simplex );
161     totaliter = optimbase_get ( this.optbase , "-iterations" );
162     funevals = optimbase_get ( this.optbase , "-funevals" );
163     ssize = optimsimplex_size ( simplex )
164     this = neldermead_log (this,sprintf("================================================================="));
165     this = neldermead_log (this,sprintf("Iteration #%d (total = %d)",iter,totaliter));
166     this = neldermead_log (this,sprintf("Function Eval #%d",funevals));
167     this = neldermead_log (this,sprintf("Xopt : %s",strcat(string(xlow)," ")));
168     this = neldermead_log (this,sprintf("Fopt : %e",flow));
169     this = neldermead_log (this,sprintf("DeltaFv : %e",deltafv));
170     this = neldermead_log (this,sprintf("Center : %s",strcat(string(currentcenter)," ")));
171     this = neldermead_log (this,sprintf("Size : %e",ssize));
172     str = optimsimplex_tostring ( simplex )
173     for i = 1:n+1
174       this = neldermead_log (this,str(i));
175     end
176     this.optbase = optimbase_set ( this.optbase , "-xopt" , xlow );
177     this.optbase = optimbase_set ( this.optbase , "-fopt" , flow );
178     neldermead_outputcmd ( this, "iter" , simplex , step )
179
180     //
181     // Update termination flag
182     //
183     if ( iter > 1 ) then
184       [this , terminate , status] = neldermead_termination (this , ...
185         fvinitial , oldfvmean , newfvmean , previouscenter , currentcenter , simplex );
186       if ( terminate ) then
187         this = neldermead_log (this,sprintf("Terminate with status : %s",status));
188         break
189       end
190     end
191     //
192     // Compute xbar, center of better vertices
193     //
194     this = neldermead_log (this,sprintf("Reflect"));
195     xbar   = optimsimplex_xbar ( simplex );
196     this = neldermead_log (this,sprintf("xbar="+strcat(string(xbar)," ")+""));
197     //
198     // Reflect the worst point with respect to center
199     //
200     xr = neldermead_interpolate ( xbar , xhigh , this.rho );
201     [ this.optbase , fr , index ] = optimbase_function ( this.optbase , xr , 2 );
202     this = neldermead_log (this,sprintf("xr=["+strcat(string(xr)," ")+"], f(xr)=%f",fr));
203     if ( fr >= flow & fr < fn ) then
204       this = neldermead_log (this,sprintf("  > Perform reflection"));
205       simplex = optimsimplex_setve ( simplex , n+1 , fr , xr )
206       step = "reflection";
207     elseif ( fr < flow ) then
208       // Expand
209       this = neldermead_log (this,sprintf("Expand"));
210       xe = neldermead_interpolate ( xbar , xhigh , this.rho*this.chi );
211       [ this.optbase , fe , index ] = optimbase_function ( this.optbase , xe , 2 );
212       this = neldermead_log (this,sprintf("xe="+strcat(string(xe)," ")+", f(xe)=%f",fe));
213       if (fe < fr) then
214         this = neldermead_log (this,sprintf("  > Perform Expansion"));
215         simplex = optimsimplex_setve ( simplex , n+1 , fe , xe )
216         step = "expansion";
217       else
218         this = neldermead_log (this,sprintf("  > Perform reflection"));
219         simplex = optimsimplex_setve ( simplex , n+1 , fr , xr )
220         step = "reflection";
221       end
222     elseif ( fr >= fn & fr < fhigh ) then
223       // Outside contraction
224       this = neldermead_log (this,sprintf("Contract - outside"));
225       xc = neldermead_interpolate ( xbar , xhigh , this.rho*this.gamma );
226       [ this.optbase , fc , index ] = optimbase_function ( this.optbase , xc , 2 );
227       this = neldermead_log (this,sprintf("xc="+strcat(string(xc)," ")+", f(xc)=%f",fc));
228       if ( fc <= fr ) then
229         this = neldermead_log (this,sprintf("  > Perform Outside Contraction"));
230         simplex = optimsimplex_setve ( simplex , n+1 , fc , xc )
231         step = "outsidecontraction";
232       else
233         //  Shrink
234         this = neldermead_log (this,sprintf("  > Perform Shrink"));
235         [ simplex , this ] = optimsimplex_shrink ( simplex , neldermead_costf , this.sigma , this );
236         step = "shrink";
237       end
238     else
239       // ( fr >= fn & fr >= fhigh )  
240       // Inside contraction
241       this = neldermead_log (this,sprintf("Contract - inside"));
242       xc = neldermead_interpolate ( xbar , xhigh , -this.gamma );
243       [ this.optbase , fc , index ] = optimbase_function ( this.optbase , xc , 2 );
244       this = neldermead_log (this,sprintf("xc="+strcat(string(xc)," ")+", f(xc)=%f",fc));
245       if ( fc < fhigh ) then
246         this = neldermead_log (this,sprintf("  > Perform Inside Contraction"));
247         simplex = optimsimplex_setve ( simplex , n+1 , fc , xc )
248         step = "insidecontraction";
249       else
250         //  Shrink
251         this = neldermead_log (this,sprintf("  > Perform Shrink"));
252         [ simplex , this ] = optimsimplex_shrink ( simplex , neldermead_costf , this.sigma , this )
253         step = "shrink";
254       end
255     end
256     //
257     // Sort simplex
258     //
259     this = neldermead_log (this,sprintf("Sort"));
260     simplex  = optimsimplex_sort ( simplex );
261   end
262   this.optbase = optimbase_set ( this.optbase , "-xopt" , xlow.' );
263   this.optbase = optimbase_set ( this.optbase , "-fopt" , flow );
264   this.optbase = optimbase_set ( this.optbase , "-status" , status );
265   this.simplexopt = simplex;
266 endfunction
267
268 //
269 // neldermead_fixed --
270 //   The simplex algorithm with fixed size simplex.
271 // Implementation note:
272 //   We implement the following "rules" of the Spendley et al.
273 //   method.
274 //   Rule 1 is strictly applied, but the reflection is done
275 //   by reflection the high point, since we minimize a function
276 //   instead of maximizing it, like Spendley.
277 //   The Rule 2 is NOT implemented, as we expect that the
278 //   function evaluation is not subject to errors.
279 //   The Rule 3 is applied, ie reflection with respect
280 //   to next to high point.
281 //   A shrink step is included, with shrinkage factor sigma.
282 //
283 //   "Rule 1. Ascertain the lowest reading y, of yi ... Yk+1
284 //   Complete a new simplex Sp by excluding the point Vp corresponding to
285 //   y, and replacing it by V* defined as above."
286 //   "Rule 2. If a result has occurred in (k + 1) successive simplexes, and is not
287 //   then eliminated by application of Rule 1, do not move in the direction
288 //   indicated by Rule 1, or at all, but discard the result and replace it by a new
289 //   observation at the same point."
290 //   "Rule 3. If y is the lowest reading in So , and if the next observation made,
291 //   y* , is the lowest reading in the new simplex S , do not apply Rule 1 and
292 //   return to So from Sp . Move out of S, by rejecting the second lowest reading
293 //   (which is also the second lowest reading in So)."
294 //
295 function this = neldermead_fixed (this)
296   // Check settings correspond to algo
297   [ this.optbase , hascons ] = optimbase_hasnlcons ( this.optbase );
298   if ( hascons ) then
299       errmsg = msprintf(gettext("%s: Problem has constraints, but fixed algorithm ignores them."), "neldermead_fixed")
300       error(errmsg)
301   end
302   //
303   // Order the vertices for the first time
304   //
305   simplex = this.simplex0;
306   n = optimbase_cget ( this.optbase , "-numberofvariables" );
307   fvinitial = optimbase_get ( this.optbase , "-fx0" );
308   // Sort function values and x points by increasing function value order
309   this = neldermead_log (this,sprintf("Sort"));
310   simplex = optimsimplex_sort ( simplex );
311   neldermead_outputcmd ( this, "init" , simplex , "init" )
312   //
313   // Compute center of simplex
314   //
315   currentcenter = optimsimplex_center ( simplex );
316   newfvmean = optimsimplex_fvmean ( simplex );
317   currentxopt = optimbase_cget ( this.optbase , "-x0" );
318   //
319   // Set indices for "clarity"
320   //
321   ilow = 1
322   ihigh = n + 1
323   inext = n
324   //
325   // Initialize
326   //
327   terminate = %f;
328   iter = 0;
329   step = "init";
330   //
331   // main N-M loop
332   //
333   while ( ~terminate )
334     this.optbase = optimbase_incriter ( this.optbase );
335     iter = iter + 1;
336     xlow = optimsimplex_getx ( simplex , ilow )
337     flow = optimsimplex_getfv ( simplex , ilow )
338     xhigh = optimsimplex_getx ( simplex , ihigh )
339     fhigh = optimsimplex_getfv ( simplex , ihigh )
340     //
341     // Store history
342     //
343     xcoords = optimsimplex_getallx ( simplex )
344     this = neldermead_storehistory ( this , n , flow , xlow , xcoords );
345     deltafv = abs(optimsimplex_deltafvmax ( simplex ));
346     currentfopt = flow;
347     previousxopt = currentxopt;
348     currentxopt = xlow;
349     previouscenter = currentcenter;
350     currentcenter = optimsimplex_center ( simplex );
351     oldfvmean = newfvmean;
352     newfvmean = optimsimplex_fvmean ( simplex );
353     totaliter = optimbase_get ( this.optbase , "-iterations" );
354     funevals = optimbase_get ( this.optbase , "-funevals" );
355     ssize = optimsimplex_size ( simplex )
356     this = neldermead_log (this,sprintf("================================================================="));
357     this = neldermead_log (this,sprintf("Iteration #%d (total = %d)",iter,totaliter));
358     this = neldermead_log (this,sprintf("Function Eval #%d",funevals));
359     this = neldermead_log (this,sprintf("Xopt : %s",strcat(string(xlow)," ")));
360     this = neldermead_log (this,sprintf("Fopt : %e",flow));
361     this = neldermead_log (this,sprintf("DeltaFv : %e",deltafv));
362     this = neldermead_log (this,sprintf("Center : %s",strcat(string(currentcenter)," ")));
363     this = neldermead_log (this,sprintf("Size : %e",ssize));
364     str = optimsimplex_tostring ( simplex )
365     for i = 1:n+1
366       this = neldermead_log (this,str(i));
367     end
368     this.optbase = optimbase_set ( this.optbase , "-xopt" , xlow );
369     this.optbase = optimbase_set ( this.optbase , "-fopt" , flow );
370     neldermead_outputcmd ( this, "iter" , simplex , step )
371     //
372     // Update termination flag
373     //
374     if ( iter > 1 ) then
375       [this , terminate , status] = neldermead_termination (this , ...
376         fvinitial , oldfvmean , newfvmean , previouscenter , currentcenter , simplex );
377       if ( terminate ) then
378         this = neldermead_log (this,sprintf("Terminate with status : %s",status));
379         break;
380       end
381     end
382     //
383     // Compute xbar, center of better vertices
384     //
385     this = neldermead_log (this,sprintf("Reflect"));
386     xbar = optimsimplex_xbar ( simplex );
387     this = neldermead_log (this,sprintf("xbar="+strcat(string(xbar)," ")+""));
388     //
389     // Reflect the worst point with respect to center
390     //
391     xr = neldermead_interpolate ( xbar , xhigh , this.rho );
392     [ this.optbase , fr , index ] = optimbase_function ( this.optbase , xr , 2 );
393     this = neldermead_log (this,sprintf("xr="+strcat(string(xr)," ")+", f(xr)=%f",fr));
394     //
395     // Replace worst point by xr if it is better
396     //
397     if ( fr < fhigh ) then
398       this = neldermead_log (this,sprintf("  > Perform reflect"));
399       simplex = optimsimplex_setve ( simplex , ihigh , fr , xr )
400       step = "reflection";
401     else
402       // Reflect / xnext
403       xnext = optimsimplex_getx ( simplex , inext );
404       fnext = optimsimplex_getfv ( simplex , inext );
405       xbar2 = optimsimplex_xbar ( simplex , inext );
406       this = neldermead_log (this,sprintf("xbar2="+strcat(string(xbar2)," ")+""));
407       xr2 = neldermead_interpolate ( xbar2 , xnext , this.rho );
408       [ this.optbase , fr2 , index ] = optimbase_function ( this.optbase , xr2 , 2 );
409       this = neldermead_log (this,sprintf("xr2="+strcat(string(xr2)," ")+", f(xr2)=%f",fr2));
410       if ( fr2 < fnext ) then
411         this = neldermead_log (this,sprintf("  > Perform reflect / next"));
412         simplex = optimsimplex_setve ( simplex , inext , fr2 , xr2 )
413         step = "reflectionnext";
414       else
415         //  Shrink
416         this = neldermead_log (this,sprintf("  > Perform Shrink"));
417         [ simplex , this ] = optimsimplex_shrink ( simplex , neldermead_costf , this.sigma , this )
418         step = "shrink";
419       end
420     end
421     //
422     // Sort simplex
423     //
424     simplex = optimsimplex_sort ( simplex );
425   end
426   this.optbase = optimbase_set ( this.optbase , "-xopt" , xlow.' );
427   this.optbase = optimbase_set ( this.optbase , "-fopt" , flow );
428   this.optbase = optimbase_set ( this.optbase , "-status" , status );
429   this.simplexopt = simplex;
430 endfunction
431 //
432 // neldermead_interpolate --
433 //   Computes xi as an interpolation between x1 and x2, with factor as :
434 //     xi = (1+fac)x1 - fac x2
435 //
436 function xi = neldermead_interpolate ( x1 , x2 , fac )
437   xi = (1 + fac)*x1 - fac*x2;
438 endfunction
439
440 //
441 // neldermead_termination --
442 //   Returns 1 if the algorithm terminates.
443 //   Returns 0 if the algorithm must continue.
444 // Arguments
445 //   this : the current object
446 //   fvinitial : initial function value
447 //   newfvmean, oldfvmean : the old and new function value average on the simplex
448 //   previousxopt : the previous value of x
449 //   currentxopt : the current value of x
450 //   simplex : the simplex
451 //     The best point in the simplex is expected to be stored at 1
452 //     The worst point in the simplex is expected to be stored at n+1
453 //   terminate : %t if the algorithm terminates, %f if the algorithm must continue.
454 //   this.status : termination status
455 //     "continue"
456 //     "maxiter"
457 //     "maxfuneval"
458 //     "tolf"
459 //     "tolx"
460 //     "tolfstdev"
461 //     "tolsize"
462 //     "tolsizedeltafv"
463 //     "kelleystagnation"
464 //     "tolboxf"
465 //     "tolvariance"
466 // Notes
467 //   Use the function average on the simplex instead of the best function value.
468 //   This is because the function average changes at each iteration.
469 //   Instead, the best function value as a step-by-step evolution and may not
470 //   change in 2 iterations, leading to astop of the algorithm.
471 // TODO : set the fvinitial, oldfvmean, newfvmean.
472 //
473 function [this , terminate , status ] = neldermead_termination (this , ...
474   fvinitial , oldfvmean , newfvmean , previousxopt , currentxopt , ...
475   simplex )
476   terminate = %f;
477   status = "continue";
478   //
479   // Termination Criteria from parent optimization class
480   //
481   [ this.optbase , terminate , status] = optimbase_terminate ( this.optbase , ...
482     fvinitial , newfvmean , previousxopt , currentxopt );
483   //
484   // Criteria #5 : standard deviation of function values
485   //
486   if ( ~terminate ) then
487     if this.tolfstdeviationmethod == "enabled" then
488       fv = optimsimplex_getallfv ( simplex )
489       sd = st_deviation(fv);
490       this.optbase = optimbase_stoplog  ( this.optbase,sprintf("  > st_deviation(fv)=%e < tolfstdeviation=%e",...
491         sd, this.tolfstdeviation));
492       if sd < this.tolfstdeviation then
493         terminate = %t;
494         status = "tolfstdev";
495       end
496     end
497   end
498   //
499   // Criteria #6 : simplex absolute + relative size
500   //
501   if ( ~terminate ) then
502     if ( this.tolsimplexizemethod ) then
503       ssize = optimsimplex_size ( simplex , "sigmaplus" );
504       this.optbase = optimbase_stoplog  ( this.optbase,sprintf("  > simplex size=%e < %e + %e * %e",...
505         ssize, this.tolsimplexizeabsolute , this.tolsimplexizerelative , this.simplexsize0 ));
506       if ssize < this.tolsimplexizeabsolute + this.tolsimplexizerelative * this.simplexsize0 then
507         terminate = %t;
508         status = "tolsize";
509       end
510     end
511   end
512   //
513   // Criteria #7 : simplex absolute size + difference in function values (Matlab-like)
514   //
515   if ( ~terminate ) then
516     if ( this.tolssizedeltafvmethod ) then
517       ssize = optimsimplex_size ( simplex , "sigmaplus" );
518       this.optbase = optimbase_stoplog  ( this.optbase,sprintf("  > simplex size=%e < %e",...
519         ssize, this.tolsimplexizeabsolute));
520       shiftfv = abs(optimsimplex_deltafvmax( simplex ))
521       this.optbase = optimbase_stoplog  ( this.optbase,sprintf("  > abs(fv(n+1) - fv(1))=%e < toldeltafv=%e",...
522         shiftfv, this.toldeltafv));
523       if ( ( ssize < this.tolsimplexizeabsolute ) & ( shiftfv < this.toldeltafv ) ) then
524         terminate = %t;
525         status = "tolsizedeltafv";
526       end
527     end
528   end
529   //
530   // Criteria #8 : Kelley stagnation, based on
531   // a sufficient decrease condition
532   //
533   if ( ~terminate ) then
534     if ( this.kelleystagnationflag ) then
535       [ sg , this ] = optimsimplex_gradientfv ( simplex , neldermead_costf , "forward" , this );
536       nsg = sg.' * sg;
537       sgstr = strcat(string(sg)," ");
538       this.optbase = optimbase_stoplog ( this.optbase , sprintf ( "Test Stagnation : nsg = %e, sg = "+sgstr, nsg) );
539       this.optbase = optimbase_stoplog ( this.optbase , ...
540         sprintf ( "Test Stagnation : newfvmean=%e >= oldfvmean=%e - %e * %e" , newfvmean, oldfvmean , this.kelleyalpha , nsg ) );
541       if ( newfvmean >= oldfvmean - this.kelleyalpha * nsg ) then
542         terminate = %t;
543         status = "kelleystagnation";
544       end
545     end
546   end
547   //
548   // Criteria #9 : Box termination criteria
549   // The number of consecutive time that an absolute tolerance on
550   // function value is met.
551   // From Algorithm 454, the tolerance is the difference between the
552   // max and the min function values in the simplex.
553   //
554   if ( ~terminate ) then
555     if ( this.boxtermination ) then
556       shiftfv = abs(optimsimplex_deltafvmax( simplex ))
557       this.optbase = optimbase_stoplog ( this.optbase , ...
558         sprintf ( "Test Box : shiftfv=%e < boxtolf=%e" , shiftfv , this.boxtolf ) );
559       if ( shiftfv < this.boxtolf ) then
560         this.boxkount = this.boxkount + 1
561       this.optbase = optimbase_stoplog ( this.optbase , ...
562         sprintf ( "Test Box : boxkount=%d == boxnbmatch=%d" , this.boxkount , this.boxnbmatch ) );
563         if ( this.boxkount == this.boxnbmatch ) then
564           terminate = %t
565           status = "tolboxf"
566         end
567       else
568         this.boxkount = 0
569       end
570     end
571   end
572   //
573   // Criteria #10 : variance of function values
574   //
575   if ( ~terminate ) then
576     if ( this.tolvarianceflag ) then
577       var = optimsimplex_fvvariance ( simplex )
578       this.optbase = optimbase_stoplog ( this.optbase , ...
579         sprintf ( "Test tolvariance : %e < %e" , var , this.tolabsolutevariance ) );
580       if ( var < this.tolrelativevariance * this.variancesimplex0 + this.tolabsolutevariance ) then
581         terminate = %t
582         status = "tolvariance"
583       end
584     end
585   end
586   //
587   // Criteria #11 : user-defined criteria
588   //
589   if ( ~terminate ) then
590     if ( this.myterminateflag ) then
591       [ this , term , stat ] = this.myterminate ( this , simplex )
592       if ( term ) then 
593         terminate = term
594         status = stat
595       end
596     end
597   end
598   this.optbase = optimbase_stoplog (this.optbase,sprintf("  > Terminate = %s, status = %s",...
599     string(terminate) , status ));
600 endfunction
601   
602
603 //
604 // neldermead_outputcmd --
605 //   Calls back user's output command
606 // Arguments
607 //   this : the current object
608 //   state : the state of the algorithm,
609 //     "init", "done", "iter"
610 //   simplex : the current simplex
611 //   step : the type of step performed during the iteration
612 //     "init", "done", "reflection", "expansion", "insidecontraction", "outsidecontraction"
613 //     "reflectionnext", "shrink"
614 //
615 function  neldermead_outputcmd ( this, ...
616    state , simplex , step )
617   outputcmd = optimbase_cget ( this.optbase , "-outputcommand" );
618   if typeof(outputcmd) <> "string" then
619     brutedata = optimbase_outstruct ( this.optbase );
620     data = tlist(["T_NMDATA",...
621       "x","fval","iteration","funccount",...
622       "simplex" , "step" ]);
623     data.x = brutedata.x;
624     data.fval = brutedata.fval;
625     data.iteration = brutedata.iteration;
626     data.funccount = brutedata.funccount;
627     data.simplex = simplex;
628     data.step = step;
629     optimbase_outputcmd ( this.optbase , state , data );
630   end
631 endfunction
632 //
633 // neldermead_storehistory --
634 //   Stores the history into the data structure.
635 // Arguments, input
636 //   this : current object
637 //   n : the number of unknown parameters
638 //   fopt : the current value of the function at optimum
639 //   xopt : arrary with size n, current optimum
640 //   xcoords : array with size n x n+1, coordinates of the n+1 vertices
641 //
642 function this = neldermead_storehistory ( this , n , fopt , xopt , xcoords )
643   storehistory = optimbase_cget ( this.optbase , "-storehistory" );
644   iterations = optimbase_get ( this.optbase , "-iterations" );
645   if ( storehistory ) then
646     this.optbase = optimbase_histset ( this.optbase , iterations , "-fopt" , fopt );
647     this.optbase = optimbase_histset ( this.optbase , iterations , "-xopt" , xopt(1:n).' );
648     this.historysimplex ( iterations , 1:n+1,1:n) = xcoords(1:n+1,1:n);
649   end
650 endfunction
651
652 //
653 // neldermead_istorestart --
654 //   Returns 1 if the optimization is to restart.
655 // Arguments
656 //   istorestart : 1 of the the optimization is to restart.
657 //
658 function [ this , istorestart ] = neldermead_istorestart ( this )
659   select this.restartdetection
660   case "oneill"
661     [ this , istorestart ] = neldermead_isroneill ( this )
662   case "kelley"
663     [ this , istorestart ] = neldermead_isrkelley ( this )
664   else
665     errmsg = msprintf(gettext("%s: Unknown restart detection %s"),"neldermead_istorestart", this.restartdetection)
666     error(errmsg)
667   end
668 endfunction
669 //
670 // neldermead_isrkelley--
671 //   Returns 1 if the optimization is to restart.
672 //   Use kelleystagnation as a criteria for restart.
673 // Arguments
674 //   istorestart : 1 of the the optimization is to restart.
675 //
676 function [ this , istorestart ] = neldermead_isrkelley ( this )
677   istorestart = 0
678   if ( this.kelleystagnationflag ) then
679     status = optimbase_get ( this.optbase , "-status" );
680     if ( status =="kelleystagnation" ) then
681        istorestart = 1
682     end
683   end
684 endfunction
685 //
686 // neldermead_isroneill --
687 //   Returns 1 if the optimization is to restart.
688 //   Use O'Neill method as a criteria for restart.
689 //   It is a axis by axis search for optimality.
690 // Arguments
691 //   xopt : the optimum, as a st of n values
692 //   fopt : function value at optimum
693 //   eps : a small value
694 //   step : a list of n values, representing
695 //     the "size" of each parameter
696 //   istorestart : 1 of the the optimization is to restart.
697 //
698 function [ this , istorestart ] = neldermead_isroneill ( this )
699   n = optimbase_cget ( this.optbase , "-numberofvariables" );
700   //
701   // If required, make a vector step from the scalar step
702   //
703   defaultstep = this.restartstep;
704   stepn = length ( defaultstep );
705   if ( stepn <> n ) then
706     step = defaultstep * ones(n,1);
707   else
708     step = defaultstep;
709   end
710
711   xopt = optimbase_get ( this.optbase , "-xopt" );
712   fopt = optimbase_get ( this.optbase , "-fopt" );
713
714     istorestart = 0
715     for ix = 1:n
716       stepix = step ( ix )
717       del = stepix * this.restarteps
718       if ( del==0.0 ) then
719          del = eps
720       end
721       xoptix =  xopt ( ix )
722       xopt ( ix ) = xoptix + del
723       [ this.optbase , fv , index ] = optimbase_function ( this.optbase , xopt , 2 )
724       if ( fv < fopt ) then
725         istorestart = 1
726         break
727       end
728       xopt ( ix ) = xoptix - del
729       [ this.optbase , fv , index ] = optimbase_function ( this.optbase , xopt , 2 )
730       if ( fv < fopt ) then
731         istorestart = 1
732         break
733       end
734       xopt ( ix ) = xoptix
735     end
736 endfunction
737
738 //
739 // neldermead_startup --
740 //   Startup the algorithm.
741 //   Computes the initial simplex, depending on the -simplex0method.
742 //
743 function this = neldermead_startup (this)
744   // 0. Check that the cost function is correctly connected
745   // Note: this call to the cost function is not used, but helps the
746   // user while he is tuning his object.
747   if ( this.checkcostfunction ) then
748     this.optbase = optimbase_checkcostfun ( this.optbase );
749   end
750   // 1. If the problem has bounds, check that they are consistent
751   [ this.optbase , hasbounds ] = optimbase_hasbounds ( this.optbase );
752   if ( hasbounds ) then
753     [ this.optbase , isok , errmsg ] = optimbase_checkbounds ( this.optbase );
754     if ( ~isok ) then
755       error ( msprintf(gettext("%s: %s"), "neldermead_startup" , errmsg ))
756     end
757   end
758   // 2. Get the initial guess and compute the initial simplex
759   x0 = optimbase_cget ( this.optbase , "-x0" );
760   select this.simplex0method
761   case "given" then
762     [ simplex0 , this ] = optimsimplex_new ( this.coords0 , ...
763       neldermead_costf , this );
764   case "axes" then
765     [ simplex0 , this ] = optimsimplex_new ( "axes" , ...
766       x0.' , neldermead_costf , this.simplex0length , this );
767   case "spendley" then
768     [ simplex0 , this ] = optimsimplex_new ( "spendley" , ...
769       x0.' , neldermead_costf , this.simplex0length , this );
770   case "pfeffer" then
771     [ simplex0 , this ] = optimsimplex_new ( "pfeffer" , ...
772       x0.' , neldermead_costf , this.simplex0deltausual , ...
773       this.simplex0deltazero , this );
774   case "randbounds" then
775     if ( this.boxnbpoints == "2n" ) then
776       this.boxnbpointseff = 2 * this.optbase.numberofvariables;
777     else
778       this.boxnbpointseff = this.boxnbpoints;
779     end
780     if ( ~hasbounds ) then
781       error ( msprintf(gettext("%s: Randomized bounds initial simplex is not available without bounds." ), "neldermead_startup"))
782     end
783     [ simplex0 , this ] = optimsimplex_new ( "randbounds" , x0.' , ...
784       neldermead_costf , this.optbase.boundsmin , this.optbase.boundsmax , ...
785       this.boxnbpointseff  , this );
786   else
787     errmsg = msprintf(gettext("%s: Unknown -simplex0method : %s"), "neldermead_startup", this.simplex0method);
788     error(errmsg);
789   end
790   //
791   // 3. Scale the simplex into the bounds and the nonlinear inequality constraints, if any
792   //
793   if ( hasbounds | this.optbase.nbineqconst > 0 ) then
794     this = neldermead_log (this,sprintf("Scaling initial simplex into nonlinear inequality constraints..."));
795     select this.scalingmethod
796     case "tox0" then
797       [ this , simplex0 ] = neldermead_scaletox0 ( this , simplex0 );
798     case "tocenter" then
799       [ this , simplex0 ] = neldermead_scaletocenter ( this , simplex0 );
800     else
801       errmsg = msprintf(gettext("%s: Unknown value %s for -scalingmethod option"),"neldermead_startup", this.scalingmethod );
802       error(errmsg);
803     end
804   end
805   //
806   // 4. Store the simplex
807   //
808   this.simplex0 = optimsimplex_destroy ( this.simplex0 );
809   this.simplex0 = simplex0;
810   this.simplexsize0 = optimsimplex_size ( simplex0 );
811   // 5. Store initial data into the base optimization component
812   fx0 = optimsimplex_getfv ( this.simplex0 , 1 );
813   this.optbase = optimbase_set ( this.optbase , "-fx0" , fx0 );
814   this.optbase = optimbase_set ( this.optbase , "-xopt" , x0.' );
815   this.optbase = optimbase_set ( this.optbase , "-fopt" , fx0 );
816   this.optbase = optimbase_set ( this.optbase , "-iterations" , 0 );
817   // 6. Initialize the termination criteria
818   this = neldermead_termstartup ( this );
819 endfunction
820 //
821 // neldermead_scaletox0 --
822 //   Scale the simplex into the 
823 //   nonlinear inequality constraints, if any.
824 //   Scale toward x0, which is feasible.
825 // Arguments
826 //   simplex0 : the initial simplex
827 //
828 function [ this , simplex0 ] = neldermead_scaletox0 ( this , simplex0 )
829     nbve = optimsimplex_getnbve ( simplex0 );
830     x0 = optimbase_cget ( this.optbase , "-x0" );
831     for ive = 2 : nbve
832       // optimsimplex returns a row vector
833       x = optimsimplex_getx ( simplex0 , ive );
834       this = neldermead_log (this,sprintf("Scaling vertex #%d/%d at ["+...
835         strcat(string(x)," ")+"]... " , ...
836         ive , nbve ));
837       // Transpose x into a row vector
838       [ this , status , xp ] = _scaleinconstraints ( this , x.' , x0 );
839       if ( ~status ) then
840         errmsg = msprintf(gettext("%s: Impossible to scale the vertex #%d/%d at [%s] into inequality constraints"), ...
841           "neldermead_startup", ive , nbve , strcat(string(x)," "));
842         error(errmsg);
843       end
844       if ( or ( x <> xp ) ) then
845         [ this.optbase , fv , c , index ] = optimbase_function ( this.optbase , xp , 2 );
846         // Transpose xp, which is a column vector
847         simplex0 = optimsimplex_setve ( simplex0 , ive , fv , xp.' );
848       end
849     end
850 endfunction
851 //
852 // neldermead_scaletocenter --
853 //   Scale the simplex into the 
854 //   nonlinear inequality constraints, if any.
855 //   Scale to the centroid of the points
856 //   which satisfy the constraints.
857 // Notes
858 //   This is Box's method for scaling.
859 //   It is unsure, since the centroid of the points
860 //   which satisfy the constraints may not be feasible.
861 // TODO : test !
862 // TODO : insert a note for that specific point
863 // Arguments
864 //   
865 //
866 function [ this , simplex0 ] = neldermead_scaletocenter ( this , simplex0 , x0 )
867     nbve = optimsimplex_getnbve ( simplex0 );
868     xref = optimsimplex_getx ( simplex0 , 1 );
869     for ive = 2 : nbve
870       xref = optimsimplex_xbar ( simplex0 , ive:nbve );
871       // optimsimplex returns a row vector
872       x = optimsimplex_getx ( simplex0 , ive );
873       this = neldermead_log (this,sprintf("Scaling vertex #%d/%d at ["+...
874         strcat(string(x)," ")+"]... " , ...
875         ive , nbve ));
876       // Transpose x into a row vector
877       [ this , status , xp ] = _scaleinconstraints ( this , x.' , xref );
878       if ( ~status ) then
879         errmsg = msprintf(gettext("%s: Impossible to scale the vertex #%d/%d at [%s] into inequality constraints"), ...
880           "neldermead_startup", ive , nbve , strcat(string(x)," "));
881         error(errmsg);
882       end
883       if ( or ( x <> xp ) ) then
884         [ this.optbase , fv , c , index ] = optimbase_function ( this.optbase , xp , 2 );
885         // Transpose xp, which is a column vector
886         simplex0 = optimsimplex_setve ( simplex0 , ive , fv , xp.' );
887       end
888     end
889 endfunction
890 //
891 // neldermead_termstartup --
892 //   Initialize Kelley's stagnation detection system when normalization is required,
893 //   by computing kelleyalpha.
894 //   If the simplex gradient is zero, then
895 //   use alpha0 as alpha.
896 // Arguments
897 //   status : the status after the failing
898 //     optimization process
899 //   simplex : the simplex computed at the end of the failing
900 //     optimization process
901 //
902 function this = neldermead_termstartup ( this )
903   //
904   // Criteria #8 : Kelley stagnation, based on
905   // a sufficient decrease condition
906   //
907   if ( this.kelleystagnationflag ) then
908       if ( ~this.kelleynormalizationflag ) then
909         this.kelleyalpha = this.kelleystagnationalpha0
910       else
911         [sg,this] = optimsimplex_gradientfv ( this.simplex0 , neldermead_costf , "forward" , this )
912         nsg = sg.' * sg
913         sigma0 = optimsimplex_size ( this.simplex0 , "sigmaplus" )
914         if nsg==0.0 then 
915           this.kelleyalpha = this.kelleystagnationalpha0
916         else
917           this.kelleyalpha = this.kelleystagnationalpha0 * sigma0 / nsg
918         end
919       end
920       this = neldermead_log (this,sprintf("Test Stagnation Kelley : alpha0 = %e", this.kelleyalpha));
921   end
922   //
923   // Criteria #10 : variance of function values
924   //
925   if ( this.tolvarianceflag ) then
926       this.variancesimplex0 = optimsimplex_fvvariance ( this.simplex0 )
927   end
928 endfunction
929 //
930 // _scaleinconstraints --
931 //   Given a point to scale and a reference point which satisfies the constraints,
932 //   scale the point towards the reference point until it satisfies all the constraints.
933 //   Returns isscaled = %T if the procedure has succeded before -boxnbnlloops
934 //   Returns isscaled = %F if the procedure has failed after -boxnbnlloops
935 //   iterations.
936 // Arguments
937 //   x : the point to scale
938 //   xref : the reference point
939 //   isscaled : %T or %F
940 //   p : scaled point
941 //
942 function [ this , isscaled , p ] = _scaleinconstraints ( this , x , xref )
943   p = x
944   //
945   // Adjust point to satisfy nonlinear inequality constraints
946   //
947   nbnlc = optimbase_cget ( this.optbase , "-nbineqconst" )
948   if ( nbnlc == 0 ) then
949     isscaled = %T
950     return;
951   end
952   isscaled = %F
953   //
954   // Try the current point and see if the constraints are satisfied.
955   // If not, move the point "halfway" to the centroid,
956   // which should satisfy the constraints, if
957   // the constraints are convex.
958   // Perform this loop until the constraints are satisfied.
959   // If all loops have been performed without success, the scaling
960   // has failed.
961   //
962   alpha = 1.0
963   p0 = p
964   while ( alpha > this.guinalphamin )
965       [ this.optbase , feasible ] = optimbase_isinnonlincons ( this.optbase , p );
966       if ( feasible ) then
967         isscaled = %T;
968         break;
969       end
970       alpha = alpha / 2.0
971       this = neldermead_log (this,sprintf("Scaling inequality constraint with alpha = %e", ...
972         alpha));
973       p = ( 1.0 - alpha ) * xref + alpha * p0;
974   end
975   this = neldermead_log (this,sprintf(" > After scaling into inequality constraints p = [%s]" , ...
976     _strvec(p) ) );
977   if ( ~isscaled ) then
978     this = neldermead_log (this,sprintf(" > Impossible to scale into constraints after %d loops" , ...
979       this.optbase.nbineqconst ));
980   end
981 endfunction
982 //
983 // neldermead_box --
984 //   The Nelder-Mead algorithm, with variable-size simplex
985 //   and modifications by Box for bounds and
986 //   inequality constraints.
987 //
988 function this = neldermead_box ( this )
989   // Check settings correspond to algo
990   [ this.optbase , hascons ] = optimbase_hasconstraints ( this.optbase );
991   if ( ~hascons ) then
992       errmsg = msprintf(gettext("%s: Problem has no constraints, but Box algorithm is designed for them."), "neldermead_box")
993       error(errmsg)
994   end
995   //
996   // Order the vertices for the first time
997   //
998   simplex = this.simplex0;
999   n = optimbase_cget ( this.optbase , "-numberofvariables" );
1000   fvinitial = optimbase_get ( this.optbase , "-fx0" );
1001   // Sort function values and x points by increasing function value order
1002   this = neldermead_log (this,"Step #1 : order");
1003   simplex = optimsimplex_sort ( simplex );
1004   currentcenter = optimsimplex_center ( simplex );
1005   currentxopt = optimbase_cget ( this.optbase , "-x0" );
1006   newfvmean = optimsimplex_fvmean ( simplex );
1007   nbve = optimsimplex_getnbve ( simplex );
1008   ihigh = nbve;
1009   inext = ihigh - 1
1010   ilow = 1
1011   [ this.optbase , hasbounds ] = optimbase_hasbounds ( this.optbase );
1012   nbnlc = optimbase_cget ( this.optbase , "-nbineqconst" )
1013   rho = this.boxreflect;
1014   //
1015   // Initialize
1016   //
1017   terminate = %f;
1018   iter = 0;
1019   step = "init";
1020   //
1021   // Nelder-Mead Loop
1022   //
1023   while ( ~terminate )
1024     this.optbase = optimbase_incriter ( this.optbase );
1025     iter = iter + 1;
1026     xlow = optimsimplex_getx ( simplex , ilow )
1027     flow = optimsimplex_getfv ( simplex , ilow )
1028     xhigh = optimsimplex_getx ( simplex , ihigh )
1029     fhigh = optimsimplex_getfv ( simplex , ihigh )
1030     xn = optimsimplex_getx ( simplex , inext )
1031     fn = optimsimplex_getfv ( simplex , inext )
1032     //
1033     // Store history
1034     //
1035     xcoords = optimsimplex_getallx ( simplex )
1036     this = neldermead_storehistory ( this , n , flow , xlow , xcoords );
1037     deltafv = abs(optimsimplex_deltafvmax ( simplex ));
1038     currentfopt = flow;
1039     previousxopt = currentxopt;
1040     currentxopt = xlow;
1041     previouscenter = currentcenter;
1042     currentcenter = optimsimplex_center ( simplex );
1043     oldfvmean = newfvmean;
1044     newfvmean = optimsimplex_fvmean ( simplex );
1045     totaliter = optimbase_get ( this.optbase , "-iterations" );
1046     funevals = optimbase_get ( this.optbase , "-funevals" );
1047     ssize = optimsimplex_size ( simplex )
1048     this = neldermead_log (this,sprintf("================================================================="));
1049     this = neldermead_log (this,sprintf("Iteration #%d (total = %d)",iter,totaliter));
1050     this = neldermead_log (this,sprintf("Function Eval #%d",funevals));
1051     this = neldermead_log (this,sprintf("Xopt : [%s]",_strvec(xlow)));
1052     this = neldermead_log (this,sprintf("Fopt : %e",flow));
1053     this = neldermead_log (this,sprintf("DeltaFv : %e",deltafv));
1054     this = neldermead_log (this,sprintf("Center : [%s]",_strvec(currentcenter)));
1055     this = neldermead_log (this,sprintf("Size : %e",ssize));
1056     str = optimsimplex_tostring ( simplex )
1057     for i = 1:nbve
1058       this = neldermead_log (this,str(i));
1059     end
1060     neldermead_outputcmd ( this, "iter" , simplex , step )
1061
1062     //
1063     // Update termination flag
1064     //
1065     if ( iter > 1 ) then
1066       [this , terminate , status] = neldermead_termination (this , ...
1067         fvinitial , oldfvmean , newfvmean , previouscenter , currentcenter , simplex );
1068       if ( terminate ) then
1069         this = neldermead_log (this,sprintf("Terminate with status : %s",status));
1070         break
1071       end
1072     end
1073     //
1074     // Compute xbar, center of better vertices
1075     //
1076     this = neldermead_log (this,sprintf("Reflect"));
1077     xbar = optimsimplex_xbar ( simplex );
1078     this = neldermead_log (this,sprintf("xbar=[%s]" , _strvec(xbar)));
1079     //
1080     // Search a point improving cost function
1081     // and satisfying constraints.
1082     //
1083     [ this , scaled , xr , fr ] = _boxlinesearch ( this , n , xbar , xhigh , fhigh , rho );
1084     if ( scaled == %f ) then
1085       status = "impossibleimprovement"
1086       break
1087     end
1088     this = neldermead_log (this,sprintf("xr=[%s], f(xr)=%f", strcat(string(xr)," ") , fr));
1089     this = neldermead_log (this,sprintf("  > Perform Reflection"));
1090     simplex = optimsimplex_setve ( simplex , ihigh , fr , xr )
1091     step = "boxreflection";
1092     //
1093     // Sort simplex
1094     //
1095     this = neldermead_log (this,sprintf("Sort"));
1096     simplex  = optimsimplex_sort ( simplex );
1097   end
1098   this.optbase = optimbase_set ( this.optbase , "-xopt" , xlow.' );
1099   this.optbase = optimbase_set ( this.optbase , "-fopt" , flow );
1100   this.optbase = optimbase_set ( this.optbase , "-status" , status );
1101   this.simplexopt = simplex;
1102 endfunction
1103
1104   //
1105   // _strvec --
1106   //  Returns a string for the given vector.
1107   //
1108   function str = _strvec ( x )
1109     str = strcat(string(x)," ")
1110   endfunction
1111   //
1112   // _boxlinesearch --
1113   //   For Box's method, perform a line search
1114   //   from xbar, on the line (xhigh,xbar) and returns:
1115   //   status : %t if the search is successful
1116   //   xr : the reflected point
1117   //   fr : the function value
1118   //   The reflected point satisfies the following
1119   //   constraints :
1120   //   * fr < fhigh
1121   //   * xr satisfies the bounds constraints
1122   //   * xr satisfies the nonlinear positive inequality constraints
1123   //   * xr satisfies the linear positive inequality constraints
1124   //   The method is based on projection and
1125   //   scaling toward the centroid.
1126   //
1127   // Arguments
1128   //   n : number of variables
1129   //   xbar : the centroid
1130   //   xhigh : the high point
1131   //   fhigh : function value at xhigh
1132   //   rho : reflection factor
1133   //
1134   function [ this , status , xr , fr ] = _boxlinesearch ( this , n , xbar , xhigh , fhigh , rho )
1135     this = neldermead_log (this,"_boxlinesearch");
1136     this = neldermead_log (this, sprintf ("> xhigh=[%s], fhigh=%e",_strvec(xhigh),fhigh));
1137     this = neldermead_log (this, sprintf ( "> xbar=[%s]" , _strvec(xbar) ) );
1138     xr = neldermead_interpolate ( xbar , xhigh , rho );
1139     this = neldermead_log (this, sprintf ( "> xr = [%s]" , _strvec ( xr ) ) );
1140     status = %f
1141     alphamin = this.guinalphamin
1142     //
1143     // Scale from xr toward xbar until fr < fhigh and update xr
1144     //
1145     xr0 = xr
1146     alpha = 1.0
1147     while ( alpha > alphamin )
1148       [ this.optbase , fr , cr , index ] = optimbase_function ( this.optbase , xr , 2 );
1149       if ( fr < fhigh ) then
1150         this = neldermead_log (this, sprintf ( "fr = %e improves %e : no need for scaling for f" , fr , fhigh ) );
1151         status = %t;
1152         break
1153       end
1154       alpha = alpha / 2.0;
1155       this = neldermead_log (this, sprintf ( "Scaling for f with alpha=%e" , alpha ) );
1156       xr = ( 1.0 - alpha ) * xbar + alpha * xr0;
1157       this = neldermead_log (this, sprintf ( "> xr = %s" , _strvec ( xr ) ) );
1158     end
1159     // If the scaling for function improvement has failed,
1160     // we return.
1161     if ( ~status ) then
1162       return;
1163     end
1164     // scaledc is set to %t if xr is updated during scaling into constraints 
1165     // That implies that the function value is to update.
1166     scaledc = %f
1167     //
1168     // Project xr into bounds, with an additionnal alpha inside the bounds.
1169     // This algo is always succesful.
1170     // Note:
1171     //   If the alpha coefficient was not used, the
1172     //   projectinbounds method could be used directly.
1173     //
1174     [ this.optbase , hasbounds ] = optimbase_hasbounds ( this.optbase );
1175     if ( hasbounds ) then
1176       boxboundsalpha = this.boxboundsalpha;
1177       boundsmax = optimbase_cget ( this.optbase , "-boundsmax" );
1178       boundsmin = optimbase_cget ( this.optbase , "-boundsmin" );
1179       for ix = 1:n
1180         xmin = boundsmin ( ix );
1181         xmax = boundsmax ( ix );
1182         xrix = xr ( ix );
1183         if ( xrix > xmax ) then
1184           this = neldermead_log (this, sprintf ( "Projecting index #%d = %e on max bound %e - %e" , ix , xrix , xmax , boxboundsalpha ) );
1185           xr ( ix ) = xmax - boxboundsalpha;
1186           if ( ~scaledc ) then
1187             scaledc = %t
1188           end
1189         elseif ( xrix < xmin ) then
1190           this = neldermead_log (this, sprintf ( "Projecting index #%e = %e on min bound %e - %e" , ix , xrix , xmin , boxboundsalpha ) );
1191           xr ( ix ) = xmin + boxboundsalpha;
1192           if ( ~scaledc ) then
1193             scaledc = %t
1194           end
1195         end
1196       end
1197       this = neldermead_log (this, sprintf ( " > After projection into bounds xr = [%s]" , _strvec(xr)));
1198     end
1199     //
1200     // Scale from xr to xbar into nonlinear inequality constraints
1201     // and update xr. 
1202     // Set status to 0 if the process fails.
1203     //
1204     nbnlc = optimbase_cget ( this.optbase , "-nbineqconst" );
1205     if ( nbnlc == 0 ) then
1206       status = %t
1207     else
1208       status = %f;
1209       alpha = 1.0;
1210       xr0 = xr;
1211       while ( alpha > alphamin )
1212         [ this.optbase , feasible ] = optimbase_isinnonlincons ( this.optbase , xr );
1213         if ( feasible ) then
1214           status = %t;
1215           break
1216         end
1217         alpha = alpha / 2.0;
1218         this = neldermead_log (this, sprintf ( "Scaling for nonlinear/linear inequality constraints with alpha=%e from xbar=[%s] toward [%s]" , ...
1219           alpha , _strvec(xbar) , _strvec(xr0) ));
1220         xr = ( 1.0 - alpha ) * xbar + alpha * xr0;
1221         this = neldermead_log (this, sprintf ( "> xr = [%s]" , _strvec(xr) ));
1222         if ( ~scaledc ) then
1223           scaledc = %t;
1224         end
1225       end
1226     end
1227     // If scaling failed, returns immediately 
1228     // (there is no need to update the function value).
1229     if ( ~status ) then
1230       return
1231     end
1232     if ( scaledc ) then
1233       // Re-compute the function value at scaled point
1234       [ this.optbase , fr , cr , index ] = optimbase_function ( this.optbase , xr , 2 );
1235     end
1236     
1237   endfunction