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