Optimization: let tests pass
[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 // Copyright (C) 2009-2011 - DIGITEO - Michael Baudin
4 //
5 // Copyright (C) 2012 - 2016 - Scilab Enterprises
6 //
7 // This file is hereby licensed under the terms of the GNU GPL v2.0,
8 // pursuant to article 5.3.4 of the CeCILL v.2.1.
9 // This file was originally licensed under the terms of the CeCILL v2.1,
10 // and continues to be available under such terms.
11 // For more information, see the COPYING file which you should have received
12 // along with this program.
13
14 //
15 // References
16 //   Sequential Application of Simplex Designs in Optimisation and Evolutionary Operation
17 //   Spendley, W. and Hext, G. R. and Himsworth, F. R.
18 //   1962
19 //
20 //   Convergence Properties of the Nelder-Mead Simplex Method in Low Dimensions
21 //   Jeffrey C. Lagarias and James A. Reeds nd Margaret H. Wright and Paul E. Wright
22 //   SIAM Journal of Optimization, 1998, volume 9, pp. 112--147
23 //
24 //   A Simplex Method for Function Minimization
25 //   Nelder, J. A.  and Mead, R.
26 //   1965
27 //
28 //   Iterative Methods for Optimization
29 //   C. T. Kelley,
30 //   SIAM Frontiers in Applied Mathematics
31 //   1999
32 //
33 //   Detection and Remediation of Stagnation in the Nelder--Mead Algorithm Using a Sufficient Decrease Condition
34 //   Kelley C. T.
35 //   SIAM J. on Optimization},
36 //   1999
37 //
38
39 //
40 // neldermead_search --
41 //   Search the minimum with Nelder-Mead algorithm.
42 //
43 function this = neldermead_search ( this , warn)
44
45     rhs = argn(2);
46
47     warn_mode = warning("query"); // Saving current warning mode to revert to it at the end
48     if rhs == 2 then
49         if warn == 0 | warn == "off" then
50             warning("off"); // Turn off warnings
51         elseif warn == 1 | warn == "on" then
52             warning("on"); // Turn on warnings
53         else
54             msg = _("%s: Wrong value for input argument #%d: ""%s"", ""%s"", %d or %d expected.\n");
55             error(msprintf(msg, "neldermead_search", 2, "off", "on", 0, 1))
56         end
57     end
58
59     withderivatives = optimbase_cget ( this.optbase , "-withderivatives" );
60     if ( withderivatives ) then
61         errmsg = msprintf(gettext("%s: The -withderivatives option is true but all algorithms in neldermead are derivative-free."), "neldermead_search")
62         error(errmsg)
63     end
64     if ( ~this.startupflag ) then
65         this = neldermead_startup ( this );
66         this.startupflag = %t;
67     end
68     stop = neldermead_outputcmd ( this, "init" , this.simplex0 , "init" )
69     if ( stop ) then
70         this.optbase = optimbase_set ( this.optbase , "-status" , "userstop" );
71         verbose = optimbase_cget ( this.optbase , "-verbose" )
72         if ( verbose == 1 ) then
73             this = neldermead_log (this,sprintf("Terminate by user''s request"));
74         end
75         return
76     end
77
78     if ( this.restartflag ) then
79         this = neldermead_autorestart ( this )
80     else
81         this = neldermead_algo ( this );
82     end
83     stop = neldermead_outputcmd ( this, "done" , this.simplexopt , "done" )
84     if ( stop ) then
85         this.optbase = optimbase_set ( this.optbase , "-status" , "userstop" );
86         verbose = optimbase_cget ( this.optbase , "-verbose" )
87         if ( verbose == 1 ) then
88             this = neldermead_log (this,sprintf("Terminate by user''s request"));
89         end
90     end
91
92     warning(warn_mode); // Revert to the warning mode that was active before the call
93
94 endfunction
95 //
96 // neldermead_algo --
97 //   Performs an optimization without restart
98 //
99 function this = neldermead_algo ( this )
100     select this.method
101     case "fixed" then
102         this = neldermead_fixed (this);
103     case "variable" then
104         this = neldermead_variable (this);
105     case "box" then
106         this = neldermead_box (this);
107     case "mine" then
108         this = this.mymethod ( this );
109     else
110         errmsg = msprintf(gettext("%s: Unknown -method %s"), "neldermead_algo", this.method)
111         error(errmsg)
112     end
113 endfunction
114 //
115 // neldermead_autorestart --
116 //   Performs an optimization with automatic restart
117 // NOTE
118 //   The loop processes for i = 1 to restartmax PLUS 1
119 //   This is because a RE-start is performed after one simulation
120 //   has been performed, hence the "RE".
121 //   Hence,
122 //     * if restartmax = 0, the number of performed loops is 1.
123 //     * if restartmax = 1, the number of performed loops is 2.
124 //     * etc...
125 // Example #1
126 //   Sample session with restart max = 3 (the default) and process pass.
127 //   restartnb = 0
128 //   Try #1/4
129 //     Run
130 //     "Must restart"
131 //     restartnb = 1
132 //   Try #2/4
133 //     Run
134 //     "Must restart"
135 //     restartnb = 2
136 //   Try #3/4
137 //     Run
138 //     "Must not restart"
139 //   "Convergence reached after 2 restarts."
140 //   reached = %t
141 //   status = depending on the triggered termination criteria
142 //   restartnb = 2 (restarts are Tries #2 and #3)
143 // Example #2
144 //   Sample session with restart max = 3 (the default) and process fails.
145 //   restartnb = 0
146 //   Try #1/4
147 //     Run
148 //     "Must restart"
149 //     restartnb = 1
150 //   Try #2/4
151 //     Run
152 //     "Must restart"
153 //     restartnb = 2
154 //   Try #3/4
155 //     Run
156 //     "Must restart"
157 //     restartnb = 3
158 //   Try #4/4
159 //     Run
160 //     Must restart
161 //   "Convergence not reached after maximum 3 restarts."
162 //   reached = %f
163 //   status = "maxrestart"
164 //   restartnb = 3 (restarts are Tries #2, #3 and #4)
165 //
166 function this = neldermead_autorestart ( this )
167     restartmax = this.restartmax;
168     reached = %f;
169     for iloop = 1: restartmax + 1
170         this = neldermead_log (this,sprintf("*****************************************************************"));
171         this = neldermead_log (this,sprintf("Try #%d/%d.", iloop , restartmax + 1 ));
172         //
173         // Run algorithm
174         this = neldermead_algo ( this );
175         //
176         // Must we restart ?
177         [ this , istorestart ] = neldermead_istorestart ( this );
178         if ( istorestart ) then
179             this = neldermead_log (this,"Must restart.");
180         else
181             this = neldermead_log (this,"Must not restart.");
182         end
183         if ( ~istorestart ) then
184             reached = %t;
185             break
186         end
187         if ( iloop < restartmax + 1 ) then
188             // We are going to perform a restart
189             this.restartnb = this.restartnb + 1;
190             this = neldermead_log (this,"Updating simplex.");
191             this = neldermead_updatesimp ( this );
192         end
193     end
194     if ( reached ) then
195         this = neldermead_log (this, sprintf ( "Convergence reached after %d restarts." , this.restartnb ) );
196     else
197         this = neldermead_log (this, sprintf ( "Convergence not reached after maximum %d restarts." , this.restartnb ) );
198         this.optbase = optimbase_set ( this.optbase , "-status" , "maxrestart" );
199     end
200 endfunction
201
202
203 //
204 // neldermead_variable --
205 //   The original Nelder-Mead algorithm, with variable-size simplex.
206 //
207 function this = neldermead_variable ( this )
208     // Check settings correspond to algo
209     [ this.optbase , hascons ] = optimbase_hasconstraints ( this.optbase );
210     if ( hascons ) then
211         errmsg = msprintf(gettext("%s: Problem has constraints, but variable algorithm ignores them."), "neldermead_variable")
212         error(errmsg)
213     end
214     verbose = optimbase_cget ( this.optbase , "-verbose" )
215     //
216     // Order the vertices for the first time
217     //
218     simplex = this.simplex0;
219     n = optimbase_cget ( this.optbase , "-numberofvariables" );
220     if (n==0) then
221         errmsg = msprintf(gettext("%s: The number of variable is zero."), "neldermead_variable")
222         error(errmsg)
223     end
224     fvinitial = optimbase_get ( this.optbase , "-fx0" );
225     // Sort function values and x points by increasing function value order
226     this = neldermead_log (this,"Step #1 : order");
227     simplex = optimsimplex_sort ( simplex );
228     // Transpose, because optimsimplex returns row vectors
229     currentcenter = optimsimplex_center ( simplex ).';
230     currentxopt = optimbase_cget ( this.optbase , "-x0" );
231     newfvmean = optimsimplex_fvmean ( simplex );
232     greedy = this.greedy;
233     //
234     // Initialize
235     //
236     terminate = %f;
237     iter = 0;
238     step = "init";
239     //
240     // Nelder-Mead Loop
241     //
242     while ( ~terminate )
243         this.optbase = optimbase_incriter ( this.optbase );
244         iter = iter + 1;
245         // Transpose, because optimsimplex returns row vectors
246         xlow = optimsimplex_getx ( simplex , 1 ).'
247         flow = optimsimplex_getfv ( simplex , 1 )
248         xhigh = optimsimplex_getx ( simplex , n+1 ).'
249         fhigh = optimsimplex_getfv ( simplex , n+1 )
250         xn = optimsimplex_getx ( simplex , n ).'
251         fn = optimsimplex_getfv ( simplex , n )
252         //
253         // Store history
254
255
256
257
258
259         //
260         xcoords = optimsimplex_getallx ( simplex )
261         this = neldermead_storehistory ( this , n , flow , xlow , xcoords );
262         currentfopt = flow;
263         previousxopt = currentxopt;
264         currentxopt = xlow;
265         previouscenter = currentcenter;
266         currentcenter = optimsimplex_center ( simplex ).';
267         oldfvmean = newfvmean;
268         newfvmean = optimsimplex_fvmean ( simplex );
269         if ( verbose == 1 ) then
270             this = neldermead_logsummary ( this, iter,xlow,flow,currentcenter,simplex )
271         end
272         this.optbase = optimbase_set ( this.optbase , "-xopt" , xlow );
273         this.optbase = optimbase_set ( this.optbase , "-fopt" , flow );
274         stop = neldermead_outputcmd ( this, "iter" , simplex , step )
275         if ( stop ) then
276             status = "userstop"
277             if ( verbose == 1 ) then
278                 this = neldermead_log (this,sprintf("Terminate by user''s request"));
279             end
280             break
281         end
282
283         //
284         // Update termination flag
285         //
286         if ( iter > 1 ) then
287             [ this , terminate , status ] = neldermead_termination (this , ..
288             fvinitial , oldfvmean , newfvmean , previouscenter , currentcenter , simplex );
289             if ( terminate ) then
290                 this = neldermead_log (this,sprintf("Terminate with status : %s",status));
291                 break
292             end
293         end
294         //
295         // Compute xbar, center of better vertices
296         //
297         if ( verbose == 1 ) then
298             this = neldermead_log (this,sprintf("Reflect"));
299         end
300         // Transpose, because optimsimplex returns row vectors
301         xbar   = optimsimplex_xbar ( simplex ).';
302         if ( verbose == 1 ) then
303             this = neldermead_log (this,sprintf("xbar="+_strvec(xbar)+""));
304         end
305         //
306         // Reflect the worst point with respect to center
307         //
308         xr = neldermead_interpolate ( xbar , xhigh , this.rho );
309         [ this.optbase , fr , index ] = optimbase_function ( this.optbase , xr , 2 );
310         if ( verbose == 1 ) then
311             this = neldermead_log (this,sprintf("xr=["+_strvec(xr)+"], f(xr)=%f",fr));
312         end
313         if ( fr >= flow & fr < fn ) then
314             if ( verbose == 1 ) then
315                 this = neldermead_log (this,sprintf("  > Perform reflection"));
316             end
317             simplex = optimsimplex_setve ( simplex , n+1 , fr , xr.' )
318             step = "reflection";
319         elseif ( fr < flow ) then
320             // Expand
321             if ( verbose == 1 ) then
322                 this = neldermead_log (this,sprintf("Expand"));
323             end
324             xe = neldermead_interpolate ( xbar , xhigh , this.rho*this.chi );
325             [ this.optbase , fe , index ] = optimbase_function ( this.optbase , xe , 2 );
326             if ( verbose == 1 ) then
327                 this = neldermead_log (this,sprintf("xe="+_strvec(xe)+", f(xe)=%f",fe));
328             end
329             if ( greedy ) then
330                 if ( fe < flow ) then
331                     if ( verbose == 1 ) then
332                         this = neldermead_log (this,sprintf("  > Perform Greedy Expansion"));
333                     end
334                     simplex = optimsimplex_setve ( simplex , n+1 , fe , xe.' )
335                     step = "expansion";
336                 else
337                     if ( verbose == 1 ) then
338                         this = neldermead_log (this,sprintf("  > Perform Greedy Reflection"));
339                     end
340                     simplex = optimsimplex_setve ( simplex , n+1 , fr , xr.' )
341                     step = "reflection";
342                 end
343             else
344                 if ( fe < fr ) then
345                     if ( verbose == 1 ) then
346                         this = neldermead_log (this,sprintf("  > Perform Expansion"));
347                     end
348                     simplex = optimsimplex_setve ( simplex , n+1 , fe , xe.' )
349                     step = "expansion";
350                 else
351                     if ( verbose == 1 ) then
352                         this = neldermead_log (this,sprintf("  > Perform Reflection"));
353                     end
354                     simplex = optimsimplex_setve ( simplex , n+1 , fr , xr.' )
355                     step = "reflection";
356                 end
357             end
358         elseif ( fr >= fn & fr < fhigh ) then
359             // Outside contraction
360             if ( verbose == 1 ) then
361                 this = neldermead_log (this,sprintf("Contract - outside"));
362             end
363             xc = neldermead_interpolate ( xbar , xhigh , this.rho*this.gamma );
364             [ this.optbase , fc , index ] = optimbase_function ( this.optbase , xc , 2 );
365             if ( verbose == 1 ) then
366                 this = neldermead_log (this,sprintf("xc="+_strvec(xc)+", f(xc)=%f",fc));
367             end
368             if ( fc <= fr ) then
369                 if ( verbose == 1 ) then
370                     this = neldermead_log (this,sprintf("  > Perform Outside Contraction"));
371                 end
372                 simplex = optimsimplex_setve ( simplex , n+1 , fc , xc.' )
373                 step = "outsidecontraction";
374             else
375                 //  Shrink
376                 if ( verbose == 1 ) then
377                     this = neldermead_log (this,sprintf("  > Perform Shrink"));
378                 end
379                 [ simplex , this ] = optimsimplex_shrink ( simplex , costf_transposex , this.sigma , this );
380                 step = "shrink";
381             end
382         else
383             // ( fr >= fn & fr >= fhigh )
384             // Inside contraction
385             if ( verbose == 1 ) then
386                 this = neldermead_log (this,sprintf("Contract - inside"));
387             end
388             xc = neldermead_interpolate ( xbar , xhigh , -this.gamma );
389             [ this.optbase , fc , index ] = optimbase_function ( this.optbase , xc , 2 );
390             if ( verbose == 1 ) then
391                 this = neldermead_log (this,sprintf("xc="+_strvec(xc)+", f(xc)=%f",fc));
392             end
393             if ( fc < fhigh ) then
394                 if ( verbose == 1 ) then
395                     this = neldermead_log (this,sprintf("  > Perform Inside Contraction"));
396                 end
397                 simplex = optimsimplex_setve ( simplex , n+1 , fc , xc.' )
398                 step = "insidecontraction";
399             else
400                 //  Shrink
401                 if ( verbose == 1 ) then
402                     this = neldermead_log (this,sprintf("  > Perform Shrink"));
403                 end
404                 [ simplex , this ] = optimsimplex_shrink ( simplex , costf_transposex , this.sigma , this )
405                 step = "shrink";
406             end
407         end
408         //
409         // Sort simplex
410         //
411         if ( verbose == 1 ) then
412             this = neldermead_log (this,sprintf("Sort"));
413         end
414         simplex  = optimsimplex_sort ( simplex );
415     end
416     this.optbase = optimbase_set ( this.optbase , "-xopt" , xlow );
417     this.optbase = optimbase_set ( this.optbase , "-fopt" , flow );
418     this.optbase = optimbase_set ( this.optbase , "-status" , status );
419     this.simplexopt = simplex;
420 endfunction
421
422 //
423 // neldermead_fixed --
424 //   The simplex algorithm with fixed size simplex.
425 // Implementation note:
426 //   We implement the following "rules" of the Spendley et al.
427 //   method.
428 //   Rule 1 is strictly applied, but the reflection is done
429 //   by reflection the high point, since we minimize a function
430 //   instead of maximizing it, like Spendley.
431 //   The Rule 2 is NOT implemented, as we expect that the
432 //   function evaluation is not subject to errors.
433 //   The Rule 3 is applied, ie reflection with respect
434 //   to next to high point.
435 //   A shrink step is included, with shrinkage factor sigma.
436 //
437 //   "Rule 1. Ascertain the lowest reading y, of yi ... Yk+1
438 //   Complete a new simplex Sp by excluding the point Vp corresponding to
439 //   y, and replacing it by V* defined as above."
440 //   "Rule 2. If a result has occurred in (k + 1) successive simplexes, and is not
441 //   then eliminated by application of Rule 1, do not move in the direction
442 //   indicated by Rule 1, or at all, but discard the result and replace it by a new
443 //   observation at the same point."
444 //   "Rule 3. If y is the lowest reading in So , and if the next observation made,
445 //   y* , is the lowest reading in the new simplex S , do not apply Rule 1 and
446 //   return to So from Sp . Move out of S, by rejecting the second lowest reading
447 //   (which is also the second lowest reading in So)."
448 //
449 function this = neldermead_fixed (this)
450     // Check settings correspond to algo
451     [ this.optbase , hascons ] = optimbase_hasnlcons ( this.optbase );
452     if ( hascons ) then
453         errmsg = msprintf(gettext("%s: Problem has constraints, but fixed algorithm ignores them."), "neldermead_fixed")
454         error(errmsg)
455     end
456     verbose = optimbase_cget ( this.optbase , "-verbose" )
457     //
458     // Order the vertices for the first time
459     //
460     simplex = this.simplex0;
461     n = optimbase_cget ( this.optbase , "-numberofvariables" );
462     fvinitial = optimbase_get ( this.optbase , "-fx0" );
463     // Sort function values and x points by increasing function value order
464     this = neldermead_log (this,sprintf("Sort"));
465     simplex = optimsimplex_sort ( simplex );
466     //
467     // Compute center of simplex
468     //
469     // Transpose, because optimsimplex returns row vectors
470     currentcenter = optimsimplex_center ( simplex ).';
471     newfvmean = optimsimplex_fvmean ( simplex );
472     currentxopt = optimbase_cget ( this.optbase , "-x0" );
473     //
474     // Set indices for "clarity"
475     //
476     ilow = 1
477     ihigh = n + 1
478     inext = n
479     //
480     // Initialize
481     //
482     terminate = %f;
483     iter = 0;
484     step = "init";
485     //
486     // main N-M loop
487     //
488     while ( ~terminate )
489         this.optbase = optimbase_incriter ( this.optbase );
490         iter = iter + 1;
491         xlow = optimsimplex_getx ( simplex , ilow ).'
492         flow = optimsimplex_getfv ( simplex , ilow )
493         xhigh = optimsimplex_getx ( simplex , ihigh ).'
494         fhigh = optimsimplex_getfv ( simplex , ihigh )
495         //
496         // Store history
497         //
498         xcoords = optimsimplex_getallx ( simplex )
499         this = neldermead_storehistory ( this , n , flow , xlow , xcoords );
500         currentfopt = flow;
501         previousxopt = currentxopt;
502         currentxopt = xlow;
503         previouscenter = currentcenter;
504         currentcenter = optimsimplex_center ( simplex ).';
505         oldfvmean = newfvmean;
506         newfvmean = optimsimplex_fvmean ( simplex );
507         if ( verbose == 1 ) then
508             this = neldermead_logsummary ( this, iter,xlow,flow,currentcenter,simplex )
509         end
510         this.optbase = optimbase_set ( this.optbase , "-xopt" , xlow );
511         this.optbase = optimbase_set ( this.optbase , "-fopt" , flow );
512         stop = neldermead_outputcmd ( this, "iter" , simplex , step )
513         if ( stop ) then
514             status = "userstop"
515             if ( verbose == 1 ) then
516                 this = neldermead_log (this,sprintf("Terminate by user''s request"));
517             end
518             break
519         end
520         //
521         // Update termination flag
522         //
523         if ( iter > 1 ) then
524             [ this , terminate , status] = neldermead_termination (this , ..
525             fvinitial , oldfvmean , newfvmean , previouscenter , currentcenter , simplex );
526             if ( terminate ) then
527                 if ( verbose == 1 ) then
528                     this = neldermead_log (this,sprintf("Terminate with status : %s",status));
529                 end
530                 break;
531             end
532         end
533         //
534         // Compute xbar, center of better vertices
535         //
536         if ( verbose == 1 ) then
537             this = neldermead_log (this,sprintf("Reflect"));
538         end
539         xbar = optimsimplex_xbar ( simplex ).';
540         if ( verbose == 1 ) then
541             this = neldermead_log (this,sprintf("xbar="+_strvec(xbar)+""));
542         end
543         //
544         // Reflect the worst point with respect to center
545         //
546         xr = neldermead_interpolate ( xbar , xhigh , this.rho );
547         [ this.optbase , fr , index ] = optimbase_function ( this.optbase , xr , 2 );
548         if ( verbose == 1 ) then
549             this = neldermead_log (this,sprintf("xr="+_strvec(xr)+", f(xr)=%f",fr));
550         end
551         //
552         // Replace worst point by xr if it is better
553         //
554         if ( fr < fhigh ) then
555             if ( verbose == 1 ) then
556                 this = neldermead_log (this,sprintf("  > Perform reflect"));
557             end
558             simplex = optimsimplex_setve ( simplex , ihigh , fr , xr.' )
559             step = "reflection";
560         else
561             // Reflect / xnext
562             xnext = optimsimplex_getx ( simplex , inext ).';
563             fnext = optimsimplex_getfv ( simplex , inext );
564             xbar2 = optimsimplex_xbar ( simplex , inext ).';
565             if ( verbose == 1 ) then
566                 this = neldermead_log (this,sprintf("xbar2="+_strvec(xbar2)+""));
567             end
568             xr2 = neldermead_interpolate ( xbar2 , xnext , this.rho );
569             [ this.optbase , fr2 , index ] = optimbase_function ( this.optbase , xr2 , 2 );
570             if ( verbose == 1 ) then
571                 this = neldermead_log (this,sprintf("xr2="+_strvec(xr2)+", f(xr2)=%f",fr2));
572             end
573             if ( fr2 < fnext ) then
574                 if ( verbose == 1 ) then
575                     this = neldermead_log (this,sprintf("  > Perform reflect / next"));
576                 end
577                 simplex = optimsimplex_setve ( simplex , inext , fr2 , xr2.' )
578                 step = "reflectionnext";
579             else
580                 //  Shrink
581                 if ( verbose == 1 ) then
582                     this = neldermead_log (this,sprintf("  > Perform Shrink"));
583                 end
584                 [ simplex , this ] = optimsimplex_shrink ( simplex , costf_transposex , this.sigma , this )
585                 step = "shrink";
586             end
587         end
588         //
589         // Sort simplex
590         //
591         simplex = optimsimplex_sort ( simplex );
592     end
593     this.optbase = optimbase_set ( this.optbase , "-xopt" , xlow );
594     this.optbase = optimbase_set ( this.optbase , "-fopt" , flow );
595     this.optbase = optimbase_set ( this.optbase , "-status" , status );
596     this.simplexopt = simplex;
597 endfunction
598 //
599 // neldermead_interpolate --
600 //   Computes xi as an interpolation between x1 and x2, with factor as :
601 //     xi = (1+fac)x1 - fac x2
602 //
603 function xi = neldermead_interpolate ( x1 , x2 , fac )
604     xi = (1 + fac)*x1 - fac*x2;
605 endfunction
606
607 //
608 // neldermead_termination --
609 //   Returns 1 if the algorithm terminates.
610 //   Returns 0 if the algorithm must continue.
611 // Arguments
612 //   this : the current object
613 //   fvinitial : initial function value
614 //   newfvmean, oldfvmean : the old and new function value average on the simplex
615 //   previousxopt : the previous value of x
616 //   currentxopt : the current value of x
617 //   simplex : the simplex
618 //     The best point in the simplex is expected to be stored at 1
619 //     The worst point in the simplex is expected to be stored at n+1
620 //   terminate : %t if the algorithm terminates, %f if the algorithm must continue.
621 //   this.status : termination status
622 //     "continue"
623 //     "maxiter"
624 //     "maxfuneval"
625 //     "tolf"
626 //     "tolx"
627 //     "tolsize"
628 //     "tolsizedeltafv"
629 //     "kelleystagnation"
630 //     "tolboxf"
631 //     "tolvariance"
632 // Notes
633 //   Use the function average on the simplex instead of the best function value.
634 //   This is because the function average changes at each iteration.
635 //   Instead, the best function value as a step-by-step evolution and may not
636 //   change in 2 iterations, leading to astop of the algorithm.
637 // TODO : set the fvinitial, oldfvmean, newfvmean.
638 //
639 function [ this , terminate , status ] = neldermead_termination (this , ..
640     fvinitial , oldfvmean , newfvmean , previousxopt , currentxopt , ..
641     simplex )
642     terminate = %f;
643     status = "continue";
644     verbose = optimbase_cget ( this.optbase , "-verbose" )
645     //
646     // Termination Criteria from parent optimization class
647     //
648     [ this.optbase , terminate , status ] = optimbase_terminate ( this.optbase , ..
649     fvinitial , newfvmean , previousxopt , currentxopt );
650     //
651     // Criteria #6 : simplex absolute + relative size
652     //
653     if ( ~terminate ) then
654         if ( this.tolsimplexizemethod ) then
655             ssize = optimsimplex_size ( simplex , "sigmaplus" );
656             tolsa = this.tolsimplexizeabsolute;
657             tolsr = this.tolsimplexizerelative;
658             ssize0 = this.simplexsize0;
659             if ( verbose == 1 ) then
660                 this.optbase = optimbase_stoplog  ( this.optbase,sprintf("  > simplex size=%s < %s + %s * %s",..
661                 string(ssize), string(tolsa) , string(tolsr) , string(ssize0) ));
662             end
663             if ( ssize < tolsa + tolsr * ssize0 ) then
664                 terminate = %t;
665                 status = "tolsize";
666             end
667         end
668     end
669     //
670     // Criteria #7 : simplex absolute size + difference in function values (Matlab-like)
671     //
672     if ( ~terminate ) then
673         if ( this.tolssizedeltafvmethod ) then
674             ssize = optimsimplex_size ( simplex , "sigmaplus" );
675             if ( verbose == 1 ) then
676                 this.optbase = optimbase_stoplog  ( this.optbase,sprintf("  > simplex size=%s < %s",..
677                 string(ssize), string(this.tolsimplexizeabsolute)));
678             end
679             shiftfv = abs(optimsimplex_deltafvmax( simplex ))
680             if ( verbose == 1 ) then
681                 this.optbase = optimbase_stoplog  ( this.optbase,sprintf("  > abs(fv(n+1) - fv(1))=%s < toldeltafv=%s",..
682                 string(shiftfv), string(this.toldeltafv)));
683             end
684             if ( ( ssize < this.tolsimplexizeabsolute ) & ( shiftfv < this.toldeltafv ) ) then
685                 terminate = %t;
686                 status = "tolsizedeltafv";
687             end
688         end
689     end
690     //
691     // Criteria #8 : Kelley stagnation, based on
692     // a sufficient decrease condition
693     //
694     if ( ~terminate ) then
695         if ( this.kelleystagnationflag ) then
696             [ sg , this ] = optimsimplex_gradientfv ( simplex , neldermead_costf , "forward" , this );
697             nsg = sg.' * sg;
698             if ( verbose == 1 ) then
699                 sgstr = _strvec(sg);
700                 this.optbase = optimbase_stoplog ( this.optbase , sprintf ( "Test Stagnation : nsg = %s, sg = [%s]", ..
701                 string(nsg) , sgstr ) );
702                 this.optbase = optimbase_stoplog ( this.optbase , ..
703                 sprintf ( "Test Stagnation : newfvmean=%s >= oldfvmean=%s - %s * %s" , ..
704                 string(newfvmean), string(oldfvmean) , string(this.kelleyalpha) , string(nsg) ) );
705             end
706             if ( newfvmean >= oldfvmean - this.kelleyalpha * nsg ) then
707                 terminate = %t;
708                 status = "kelleystagnation";
709             end
710         end
711     end
712     //
713     // Criteria #9 : Box termination criteria
714     // The number of consecutive time that an absolute tolerance on
715     // function value is met.
716     // From Algorithm 454, the tolerance is the difference between the
717     // max and the min function values in the simplex.
718     //
719     if ( ~terminate ) then
720         if ( this.boxtermination ) then
721             shiftfv = abs(optimsimplex_deltafvmax( simplex ))
722             if ( verbose == 1 ) then
723                 this.optbase = optimbase_stoplog ( this.optbase , ..
724                 sprintf ( "Test Box : shiftfv=%s < boxtolf=%s" , ..
725                 string(shiftfv) , string(this.boxtolf) ) );
726             end
727             if ( shiftfv < this.boxtolf ) then
728                 this.boxkount = this.boxkount + 1
729                 if ( verbose == 1 ) then
730                     this.optbase = optimbase_stoplog ( this.optbase , ..
731                     sprintf ( "Test Box : boxkount=%d == boxnbmatch=%d" , ..
732                     this.boxkount , this.boxnbmatch ) );
733                 end
734                 if ( this.boxkount == this.boxnbmatch ) then
735                     terminate = %t
736                     status = "tolboxf"
737                 end
738             else
739                 this.boxkount = 0
740             end
741         end
742     end
743     //
744     // Obsolete option: maintain for backward compatibility
745     // Criteria #10 : variance of function values
746     //
747     if ( ~terminate ) then
748         if ( this.tolvarianceflag ) then
749             var = optimsimplex_fvvariance ( simplex )
750             if ( verbose == 1 ) then
751                 this.optbase = optimbase_stoplog ( this.optbase , ..
752                 sprintf ( "Test tolvariance : %s < %s" , ..
753                 string(var) , string(this.tolabsolutevariance) ) );
754             end
755             if ( var < this.tolrelativevariance * this.variancesimplex0 + this.tolabsolutevariance ) then
756                 terminate = %t
757                 status = "tolvariance"
758             end
759         end
760     end
761     //
762     // Obsolete option: maintain for backward compatibility
763     // Criteria #11 : user-defined criteria
764     //
765     if ( ~terminate ) then
766         if ( this.myterminateflag ) then
767             [ this , term , stat ] = this.myterminate ( this , simplex )
768             if ( term ) then
769                 terminate = term
770                 status = stat
771             end
772         end
773     end
774     //
775     // A verbose message
776     //
777     if ( verbose == 1 ) then
778         this.optbase = optimbase_stoplog (this.optbase,sprintf("  > Terminate = %s, status = %s",..
779         string(terminate) , status ));
780     end
781 endfunction
782
783
784 //
785 // neldermead_outputcmd --
786 //   Calls back user's output command
787 // Arguments
788 //   this : the current object
789 //   state : the state of the algorithm,
790 //     "init", "done", "iter"
791 //   simplex : the current simplex
792 //   step : the type of step performed during the iteration
793 //     "init", "done", "reflection", "expansion", "insidecontraction", "outsidecontraction"
794 //     "reflectionnext", "shrink"
795 //   stop : set to true to stop algorithm
796 //
797 function stop = neldermead_outputcmd ( this, ..
798     state , simplex , step )
799     outputcmd = optimbase_cget ( this.optbase , "-outputcommand" );
800     if typeof(outputcmd) <> "string" then
801         brutedata = optimbase_outstruct ( this.optbase );
802         data = tlist(["T_NMDATA",..
803         "x","fval","iteration","funccount",..
804         "simplex" , "step" ]);
805         data.x = brutedata.x;
806         data.fval = brutedata.fval;
807         data.iteration = brutedata.iteration;
808         data.funccount = brutedata.funccount;
809         data.simplex = simplex;
810         data.step = step;
811         stop = optimbase_outputcmd ( this.optbase , state , data );
812     else
813         stop = %f
814     end
815 endfunction
816 //
817 // neldermead_storehistory --
818 //   Stores the history into the data structure.
819 // Arguments, input
820 //   this : current object
821 //   n : the number of unknown parameters
822 //   fopt : the current value of the function at optimum
823 //   xopt : arrary with size n, current optimum
824 //   xcoords : array with size n x n+1, coordinates of the n+1 vertices
825 //
826 function this = neldermead_storehistory ( this , n , fopt , xopt , xcoords )
827     storehistory = optimbase_cget ( this.optbase , "-storehistory" );
828     iterations = optimbase_get ( this.optbase , "-iterations" );
829     if ( storehistory ) then
830         this.optbase = optimbase_histset ( this.optbase , iterations , "-fopt" , fopt );
831         this.optbase = optimbase_histset ( this.optbase , iterations , "-xopt" , xopt(1:n) );
832         this.historysimplex ( iterations , 1:n+1,1:n) = xcoords(1:n+1,1:n);
833     end
834 endfunction
835
836 //
837 // neldermead_istorestart --
838 //   Returns 1 if the optimization is to restart.
839 // Arguments
840 //   istorestart : %t of the optimization is to restart.
841 //
842 function [ this , istorestart ] = neldermead_istorestart ( this )
843     status = optimbase_get ( this.optbase , "-status" );
844     if ( status =="maxfuneval" ) then
845         istorestart = %f
846         return
847     end
848     select this.restartdetection
849     case "oneill"
850         [ this , istorestart ] = neldermead_isroneill ( this )
851     case "kelley"
852         [ this , istorestart ] = neldermead_isrkelley ( this )
853     else
854         errmsg = msprintf(gettext("%s: Unknown restart detection %s"),"neldermead_istorestart", this.restartdetection)
855         error(errmsg)
856     end
857 endfunction
858 //
859 // neldermead_isrkelley--
860 //   Returns 1 if the optimization is to restart.
861 //   Use kelleystagnation as a criteria for restart.
862 // Arguments
863 //   istorestart : %t of the optimization is to restart.
864 //
865 function [ this , istorestart ] = neldermead_isrkelley ( this )
866     istorestart = %f
867     if ( this.kelleystagnationflag ) then
868         status = optimbase_get ( this.optbase , "-status" );
869         if ( status =="kelleystagnation" ) then
870             istorestart = %t
871         end
872     end
873 endfunction
874 //
875 // neldermead_isroneill --
876 //   Returns 1 if the optimization is to restart.
877 //   Use O'Neill method as a criteria for restart.
878 //   It is an axis by axis search for optimality.
879 // Arguments
880 //   xopt : the optimum, as a st of n values
881 //   fopt : function value at optimum
882 //   eps : a small value
883 //   step : a list of n values, representing
884 //     the "size" of each parameter
885 //   istorestart : %t if the optimization is to restart.
886 //
887 function [ this , istorestart ] = neldermead_isroneill ( this )
888     n = optimbase_cget ( this.optbase , "-numberofvariables" );
889     //
890     // If required, make a vector step from the scalar step
891     //
892     defaultstep = this.restartstep;
893     steprows = size ( defaultstep , "r" );
894     if ( steprows == 1 ) then
895         step = defaultstep * ones(n,1);
896     else
897         step = defaultstep;
898     end
899     restarteps = this.restarteps;
900
901     x = optimbase_get ( this.optbase , "-xopt" );
902     fopt = optimbase_get ( this.optbase , "-fopt" );
903     verbose = optimbase_cget ( this.optbase , "-verbose" )
904
905     if ( verbose ) then
906         this = neldermead_log (this,sprintf("================================================================="));
907         this = neldermead_log (this, sprintf ( "O''Neill Restart\n") );
908         this = neldermead_log (this, sprintf ( "Using step [%s]" , _strvec(step) ) );
909     end
910
911     istorestart = %f
912     for ix = 1:n
913         stepix = step ( ix )
914         del = stepix * restarteps
915         xix =  x ( ix )
916         x ( ix ) = xix + del
917         [ this.optbase , fv , index ] = optimbase_function ( this.optbase , x , 2 )
918         if ( fv < fopt ) then
919             istorestart = %t
920             if ( verbose ) then
921                 this = neldermead_log (this, sprintf ( "Must restart because fv=%s at [%s] is lower than fopt=%s" , ..
922                 string(fv) , _strvec(x) , string(fopt)) );
923             end
924             break
925         end
926         x ( ix ) = xix - del
927         [ this.optbase , fv , index ] = optimbase_function ( this.optbase , x , 2 )
928         if ( fv < fopt ) then
929             istorestart = %t
930             if ( verbose ) then
931                 this = neldermead_log (this, sprintf( "Must restart because fv=%s at [%s] is lower than fopt=%s" , ..
932                 string(fv) , _strvec(x) , string(fopt)) );
933             end
934             break
935         end
936         x ( ix ) = xix
937     end
938 endfunction
939
940 //
941 // neldermead_startup --
942 //   Startup the algorithm.
943 //   Computes the initial simplex, depending on the -simplex0method.
944 //
945 function this = neldermead_startup (this)
946     //
947     // 0. Check that the cost function is correctly connected
948     // Note: this call to the cost function is not used, but helps the
949     // user while he is tuning his object.
950     if ( this.checkcostfunction ) then
951         this.optbase = optimbase_checkcostfun ( this.optbase );
952     end
953     //
954     // 1. If the problem has bounds, check that they are consistent
955     [ this.optbase , hasbounds ] = optimbase_hasbounds ( this.optbase );
956     if ( hasbounds ) then
957         this.optbase = optimbase_checkbounds ( this.optbase );
958     end
959     //
960     // 2. Get the initial guess and compute the initial simplex
961     x0 = optimbase_cget ( this.optbase , "-x0" );
962     select this.simplex0method
963     case "given" then
964         [ simplex0 , this ] = optimsimplex_new ( this.coords0 , [] , this );
965     case "axes" then
966         [ simplex0 , this ] = optimsimplex_new ( "axes" , ..
967         x0.' , [] , this.simplex0length , this );
968     case "spendley" then
969         [ simplex0 , this ] = optimsimplex_new ( "spendley" , ..
970         x0.' , [] , this.simplex0length , this );
971     case "pfeffer" then
972         [ simplex0 , this ] = optimsimplex_new ( "pfeffer" , ..
973         x0.' , [] , this.simplex0deltausual , ..
974         this.simplex0deltazero , this );
975     case "randbounds" then
976         if ( this.boxnbpoints == "2n" ) then
977             this.boxnbpointseff = 2 * this.optbase.numberofvariables;
978         else
979             this.boxnbpointseff = this.boxnbpoints;
980         end
981         if ( ~hasbounds ) then
982             error ( msprintf(gettext("%s: Randomized bounds initial simplex is not available without bounds." ), "neldermead_startup"))
983         end
984         [ simplex0 , this ] = optimsimplex_new ( "randbounds" , x0.' , ..
985         [] , this.optbase.boundsmin , this.optbase.boundsmax , ..
986         this.boxnbpointseff  , this );
987     else
988         errmsg = msprintf(gettext("%s: Unknown value %s for -simplex0method option"), "neldermead_startup", this.simplex0method);
989         error(errmsg);
990     end
991     //
992     // 3. Scale the initial simplex into the bounds and the nonlinear inequality constraints, if any
993     [ this.optbase , hasnlcons ] = optimbase_hasnlcons ( this.optbase );
994     if ( hasbounds | hasnlcons ) then
995         // Check that initial guess is feasible
996         [ this.optbase , isfeasible ] = optimbase_isfeasible ( this.optbase , x0 );
997         if ( isfeasible <> 1 ) then
998             error ( msprintf ( gettext ( "%s: Initial guess [%s] is not feasible." ) , "neldermead_startup" , _strvec ( x0 ) ) )
999         end
1000         this = neldermead_log (this,sprintf("Scaling initial simplex into nonlinear inequality constraints..."));
1001         select this.scalingsimplex0
1002         case "tox0" then
1003             [ this , simplex0 ] = neldermead_scaletox0 ( this , simplex0 );
1004         case "tocenter" then
1005             [ this , simplex0 ] = neldermead_scaletocenter ( this , simplex0 );
1006         else
1007             errmsg = msprintf(gettext("%s: Unknown value %s for -scalingsimplex0 option"),"neldermead_startup", this.scalingsimplex0 );
1008             error(errmsg);
1009         end
1010     end
1011     //
1012     // 4. Compute the function values
1013     nbve = optimsimplex_getnbve ( simplex0 );
1014     for ive = 1 : nbve
1015         // Transpose, because optimsimplex returns row vectors
1016         x = optimsimplex_getx ( simplex0 , ive ).';
1017         if ( hasnlcons ) then
1018             [ this.optbase , fv , c , index ] = optimbase_function ( this.optbase , x , 2 );
1019         else
1020             [ this.optbase , fv , index ] = optimbase_function ( this.optbase , x , 2 );
1021         end
1022         // Transpose xp, which is a column vector
1023         simplex0 = optimsimplex_setve ( simplex0 , ive , fv , x.' );
1024     end
1025     //
1026     // 5. Store the simplex
1027     this.simplex0 = optimsimplex_destroy ( this.simplex0 );
1028     this.simplex0 = simplex0;
1029     this.simplexsize0 = optimsimplex_size ( simplex0 );
1030     //
1031     // 6. Store initial data into the base optimization component
1032     fx0 = optimsimplex_getfv ( this.simplex0 , 1 );
1033     this.optbase = optimbase_set ( this.optbase , "-fx0" , fx0 );
1034     this.optbase = optimbase_set ( this.optbase , "-xopt" , x0 );
1035     this.optbase = optimbase_set ( this.optbase , "-fopt" , fx0 );
1036     this.optbase = optimbase_set ( this.optbase , "-iterations" , 0 );
1037     //
1038     // 7. Initialize the termination criteria
1039     this = neldermead_termstartup ( this );
1040 endfunction
1041 //
1042 // neldermead_scaletox0 --
1043 //   Scale the simplex into the
1044 //   nonlinear inequality constraints, if any.
1045 //   Scale toward x0, which is feasible.
1046 //   Do not update the function values.
1047 // Arguments
1048 //   simplex0 : the initial simplex
1049 //
1050 function [ this , simplex0 ] = neldermead_scaletox0 ( this , simplex0 )
1051     [ this.optbase , hasnlcons ] = optimbase_hasnlcons ( this.optbase );
1052     nbve = optimsimplex_getnbve ( simplex0 );
1053     x0 = optimbase_cget ( this.optbase , "-x0" );
1054     for ive = 2 : nbve
1055         // Transpose, because optimsimplex returns row vectors
1056         x = optimsimplex_getx ( simplex0 , ive ).';
1057         this = neldermead_log (this,sprintf("Scaling vertex #%d/%d at ["+..
1058         _strvec(x)+"]... " , ive , nbve ));
1059         // Transpose x into a row vector
1060         [ this , status , xp ] = _scaleinconstraints ( this , x , x0 );
1061         if ( ~status ) then
1062             lclmsg = gettext("%s: Impossible to scale the vertex #%d/%d at [%s] into inequality constraints")
1063             errmsg = msprintf(lclmsg, ..
1064             "neldermead_startup", ive , nbve , _strvec(x));
1065             error(errmsg);
1066         end
1067         if ( or ( x <> xp ) ) then
1068             // Transpose xp, which is a column vector
1069             simplex0 = optimsimplex_setx ( simplex0 , ive , xp.' );
1070         end
1071     end
1072 endfunction
1073 //
1074 // neldermead_scaletocenter --
1075 //   Scale the simplex into the
1076 //   nonlinear inequality constraints, if any.
1077 //   Scale to the centroid of the points
1078 //   which satisfy the constraints.
1079 //   Do not update the function values.
1080 // Notes
1081 //   This is Box's method for scaling.
1082 //   It is unsure, since the centroid of the points
1083 //   which satisfy the constraints may not be feasible.
1084 // Arguments
1085 //
1086 //
1087 function [ this , simplex0 ] = neldermead_scaletocenter ( this , simplex0 , x0 )
1088     [ this.optbase , hasnlcons ] = optimbase_hasnlcons ( this.optbase );
1089     nbve = optimsimplex_getnbve ( simplex0 );
1090     for ive = 2 : nbve
1091         xref = optimsimplex_xbar ( simplex0 , ive:nbve ).';
1092         // Transpose, because optimsimplex returns row vectors
1093         x = optimsimplex_getx ( simplex0 , ive ).';
1094         this = neldermead_log (this,sprintf("Scaling vertex #%d/%d at ["+..
1095         _strvec(x)+"]... " , ive , nbve ));
1096         // Transpose x into a row vector
1097         [ this , status , xp ] = _scaleinconstraints ( this , x , xref );
1098         if ( ~status ) then
1099             errmsg = msprintf(gettext("%s: Impossible to scale the vertex #%d/%d at [%s] into inequality constraints"), ..
1100             "neldermead_startup", ive , nbve , _strvec(x));
1101             error(errmsg);
1102         end
1103         if ( or ( x <> xp ) ) then
1104             // Transpose xp, which is a column vector
1105             simplex0 = optimsimplex_setx ( simplex0 , ive , xp.' );
1106         end
1107     end
1108 endfunction
1109 //
1110 // neldermead_termstartup --
1111 //   Initialize Kelley's stagnation detection system when normalization is required,
1112 //   by computing kelleyalpha.
1113 //   If the simplex gradient is zero, then
1114 //   use alpha0 as alpha.
1115 // Arguments
1116 //   status : the status after the failing
1117 //     optimization process
1118 //   simplex : the simplex computed at the end of the failing
1119 //     optimization process
1120 //
1121 function this = neldermead_termstartup ( this )
1122     //
1123     // Criteria #8 : Kelley stagnation, based on
1124     // a sufficient decrease condition
1125     //
1126     if ( this.kelleystagnationflag ) then
1127         if ( ~this.kelleynormalizationflag ) then
1128             this.kelleyalpha = this.kelleystagnationalpha0
1129         else
1130             [sg,this] = optimsimplex_gradientfv ( this.simplex0 , neldermead_costf , "forward" , this )
1131             nsg = sg.' * sg
1132             sigma0 = optimsimplex_size ( this.simplex0 , "sigmaplus" )
1133             if nsg==0.0 then
1134                 this.kelleyalpha = this.kelleystagnationalpha0
1135             else
1136                 this.kelleyalpha = this.kelleystagnationalpha0 * sigma0 / nsg
1137             end
1138         end
1139         this = neldermead_log (this,sprintf("Test Stagnation Kelley : alpha0 = %s", string(this.kelleyalpha)));
1140     end
1141     //
1142     // Criteria #10 : variance of function values
1143     //
1144     if ( this.tolvarianceflag ) then
1145         this.variancesimplex0 = optimsimplex_fvvariance ( this.simplex0 )
1146     end
1147 endfunction
1148 //
1149 // _scaleinconstraints --
1150 //   Given a point to scale and a reference point which satisfies the constraints,
1151 //   scale the point towards the reference point until it satisfies all the constraints.
1152 //   Returns isscaled = %T if the procedure has succeded before -boxnbnlloops
1153 //   Returns isscaled = %F if the procedure has failed after -boxnbnlloops
1154 //   iterations.
1155 // Arguments
1156 //   x : the point to scale
1157 //   xref : the reference point
1158 //   isscaled : %T or %F
1159 //   p : scaled point
1160 //
1161 function [ this , isscaled , p ] = _scaleinconstraints ( this , x , xref )
1162     p = x
1163     [ this.optbase , hasbounds ] = optimbase_hasbounds ( this.optbase );
1164     nbnlc = optimbase_cget ( this.optbase , "-nbineqconst" )
1165     //
1166     // 1. No bounds, no nonlinear inequality constraints
1167     // => no problem
1168     //
1169     if ( ( hasbounds == %f ) & ( nbnlc == 0 ) ) then
1170         isscaled = %T
1171         return;
1172     end
1173     //
1174     // 2. Scale into bounds
1175     //
1176     if ( hasbounds ) then
1177         [ this.optbase , p ] = optimbase_proj2bnds ( this.optbase ,  p );
1178         this = neldermead_log (this,sprintf(" > After projection into bounds p = [%s]" , ..
1179         _strvec(p)));
1180     end
1181     //
1182     // 3. Scale into non linear constraints
1183     // Try the current point and see if the constraints are satisfied.
1184     // If not, move the point "halfway" to the centroid,
1185     // which should satisfy the constraints, if
1186     // the constraints are convex.
1187     // Perform this loop until the constraints are satisfied.
1188     // If all loops have been performed without success, the scaling
1189     // has failed.
1190     //
1191     isscaled = %F
1192     alpha = 1.0
1193     p0 = p
1194     while ( alpha > this.guinalphamin )
1195         [ this.optbase , feasible ] = optimbase_isinnonlincons ( this.optbase , p );
1196         if ( feasible ) then
1197             isscaled = %T;
1198             break;
1199         end
1200         alpha = alpha * this.boxineqscaling
1201         this = neldermead_log (this,sprintf("Scaling inequality constraint with alpha = %s", ..
1202         string(alpha)));
1203         p = ( 1.0 - alpha ) * xref + alpha * p0;
1204     end
1205     this = neldermead_log (this,sprintf(" > After scaling into inequality constraints p = [%s]" , ..
1206     _strvec(p) ) );
1207     if ( ~isscaled ) then
1208         this = neldermead_log (this,sprintf(" > Impossible to scale into constraints after %d loops" , ..
1209         this.optbase.nbineqconst ));
1210     end
1211 endfunction
1212 //
1213 // neldermead_logsummary --
1214 //   At each iteration, prints this iteration summary.
1215 //
1216 function this = neldermead_logsummary ( this, iter,xlow,flow,currentcenter,simplex )
1217     deltafv = abs(optimsimplex_deltafvmax ( simplex ));
1218     totaliter = optimbase_get ( this.optbase , "-iterations" );
1219     funevals = optimbase_get ( this.optbase , "-funevals" );
1220     ssize = optimsimplex_size ( simplex )
1221     this = neldermead_log (this,sprintf("================================================================="));
1222     this = neldermead_log (this,sprintf("Iteration #%d (total = %d)",iter,totaliter));
1223     this = neldermead_log (this,sprintf("Function Eval #%d",funevals));
1224     this = neldermead_log (this,sprintf("Xopt : [%s]",_strvec(xlow)));
1225     this = neldermead_log (this,sprintf("Fopt : %s",string(flow)));
1226     this = neldermead_log (this,sprintf("DeltaFv : %s",string(deltafv)));
1227     this = neldermead_log (this,sprintf("Center : [%s]",_strvec(currentcenter)));
1228     this = neldermead_log (this,sprintf("Size : %s",string(ssize)));
1229     str = string ( simplex )
1230     for i = 1:size(str,"*")
1231         this = neldermead_log (this,str(i));
1232     end
1233 endfunction
1234
1235 //
1236 // neldermead_box --
1237 //   The Nelder-Mead algorithm, with variable-size simplex
1238 //   and modifications by Box for bounds and
1239 //   inequality constraints.
1240 //
1241 function this = neldermead_box ( this )
1242     // Check settings correspond to algo
1243     [ this.optbase , hascons ] = optimbase_hasconstraints ( this.optbase );
1244     if ( ~hascons ) then
1245         errmsg = msprintf(gettext("%s: Problem has no constraints, but Box algorithm is designed for them."), "neldermead_box")
1246         error(errmsg)
1247     end
1248     verbose = optimbase_cget ( this.optbase , "-verbose" )
1249     //
1250     // Order the vertices for the first time
1251     //
1252     simplex = this.simplex0;
1253     n = optimbase_cget ( this.optbase , "-numberofvariables" );
1254     fvinitial = optimbase_get ( this.optbase , "-fx0" );
1255     // Sort function values and x points by increasing function value order
1256     this = neldermead_log (this,"Step #1 : order");
1257     simplex = optimsimplex_sort ( simplex );
1258     // Transpose, because optimsimplex returns row vectors
1259     currentcenter = optimsimplex_center ( simplex ).';
1260     currentxopt = optimbase_cget ( this.optbase , "-x0" );
1261     newfvmean = optimsimplex_fvmean ( simplex );
1262     nbve = optimsimplex_getnbve ( simplex );
1263     ihigh = nbve;
1264     inext = ihigh - 1
1265     ilow = 1
1266     [ this.optbase , hasbounds ] = optimbase_hasbounds ( this.optbase );
1267     nbnlc = optimbase_cget ( this.optbase , "-nbineqconst" )
1268     rho = this.boxreflect;
1269     //
1270     // Initialize
1271     //
1272     terminate = %f;
1273     iter = 0;
1274     step = "init";
1275     //
1276     // Nelder-Mead Loop
1277     //
1278     while ( ~terminate )
1279         this.optbase = optimbase_incriter ( this.optbase );
1280         iter = iter + 1;
1281         xlow = optimsimplex_getx ( simplex , ilow ).'
1282         flow = optimsimplex_getfv ( simplex , ilow )
1283         xhigh = optimsimplex_getx ( simplex , ihigh ).'
1284         fhigh = optimsimplex_getfv ( simplex , ihigh )
1285         xn = optimsimplex_getx ( simplex , inext ).'
1286         fn = optimsimplex_getfv ( simplex , inext )
1287         //
1288         // Store history
1289         //
1290         xcoords = optimsimplex_getallx ( simplex )
1291         this = neldermead_storehistory ( this , n , flow , xlow , xcoords );
1292         currentfopt = flow;
1293         previousxopt = currentxopt;
1294         currentxopt = xlow;
1295         previouscenter = currentcenter;
1296         currentcenter = optimsimplex_center ( simplex ).';
1297         oldfvmean = newfvmean;
1298         newfvmean = optimsimplex_fvmean ( simplex );
1299         if ( verbose == 1 ) then
1300             this = neldermead_logsummary ( this, iter,xlow,flow,currentcenter,simplex);
1301         end
1302         this.optbase = optimbase_set ( this.optbase , "-xopt" , xlow );
1303         this.optbase = optimbase_set ( this.optbase , "-fopt" , flow );
1304         stop = neldermead_outputcmd ( this, "iter" , simplex , step )
1305         if ( stop ) then
1306             status = "userstop"
1307             if ( verbose == 1 ) then
1308                 this = neldermead_log (this,sprintf("Terminate by user''s request"));
1309             end
1310             break
1311         end
1312
1313         //
1314         // Update termination flag
1315         //
1316         if ( iter > 1 ) then
1317             [ this , terminate , status ] = neldermead_termination (this , ..
1318             fvinitial , oldfvmean , newfvmean , previouscenter , currentcenter , simplex );
1319             if ( terminate ) then
1320                 if ( verbose == 1 ) then
1321                     this = neldermead_log (this,sprintf("Terminate with status : %s",status));
1322                 end
1323                 break
1324             end
1325         end
1326         //
1327         // Compute xbar, center of better vertices
1328         //
1329         if ( verbose == 1 ) then
1330             this = neldermead_log (this,sprintf("Reflect"));
1331         end
1332         xbar = optimsimplex_xbar ( simplex ).';
1333         if ( verbose == 1 ) then
1334             this = neldermead_log (this,sprintf("xbar=[%s]" , _strvec(xbar)));
1335         end
1336         //
1337         // Search a point improving cost function
1338         // and satisfying constraints.
1339         //
1340         [ this , scaled , xr , fr ] = _boxlinesearch ( this , n , xbar , xhigh , fhigh , rho );
1341         if ( scaled == %f ) then
1342             status = "impossibleimprovement"
1343             break
1344         end
1345         if ( verbose == 1 ) then
1346             this = neldermead_log (this,sprintf("xr=[%s], f(xr)=%f", _strvec(xr) , fr));
1347             this = neldermead_log (this,sprintf("  > Perform Reflection"));
1348         end
1349         simplex = optimsimplex_setve ( simplex , ihigh , fr , xr.' )
1350         step = "boxreflection";
1351         //
1352         // Sort simplex
1353         //
1354         if ( verbose == 1 ) then
1355             this = neldermead_log (this,sprintf("Sort"));
1356         end
1357         simplex  = optimsimplex_sort ( simplex );
1358     end
1359     this.optbase = optimbase_set ( this.optbase , "-xopt" , xlow );
1360     this.optbase = optimbase_set ( this.optbase , "-fopt" , flow );
1361     this.optbase = optimbase_set ( this.optbase , "-status" , status );
1362     this.simplexopt = simplex;
1363 endfunction
1364
1365 //
1366 // _strvec --
1367 //  Returns a string for the given vector.
1368 //
1369 function str = _strvec ( x )
1370     if isempty(x) then
1371         str = "";
1372         return
1373     end
1374     str = strcat(string(x)," ")
1375 endfunction
1376 //
1377 // _boxlinesearch --
1378 //   For Box's method, perform a line search
1379 //   from xbar, on the line (xhigh,xbar) and returns:
1380 //   status : %t if the search is successful,
1381 //            status=%f if the search failed.
1382 //   xr : the reflected, scaled point
1383 //   fr : the function value. fr=f(xr) if the search successful,
1384 //        fr is Nan if the search failed.
1385 //   The reflected point satisfies the following
1386 //   constraints :
1387 //   * fr < fhigh
1388 //   * xr satisfies the bounds constraints
1389 //   * xr satisfies the nonlinear positive inequality constraints
1390 //   * xr satisfies the linear positive inequality constraints
1391 //   The method is based on projection and scaling toward the centroid xbar.
1392 //
1393 // Arguments
1394 //   n : number of variables
1395 //   xbar : the centroid
1396 //   xhigh : the high point
1397 //   fhigh : function value at xhigh
1398 //   rho : reflection factor
1399 //
1400 function [ this , status , xr , fr ] = _boxlinesearch ( this , n , xbar , xhigh , fhigh , rho )
1401     if ( verbose == 1 ) then
1402         this = neldermead_log (this,"_boxlinesearch");
1403         this = neldermead_log (this, sprintf ("> xhigh=[%s], fhigh=%s",_strvec(xhigh),string(fhigh)));
1404         this = neldermead_log (this, sprintf ("> xbar=[%s]" , _strvec(xbar) ) );
1405     end
1406     xr = neldermead_interpolate ( xbar , xhigh , rho );
1407     if ( verbose == 1 ) then
1408         this = neldermead_log (this, sprintf ( "> xr = [%s]" , _strvec ( xr ) ) );
1409     end
1410     alphamin = this.guinalphamin
1411     [ this.optbase , hasnlcons ] = optimbase_hasnlcons ( this.optbase );
1412     //
1413     // The algorithm has 3 steps:
1414     // 1. scale for bounds (cannot fail)
1415     // 2. scale for nonlinear constraints (may fail)
1416     // 3. scale for function improvement (may fail)
1417     //
1418     // 1. Project xr into bounds, with an additional alpha inside the bounds.
1419     // This algo is always succesful.
1420     // Note:
1421     //   If the alpha coefficient was not used, the
1422     //   projectinbounds method could be used directly.
1423     //
1424     [ this.optbase , hasbounds ] = optimbase_hasbounds ( this.optbase );
1425     if ( hasbounds ) then
1426         boxboundsalpha = this.boxboundsalpha;
1427         boundsmax = optimbase_cget ( this.optbase , "-boundsmax" );
1428         boundsmin = optimbase_cget ( this.optbase , "-boundsmin" );
1429         for ix = 1:n
1430             xmin = boundsmin ( ix );
1431             xmax = boundsmax ( ix );
1432             xrix = xr ( ix );
1433             if ( xrix > xmax ) then
1434                 xr ( ix ) = xmax - boxboundsalpha;
1435             elseif ( xrix < xmin ) then
1436                 xr ( ix ) = xmin + boxboundsalpha;
1437             end
1438         end
1439     end
1440     //
1441     // 2. Scale from xr to xbar into nonlinear inequality constraints
1442     // and update xr.
1443     // Set status to %f if the process fails, set status=%t if it succeeds.
1444     //
1445     nbnlc = optimbase_cget ( this.optbase , "-nbineqconst" );
1446     if ( nbnlc == 0 ) then
1447         status = %t
1448     else
1449         status = %f;
1450         alpha = 1.0;
1451         xr0 = xr;
1452         while ( alpha > alphamin )
1453             [ this.optbase , feasible ] = optimbase_isinnonlincons ( this.optbase , xr );
1454             if ( feasible ) then
1455                 status = %t;
1456                 break
1457             end
1458             alpha = alpha * this.boxineqscaling;
1459             xr = ( 1.0 - alpha ) * xbar + alpha * xr0;
1460         end
1461     end
1462     // If scaling failed, returns immediately
1463     // (there is no need to update the function value).
1464     if ( ~status ) then
1465         fr = %nan
1466         return
1467     end
1468     //
1469     // 3. Scale from xr toward xbar until fr < fhigh.
1470     //
1471     status = %f;
1472     xr0 = xr
1473     alpha = 1.0
1474     while ( alpha > alphamin )
1475         if ( hasnlcons ) then
1476             [ this.optbase , fr , cr , index ] = optimbase_function ( this.optbase , xr , 2 );
1477         else
1478             [ this.optbase , fr , index ] = optimbase_function ( this.optbase , xr , 2 );
1479         end
1480         if ( fr < fhigh ) then
1481             status = %t;
1482             break
1483         end
1484         alpha = alpha * this.boxineqscaling;
1485         xr = ( 1.0 - alpha ) * xbar + alpha * xr0;
1486     end
1487     // If the scaling for function improvement has failed,
1488     // we return.
1489     if ( ~status ) then
1490         fr = %nan
1491         return;
1492     end
1493 endfunction
1494 //
1495 // costf_transposex --
1496 //   Call the cost function and return the value.
1497 //   Transpose the value of x, so that the input row vector,
1498 //   given by optimsimplex, is transposed into a column vector,
1499 //   as required by the cost function.
1500 //
1501 function [ f , this ] = costf_transposex ( x , this )
1502     xt = x.'
1503     [ f , this ] = neldermead_costf ( xt , this )
1504 endfunction