Added logging to file mode. Fixed Box algorithm and separated from nmconstraints.
Michaƫl Baudin [Fri, 2 Oct 2009 08:55:31 +0000 (10:55 +0200)]
31 files changed:
scilab/modules/optimization/demos/neldermead/neldermead_boxproblemA.sce
scilab/modules/optimization/demos/neldermead/neldermead_boxproblemB.sce [new file with mode: 0644]
scilab/modules/optimization/help/en_US/neldermead/neldermead.xml
scilab/modules/optimization/help/en_US/optimbase/optimbase.xml
scilab/modules/optimization/macros/neldermead/neldermead_cget.sci
scilab/modules/optimization/macros/neldermead/neldermead_configure.sci
scilab/modules/optimization/macros/neldermead/neldermead_new.sci
scilab/modules/optimization/macros/neldermead/neldermead_search.sci
scilab/modules/optimization/macros/neldermead/neldermead_updatesimp.sci
scilab/modules/optimization/macros/optimbase/optimbase_cget.sci
scilab/modules/optimization/macros/optimbase/optimbase_checkx0.sci
scilab/modules/optimization/macros/optimbase/optimbase_configure.sci
scilab/modules/optimization/macros/optimbase/optimbase_destroy.sci
scilab/modules/optimization/macros/optimbase/optimbase_function.sci
scilab/modules/optimization/macros/optimbase/optimbase_isfeasible.sci
scilab/modules/optimization/macros/optimbase/optimbase_isinbounds.sci
scilab/modules/optimization/macros/optimbase/optimbase_isinnonlincons.sci
scilab/modules/optimization/macros/optimbase/optimbase_log.sci
scilab/modules/optimization/macros/optimbase/optimbase_logshutdown.sci [new file with mode: 0644]
scilab/modules/optimization/macros/optimbase/optimbase_logstartup.sci [new file with mode: 0644]
scilab/modules/optimization/macros/optimbase/optimbase_new.sci
scilab/modules/optimization/macros/optimbase/optimbase_proj2bnds.sci
scilab/modules/optimization/macros/optimbase/optimbase_stoplog.sci
scilab/modules/optimization/macros/optimbase/optimbase_terminate.sci
scilab/modules/optimization/tests/unit_tests/neldermead/neldermead_rosensuzuki.dia.ref
scilab/modules/optimization/tests/unit_tests/neldermead/neldermead_rosensuzuki.tst
scilab/modules/optimization/tests/unit_tests/neldermead/neldermead_search.dia.ref
scilab/modules/optimization/tests/unit_tests/neldermead/neldermead_search.tst
scilab/modules/optimization/tests/unit_tests/neldermead/neldermead_searchaddarg2.tst
scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_function.tst
scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_log.tst [new file with mode: 0644]

index e047983..d70595f 100644 (file)
 //   numerical experiment presented in Box's paper.
 //
 
-mprintf("Illustrates Box''' algorithm on Box problem A.\n");
+mprintf("Illustrates Box'' algorithm on Box problem A.\n");
 mprintf("Defining Box Problem A function...\n");
 
 //
+//  Reference:
+//
+//   M.J.Box, 
+//   "A new method of constrained optimization 
+//   and a comparison with other methods".
+//   The Computer Journal, Volume 8, Number 1, 1965, 42--52
+//   Problem A
+//
+//   Algorithm 454: the complex method for constrained optimization [E4]
+//   Communications of the ACM, Volume 16 ,  Issue 8  (August 1973)
+//   Pages: 487 - 489   
+//
+
+//
 // boxproblemA --
 //   Computes the Box problem A cost function and 
 //   inequality constraints.
@@ -37,7 +51,7 @@ mprintf("Defining Box Problem A function...\n");
 //    constraints
 //  The inequality constraints are expected to be positive.
 //
-function result = boxproblemA ( x , index , data )
+function [ result , data ] = boxproblemA ( x , index , data )
   if (~isdef('index','local')) then
     index = 1
   end
@@ -130,18 +144,21 @@ boxparams.k33 = 23.3088196;
 boxparams.k34 = -27097.648;
 boxparams.k35 = -50843.766;
 
-    //
-    // Initialize the random number generator, so that the results are always the
-    // same.
-    //
-    rand("seed" , 0)
+//
+// Initialize the random number generator, so that the results are always the
+// same.
+//
+rand("seed" , 0)
 
 x0 = [2.52 2.0 37.5 9.25 6.8].';
 // Compute f(x0) : should be close to -2351244.0
 fx0 = boxproblemA ( x0 , data = boxparams );
+mprintf("Computed fx0 = %e (expected = %e)\n",fx0 , -2351244. );
+
 xopt = [4.53743 2.4 60.0 9.3 7.0].';
 // Compute f(xopt) : should be -5280334.0
 fopt = boxproblemA ( xopt , data = boxparams );
+mprintf("Computed fopt = %e (expected = %e)\n", fopt , -5280334.0 );
 
 nm = neldermead_new ();
 nm = neldermead_configure(nm,"-numberofvariables",5);
@@ -152,7 +169,8 @@ nm = neldermead_configure(nm,"-maxiter",300);
 nm = neldermead_configure(nm,"-maxfunevals",1000);
 nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-4);
 nm = neldermead_configure(nm,"-method","box");
-//nm = neldermead_configure(nm,"-verbose",1);
+nm = neldermead_configure(nm,"-verbose",1);
+nm = neldermead_configure(nm,"-logfile" , "boxproblemA.txt" );
 nm = neldermead_configure(nm,"-verbosetermination",1);
 // Configure like Box
 nm = neldermead_configure(nm,"-boundsmin",[0.0 1.2 20.0 9.0 6.0]);
@@ -181,5 +199,5 @@ mprintf("f expected=%f\n",fopt);
 shift = abs(fcomp-fopt)/abs(fopt);
 mprintf("Shift =%f\n",shift);
 nm = neldermead_destroy(nm);
-
+deletefile ( "boxproblemA.txt" )
 
diff --git a/scilab/modules/optimization/demos/neldermead/neldermead_boxproblemB.sce b/scilab/modules/optimization/demos/neldermead/neldermead_boxproblemB.sce
new file mode 100644 (file)
index 0000000..e7111d0
--- /dev/null
@@ -0,0 +1,148 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008-2009 - INRIA - Michael Baudin
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+
+//
+// nmplot_boxproblemA.sce --
+//   Show that the Box algorithm is able to reproduce the 
+//   numerical experiment presented in Box's paper.
+//
+
+mprintf("Illustrates Box'' algorithm on Box problem A.\n");
+mprintf("Defining Box Problem A function...\n");
+
+//
+//  Reference:
+//
+//   M.J.Box, 
+//   "A new method of constrained optimization 
+//   and a comparison with other methods".
+//   The Computer Journal, Volume 8, Number 1, 1965, 42--52
+//   Problem A
+//
+//   Algorithm 454: the complex method for constrained optimization [E4]
+//   Communications of the ACM, Volume 16 ,  Issue 8  (August 1973)
+//   Pages: 487 - 489   
+//
+// Note 
+//    The maximum bound for the parameter x1
+//    is not indicated in Box's paper, but indicated in "Algo 454".
+//    The maximum bound for x2 is set to 100/sqrt(3) and satisfies the constraint on x2.
+//    The original problem was to maximize the cost function.
+//    Here, the problem is to minimize the cost function.
+
+//
+// boxproblemB --
+//   Computes the Box problem B cost function and 
+//   inequality constraints.
+//
+// Arguments
+//   x: the point where to compute the function
+//   index : the stuff to compute
+//   data : the parameters of Box cost function
+// Note
+//  The following protocol is used
+//  * if index=1, or no index, returns the value of the cost 
+//    function (default case)
+//  * if index=2, returns the value of the nonlinear inequality 
+//    constraints, as a row array
+//  * if index=3, returns an array which contains
+//    at index #0, the value of the cost function  
+//    at index #1 to the end is the list of the values of the nonlinear 
+//    constraints
+//  The inequality constraints are expected to be positive.
+//
+function result = boxproblemB ( x , index )
+  if (~isdef('index','local')) then
+    index = 1
+  end
+  x3 = x(1) + sqrt(3.0) * x(2)
+  if ( index==1 | index==3 ) then
+    f = -(9.0 - (x(1) - 3.0) ^ 2) * x(2) ^ 3 / 27.0 / sqrt(3.0)
+  end
+  if ( index==2 | index==3 ) then
+      c1 = x(1) / sqrt(3.0) - x(2)
+      c2 = x3
+      c3 = 6.0 - x3
+  end
+  select index
+  case 1 then
+    result = f
+  case 2 then
+    result = [c1 c2 c3]
+  case 3 then
+    result = [f c1 c2 c3]
+  else
+    errmsg = sprintf("Unexpected index %d" , index);
+    error(errmsg);
+  end
+endfunction
+
+
+//
+// Initialize the random number generator, so that the results are always the
+// same.
+//
+rand("seed" , 0)
+
+x0 = [1.0 0.5].';
+// Compute f(x0) : should be close to -0.0133645895646
+fx0 = boxproblemB ( x0 );
+mprintf("Computed fx0 = %e (expected = %e)\n",fx0 , -0.0133645895646 );
+result = boxproblemB ( x0 , 2 );
+mprintf("Computed Constraints(x0) = [%e %e %e]\n", ...
+  result(1), result(2), result(3) );
+mprintf("Expected Constraints(x0) = [%e %e %e]\n", ...
+  0.0773503 , 1.8660254 , 4.1339746 );
+
+
+xopt = [3.0 1.7320508075688774].';
+// Compute f(xopt) : should be -1.0
+fopt = boxproblemB ( xopt );
+mprintf("Computed fopt = %e (expected = %e)\n", fopt , -1.0 );
+
+nm = neldermead_new ();
+nm = neldermead_configure(nm,"-numberofvariables",2);
+nm = neldermead_configure(nm,"-function",boxproblemB);
+nm = neldermead_configure(nm,"-x0",x0);
+nm = neldermead_configure(nm,"-maxiter",300);
+nm = neldermead_configure(nm,"-maxfunevals",300);
+nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-4);
+nm = neldermead_configure(nm,"-method","box");
+nm = neldermead_configure(nm,"-verbose",1);
+nm = neldermead_configure(nm,"-logfile" , "boxproblemB.txt" );
+nm = neldermead_configure(nm,"-verbosetermination",1);
+// Configure like Box
+nm = neldermead_configure(nm,"-boundsmin",[0.0 0.0]);
+nm = neldermead_configure(nm,"-boundsmax",[100.0 57.735026918962582]);
+nm = neldermead_configure(nm,"-simplex0method","randbounds");
+nm = neldermead_configure(nm,"-nbineqconst",6);
+//
+// Check that the cost function is correctly connected.
+// The index must be provided, because the additionnal argument "data"
+// comes after.
+//
+[ nm , result ] = neldermead_function ( nm , x0 , 1 );
+//
+// Perform optimization
+//
+nm = neldermead_search(nm);
+neldermead_display(nm);
+xcomp = neldermead_get(nm,"-xopt");
+mprintf("x computed=%s\n",strcat(string(xcomp)," "));
+mprintf("x expected=%s\n",strcat(string(xopt)," "));
+shift = norm(xcomp-xopt)/norm(xopt);
+mprintf("Shift =%f\n",shift);
+fcomp = neldermead_get(nm,"-fopt");
+mprintf("f computed=%f\n",fcomp);
+mprintf("f expected=%f\n",fopt);
+shift = abs(fcomp-fopt)/abs(fopt);
+mprintf("Shift =%f\n",shift);
+nm = neldermead_destroy(nm);
+deletefile ( "boxproblemB.txt" )
+
index d1565b8..aa28b47 100644 (file)
@@ -1991,7 +1991,7 @@ Iteration #157, Feval #300, Fval = 6.023002e-027, x = 1 1, Size = 2.814328e-013
 
     <itemizedlist>
       <listitem>
-        <para>add the Box-Guin algoritm as a 4th method.</para>
+        <para>Add the Box-Guin algoritm as a 4th method.</para>
 
         <para>This algorithm solves an constrained optimization problem with a
         variable sized simplex made of an arbitrary k number of vertices. This
@@ -1999,8 +1999,21 @@ Iteration #157, Feval #300, Fval = 6.023002e-027, x = 1 1, Size = 2.814328e-013
       </listitem>
 
       <listitem>
-        <para>add the optimization of the Rosenbrock test case, with the
-        interactive plot thanks to the -outputcommand option</para>
+        <para>Add the optimization of the Rosenbrock test case, with the
+        interactive plot thanks to the -outputcommand option.</para>
+      </listitem>
+
+      <listitem>
+        <para>Add an option to let the user setup his own termination
+        criteria.</para>
+
+        <para>This might be handy in case where the user try to reproduce the
+        behaviour of another algorithm (be it a software or a paper), where
+        none of the provided termination criteria match the required one. This
+        feature would allow to use our own termination criteria, without
+        modifying the internal source code of the component. This might be
+        dangerous if the user does not clearly understand what is going on.
+        But that might be very handy as well.</para>
       </listitem>
     </itemizedlist>
   </refsection>
index 5bf1cf6..28c77f0 100644 (file)
@@ -1112,7 +1112,7 @@ g_i(x) &lt;= 0, i = 1,nbineq
               <term>terminate</term>
 
               <listitem>
-                <para>1 if the algorithm must terminate, 0 if the algorithm
+                <para>%t if the algorithm must terminate, %f if the algorithm
                 must continue</para>
               </listitem>
             </varlistentry>
@@ -1121,8 +1121,8 @@ g_i(x) &lt;= 0, i = 1,nbineq
               <term>status</term>
 
               <listitem>
-                <para>if terminate = 1, the detailed status of the
-                termination, as a string. If terminate = 0, the status is
+                <para>if terminate = %t, the detailed status of the
+                termination, as a string. If terminate = %f, the status is
                 "continue".</para>
 
                 <para>The following status are available :</para>
@@ -1179,8 +1179,8 @@ g_i(x) &lt;= 0, i = 1,nbineq
           <para>Check that the cost function is correctly connected. Generate
           an error if there is one. Takes into account for the cost function
           at the initial guess x0 only. If there are nonlinear constraints,
-          check that the index flag is correctly used for index=1,2 or 3.
-          </para>
+          check that the index flag is correctly used for index=1,2 or
+          3.</para>
 
           <para>This function requires at least one call to the cost function
           to make the necessary checks.</para>
index 60d8803..8560a69 100644 (file)
@@ -69,14 +69,22 @@ function value = neldermead_cget (this,key)
     value = this.restartsimplexmethod;
   case "-boxnbpoints" then
     value = this.boxnbpoints;
-  case "-nbineqloops" then
-    value = this.nbineqloops;
   case "-ineqscaling" then
     value = this.ineqscaling;
   case "-checkcostfunction" then
     value = this.checkcostfunction;
   case "-scalingmethod" then
     value = this.scalingmethod;
+  case "-guinalphamin" then
+    value = this.guinalphamin;
+  case "-boxtermination" then
+    value = this.boxtermination
+  case "-boxtolf" then
+    value = this.boxtolf
+  case "-boxnbmatch" then
+    value = this.boxnbmatch
+  case "-boxreflect" then
+    value = this.boxreflect
   else
     // Delegate to the optimization object
     value = optimbase_cget ( this.optbase , key );
index 2a968a2..2cea51f 100644 (file)
@@ -21,6 +21,8 @@ function this = neldermead_configure (this,key,value)
       this.method = "variable";
     case "box" then
       this.method = "box";
+    case "nmconstraints" then
+      this.method = "nmconstraints";
     else
       errmsg = msprintf(gettext("%s: Unknown value %s for -method option"),"neldermead_configure",value);
       error(errmsg);
@@ -115,8 +117,6 @@ function this = neldermead_configure (this,key,value)
     this.restartsimplexmethod = value;
   case "-boxnbpoints" then
     this.boxnbpoints = value;
-  case "-nbineqloops" then
-    this.nbineqloops = value;
   case "-ineqscaling" then
     this.ineqscaling = value;
   case "-checkcostfunction" then
@@ -131,6 +131,22 @@ function this = neldermead_configure (this,key,value)
     end
   case "-scalingmethod" then
     this.scalingmethod = value;
+  case "-guinalphamin" then
+    if ( value <=0.0 ) then 
+      errmsg = msprintf(gettext("%s: Unexpected negative value %s for -guinalphamin option"),"neldermead_configure", value);
+      error(errmsg);
+    end
+    this.guinalphamin = value;
+  case "-boxboundsalpha" then
+    this.boxboundsalpha = value
+  case "-boxtermination" then
+    this.boxtermination = value
+  case "-boxtolf" then
+    this.boxtolf = value
+  case "-boxnbmatch" then
+    this.boxnbmatch = value
+  case "-boxreflect" then
+    this.boxreflect = value
   else
     // Delegate to the optimization object
     this.optbase = optimbase_configure ( this.optbase , key , value );
index cbacc96..6fcb2f4 100644 (file)
 //   Creates a new Nelder-Mead object.
 //
 function newobj = neldermead_new ()
-  newobj = tlist(["T_NELDERMEAD",...
-    "optbase","method",...
-    "simplex0","simplex0method","simplex0length",...
-    "rho","chi","gamma","sigma",...
-    "tolfstdeviation","tolfstdeviationmethod",...
-    "tolsimplexizeabsolute","tolsimplexizerelative","tolsimplexizemethod", "simplexsize0", ...
-    "toldeltafv","tolssizedeltafvmethod",...
-    "historysimplex", ...
-    "coords0",...
-    "simplex0deltausual","simplex0deltazero", ...
-    "restartsimplexmethod",...
-    "simplexopt","restartmax" , "restarteps", ...
-    "restartstep","kelleystagnationflag",...
-    "kelleynormalizationflag","kelleystagnationalpha0", ...
-    "kelleyalpha","restartnb","restartflag","restartdetection" , ...
-    "startupflag" , ...
-    "boxnbpoints" , "boxnbpointseff" , "nbineqloops" , "ineqscaling" , ...
-    "checkcostfunction" , "scalingmethod" ]);
+  newobj = tlist(["T_NELDERMEAD" 
+    "optbase" 
+    "method"
+    "simplex0"
+    "simplex0method"
+    "simplex0length"
+    "rho"
+    "chi"
+    "gamma"
+    "sigma"
+    "tolfstdeviation"
+    "tolfstdeviationmethod"
+    "tolsimplexizeabsolute"
+    "tolsimplexizerelative"
+    "tolsimplexizemethod"
+    "simplexsize0"
+    "toldeltafv"
+    "tolssizedeltafvmethod"
+    "historysimplex"
+    "coords0"
+    "simplex0deltausual"
+    "simplex0deltazero"
+    "restartsimplexmethod"
+    "simplexopt"
+    "restartmax"
+    "restarteps"
+    "restartstep"
+    "kelleystagnationflag"
+    "kelleynormalizationflag"
+    "kelleystagnationalpha0"
+    "kelleyalpha"
+    "restartnb"
+    "restartflag"
+    "restartdetection"
+    "startupflag"
+    "boxnbpoints"
+    "boxnbpointseff"
+    "nbineqloops"
+    "ineqscaling"
+    "checkcostfunction"
+    "scalingmethod"
+    "guinalphamin"
+    "boxboundsalpha"
+    "boxtermination"
+    "boxtolf"
+    "boxnbmatch"
+    "boxkount"
+    "boxreflect"
+    ]);
+
   newobj.optbase = optimbase_new();
   // Possible values "variable", "fixed".
   newobj.method = "variable";
@@ -70,7 +102,7 @@ function newobj = neldermead_new ()
   newobj.coords0 = [];
   // The Kelley stagnation detection in termination criteria :  0/1
   // (i.e. sufficient decrease of function value)
-  newobj.kelleystagnationflag = 0
+  newobj.kelleystagnationflag = 0// TODO : move this to %f/%t
   // The Kelley stagnation detection parameter
   newobj.kelleystagnationalpha0 = 1.e-4
   // The Kelley stagnation detection can be normalized or not.
@@ -95,22 +127,19 @@ function newobj = neldermead_new ()
   // Possible values : "oriented", "axes", "spendley", "pfeffer" 
   newobj.restartsimplexmethod = "oriented";
   // Possible values : 0/1
-  newobj.restartflag = 0;
+  newobj.restartflag = 0;// TODO : move this to %f/%t
   // Number of restarts performed
   newobj.restartnb = 0;
   // Type of restart detection method : "kelley", "oneill"
   newobj.restartdetection = "oneill";
   // Set to 1 when the startup has been performed
-  newobj.startupflag = 0;
+  newobj.startupflag = 0; // TODO : move this to %f/%t
   // Initial size of the simplex, for the tolerance on the simplex size
   newobj.simplexsize0 = 0.0
   // Number of points required in the simplex (for Box method)
   newobj.boxnbpoints = "2n"
   // Effective number of points required in the simplex (for Box method)
   newobj.boxnbpointseff = 0
-  // Number of loops performed to satisfy nonlinear
-  // inequality constraints (for Box method)
-  newobj.nbineqloops = 10
   // The scaling coefficient in nonlinear inequality constraints
   // in Box method, in (0,1) range
   newobj.ineqscaling = 0.5
@@ -118,4 +147,23 @@ function newobj = neldermead_new ()
   newobj.checkcostfunction = 1;
   // The scaling algorithm : "tox0", "tocentroid"
   newobj.scalingmethod = "tox0";
+  // Minimum alpha for constraints scaling
+  newobj.guinalphamin = 1.e-5;
+  // Box's alpha coefficient for bounds constraints.
+  // The value used in Box's paper was 0.000001 (delta in
+  // Richardson and Kuester's algorithm 454)
+  newobj.boxboundsalpha = 0.000001
+  // Set to 1 to enable Box termination criteria
+  newobj.boxtermination = 0
+  // The absolute tolerance on function value in Box termination criteria (beta in
+  // Richardson and Kuester's algorithm 454)
+  newobj.boxtolf = 1.e-5
+  // The number of consecutive match in Box termination criteria (gamma in
+  // Richardson and Kuester's algorithm 454)
+  newobj.boxnbmatch = 5
+  // Current number of consecutive match
+  newobj.boxkount = 0
+  // Box reflection/expansion factor
+  newobj.boxreflect = 1.3
 endfunction
+
index b5b8a64..15ffcba 100644 (file)
@@ -61,6 +61,8 @@ function this = neldermead_algo ( this )
       this = neldermead_variable (this);
     case "box" then
       this = neldermead_box (this);
+    case "nmconstraints" then
+      this = neldermead_constraints ( this );
     else
       errmsg = msprintf(gettext("%s: Unknown -method %s"), "neldermead_algo", this.method)
       error(errmsg)
@@ -170,7 +172,7 @@ function this = neldermead_variable ( this )
     if ( iter > 1 ) then
       [this , terminate , status] = neldermead_termination (this , ...
         fvinitial , oldfvmean , newfvmean , previouscenter , currentcenter , simplex );
-      if (terminate==1) then
+      if ( terminate ) then
         this = neldermead_log (this,sprintf("Terminate with status : %s",status));
         break
       end
@@ -355,7 +357,7 @@ function this = neldermead_fixed (this)
     if ( iter > 1 ) then
       [this , terminate , status] = neldermead_termination (this , ...
         fvinitial , oldfvmean , newfvmean , previouscenter , currentcenter , simplex );
-      if (terminate==1) then
+      if ( terminate ) then
         this = neldermead_log (this,sprintf("Terminate with status : %s",status));
         break;
       end
@@ -431,7 +433,7 @@ endfunction
 //   simplex : the simplex
 //     The best point in the simplex is expected to be stored at 1
 //     The worst point in the simplex is expected to be stored at n+1
-//   terminate : 1 if the algorithm terminates, 0 if the algorithm must continue.
+//   terminate : %t if the algorithm terminates, %f if the algorithm must continue.
 //   this.status : termination status
 //     status = "continue"
 //     status = "maxiter"
@@ -451,7 +453,7 @@ endfunction
 function [this , terminate , status ] = neldermead_termination (this , ...
   fvinitial , oldfvmean , newfvmean , previousxopt , currentxopt , ...
   simplex )
-  terminate = 0;
+  terminate = %f;
   status = "continue";
   //
   // Termination Criteria from parent optimization class
@@ -461,14 +463,14 @@ function [this , terminate , status ] = neldermead_termination (this , ...
   //
   // Criteria #5 : standard deviation of function values
   //
-  if (terminate == 0) then
+  if ( ~terminate ) then
     if this.tolfstdeviationmethod == "enabled" then
       fv = optimsimplex_getallfv ( simplex )
       sd = st_deviation(fv);
       this.optbase = optimbase_stoplog  ( this.optbase,sprintf("  > st_deviation(fv)=%e < tolfstdeviation=%e",...
         sd, this.tolfstdeviation));
       if sd < this.tolfstdeviation then
-        terminate = 1;
+        terminate = %f;
         status = "tolfstdev";
       end
     end
@@ -476,13 +478,13 @@ function [this , terminate , status ] = neldermead_termination (this , ...
   //
   // Criteria #6 : simplex absolute + relative size
   //
-  if (terminate == 0) then
+  if ( ~terminate ) then
     if this.tolsimplexizemethod == "enabled" then
       ssize = optimsimplex_size ( simplex , "sigmaplus" );
       this.optbase = optimbase_stoplog  ( this.optbase,sprintf("  > simplex size=%e < %e + %e * %e",...
         ssize, this.tolsimplexizeabsolute , this.tolsimplexizerelative , this.simplexsize0 ));
       if ssize < this.tolsimplexizeabsolute + this.tolsimplexizerelative * this.simplexsize0 then
-        terminate = 1;
+        terminate = %f;
         status = "tolsize";
       end
     end
@@ -490,7 +492,7 @@ function [this , terminate , status ] = neldermead_termination (this , ...
   //
   // Criteria #7 : simplex absolute size + difference in function values (Matlab-like)
   //
-  if (terminate == 0) then
+  if ( ~terminate ) then
     if this.tolssizedeltafvmethod == "enabled" then
       ssize = optimsimplex_size ( simplex , "sigmaplus" );
       this.optbase = optimbase_stoplog  ( this.optbase,sprintf("  > simplex size=%e < %e",...
@@ -499,7 +501,7 @@ function [this , terminate , status ] = neldermead_termination (this , ...
       this.optbase = optimbase_stoplog  ( this.optbase,sprintf("  > abs(fv(n+1) - fv(1))=%e < toldeltafv=%e",...
         shiftfv, this.toldeltafv));
       if ssize < this.tolsimplexizeabsolute & shiftfv < this.toldeltafv then
-        terminate = 1;
+        terminate = %f;
         status = "tolsizedeltafv";
       end
     end
@@ -508,7 +510,7 @@ function [this , terminate , status ] = neldermead_termination (this , ...
   // Criteria #8 : Kelley stagnation, based on
   // a sufficient decrease condition
   //
-  if ( terminate == 0 ) then
+  if ( ~terminate ) then
     if ( this.kelleystagnationflag==1 ) then
       [ sg , this ] = optimsimplex_gradientfv ( simplex , neldermead_costf , "forward" , this );
       nsg = sg.' * sg;
@@ -517,11 +519,47 @@ function [this , terminate , status ] = neldermead_termination (this , ...
       this.optbase = optimbase_stoplog ( this.optbase , ...
         sprintf ( "Test Stagnation : newfvmean=%e >= oldfvmean=%e - %e * %e" , newfvmean, oldfvmean , this.kelleyalpha , nsg ) );
       if ( newfvmean >= oldfvmean - this.kelleyalpha * nsg ) then
-        terminate = 1;
+        terminate = %f;
         status = "kelleystagnation";
       end
     end
   end
+  //
+  // Criteria #9 : Box termination criteria
+  // The number of consecutive time that an absolute tolerance on
+  // function value is met.
+  // From Algorithm 454, the tolerance is the difference between the
+  // max and the min function values in the simplex.
+  //
+  if {$terminate == 0} then {
+    if ( this.boxtermination ) then
+      shiftfv = abs(optimsimplex_deltafvmax( simplex ))
+      if ( shiftfv < this.boxtolf ) then
+        incr boxkount
+        set boxnbmatch [$self cget -boxnbmatch]
+        if {$boxkount == $boxnbmatch} then {
+          set terminate 1
+          set status "boxtolf"
+        }
+      else
+        set boxkount 0
+      end
+    end
+  end
+  //
+  // Criteria #10 : variance of function values
+  //
+  if {$terminate == 0} then {
+    if ( this.tolvarianceflag ) then
+      set tav [$self cget -tolabsolutevariance]
+      set var [optimsimplex::variance $simplex]
+      $optimpb stoplog "Test tolvariance : $var < $tav"
+      if {$var < $tav} then {
+        set terminate 1
+        set status "tolvariance"
+      }
+    end
+  end
 endfunction
   
 
@@ -671,7 +709,7 @@ function this = neldermead_startup (this)
   // user while he is tuning his object.
   checkfun = this.checkcostfunction;
   if checkfun == 1 then
-    optimbase_checkcostfun ( this.optbase );
+    this.optbase = optimbase_checkcostfun ( this.optbase );
   end
   // 1. If the problem has bounds, check that they are consistent
   [ this.optbase , hasbounds ] = optimbase_hasbounds ( this.optbase );
@@ -747,7 +785,7 @@ function this = neldermead_startup (this)
 endfunction
 //
 // neldermead_scaletox0 --
-//   Scale the simplex into the bounds and the 
+//   Scale the simplex into the 
 //   nonlinear inequality constraints, if any.
 //   Scale toward x0, which is feasible.
 // Arguments
@@ -778,7 +816,7 @@ function [ this , simplex0 ] = neldermead_scaletox0 ( this , simplex0 )
 endfunction
 //
 // neldermead_scaletocenter --
-//   Scale the simplex into the bounds and the 
+//   Scale the simplex into the 
 //   nonlinear inequality constraints, if any.
 //   Scale to the centroid of the points
 //   which satisfy the constraints.
@@ -845,86 +883,129 @@ function this = neldermead_kelleystag ( this )
     end
 endfunction
   //
-  // _scaleinconstraints --
+  // _scaleinboundsandcons --
   //   Given a point to scale and a reference point which satisfies the constraints, 
-  //   scale the point towards the reference point until it satisfies all the constraints.
-  //   Returns a list of key values with the following
-  //   keys : -status 0/1 -x x, where status
-  //   is 0 if the procedure has failed after -boxnbnlloops
+  //   scale the point towards the reference point until it satisfies all the constraints,
+  //   including boun constraints.
+  //   Returns isscaled = %T if the procedure has succeded before -boxnbnlloops
+  //   Returns isscaled = %F if the procedure has failed after -boxnbnlloops
   //   iterations.
   // Arguments
   //   x : the point to scale
   //   xref : the reference point
-  //   status : %T or %F
+  //   isscaled : %T or %F
   //   p : scaled point
   //
-function [ this , status , p ] = _scaleinconstraints ( this , x , xref )
+function [ this , isscaled , p ] = _scaleinboundsandcons ( this , x , xref )
   p = x
+  [ this.optbase , hasbounds ] = optimbase_hasbounds ( this.optbase );
+  nbnlc = optimbase_cget ( this.optbase , "-nbineqconst" )
   //
-  // Project the point into the bounds
+  // 1. No bounds, no nonlinear inequality constraints
+  // => no problem
+  //
+  if ( ( hasbounds == %f ) & ( nbnlc == 0 ) ) then
+    isscaled = %T
+    return;
+  end
+  isscaled = %F
+  //
+  // 2. Scale into bounds
   //
-  [ this.optbase , hasbounds ] = optimbase_hasbounds ( this.optbase );
   if ( hasbounds ) then
     [ this.optbase , p ] = optimbase_proj2bnds ( this.optbase ,  p );
     this = neldermead_log (this,sprintf(" > After projection into bounds p = [%s]" , ...
-      strcat(string(p)," ")));
+      _strvec(p)));
+  end
+  //
+  // 2. Scale into nonlinear constraints
+  // Try the current point and see if the constraints are satisfied.
+  // If not, move the point "halfway" to the centroid,
+  // which should satisfy the constraints, if
+  // the constraints are convex.
+  // Perform this loop until the constraints are satisfied.
+  // If all loops have been performed without success, the scaling
+  // has failed.
+  //
+  alpha = 1.0
+  p0 = p
+  while ( alpha > this.guinalphamin )
+      [ this.optbase , feasible ] = optimbase_isinnonlincons ( this.optbase , p );
+      if ( feasible ) then
+        isscaled = %T;
+        break;
+      end
+      alpha = alpha / 2.0
+      this = neldermead_log (this,sprintf("Scaling inequality constraint with alpha = %e", ...
+        alpha ));
+      p = ( 1.0 - alpha ) * xref + alpha * p0;
   end
+  this = neldermead_log (this,sprintf(" > After scaling into inequality constraints p = [%s]" , ...
+    _strvec(p) ) );
+  if ( ~isscaled ) then
+    this = neldermead_log (this,sprintf(" > Impossible to scale into constraints." ));
+  end
+endfunction
+  //
+  // _scaleinconstraints --
+  //   Given a point to scale and a reference point which satisfies the constraints, 
+  //   scale the point towards the reference point until it satisfies all the constraints.
+  //   Returns isscaled = %T if the procedure has succeded before -boxnbnlloops
+  //   Returns isscaled = %F if the procedure has failed after -boxnbnlloops
+  //   iterations.
+  // Arguments
+  //   x : the point to scale
+  //   xref : the reference point
+  //   isscaled : %T or %F
+  //   p : scaled point
+  //
+function [ this , isscaled , p ] = _scaleinconstraints ( this , x , xref )
+  p = x
   //
   // Adjust point to satisfy nonlinear inequality constraints
   //
   nbnlc = optimbase_cget ( this.optbase , "-nbineqconst" )
   if ( nbnlc == 0 ) then
-    scaled = %T
-  else
-    scaled = %F
-    //
-    // Try the current point and see if the constraints are satisfied.
-    // If not, move the point "halfway" to the centroid,
-    // which should satisfy the constraints, if
-    // the constraints are convex.
-    // Perform this loop until the constraints are satisfied.
-    // If all loops have been performed without success, the scaling
-    // has failed.
-    //
-    for i = 1 : this.nbineqloops
-      [ this , constlist ] = neldermead_function ( this , p , index=2 );
-      feasible = %T
-      for ic = 1 : this.optbase.nbineqconst;
-        const = constlist ( ic );
-        if ( const < 0.0 ) then
-          this = neldermead_log (this,sprintf("Constraint #%d/%d is not satisfied", ...
-            ic , this.optbase.nbineqconst ));
-          feasible = %F;
-          break;
-        end
-      end
+    isscaled = %T
+    return;
+  end
+  isscaled = %F
+  //
+  // Try the current point and see if the constraints are satisfied.
+  // If not, move the point "halfway" to the centroid,
+  // which should satisfy the constraints, if
+  // the constraints are convex.
+  // Perform this loop until the constraints are satisfied.
+  // If all loops have been performed without success, the scaling
+  // has failed.
+  //
+  alpha = 1.0
+  p0 = p
+  while ( alpha > this.guinalphamin )
+      [ this.optbase , feasible ] = optimbase_isinnonlincons ( this.optbase , p );
       if ( feasible ) then
-        scaled = %T;
+        isscaled = %T;
         break;
-      else
-        this = neldermead_log (this,sprintf("Scaling inequality constraint at loop #%d/%d", ...
-          i , this.nbineqloops));
-        p = ( xref + p ) * this.ineqscaling;
       end
-    end
-    this = neldermead_log (this,sprintf(" > After scaling into inequality constraints p = [%s]" , ...
-      strcat(string(p)," ") ) );
+      alpha = alpha / 2.0
+      this = neldermead_log (this,sprintf("Scaling inequality constraint with alpha = %e", ...
+        alpha));
+      p = ( 1.0 - alpha ) * xref + alpha * p0;
   end
-  if ( scaled ) then
-    status = %T
-  else
-    status = %F
+  this = neldermead_log (this,sprintf(" > After scaling into inequality constraints p = [%s]" , ...
+    _strvec(p) ) );
+  if ( ~isscaled ) then
     this = neldermead_log (this,sprintf(" > Impossible to scale into constraints after %d loops" , ...
       this.optbase.nbineqconst ));
   end
 endfunction
 //
-// neldermead_box --
+// neldermead_constraints --
 //   The Nelder-Mead algorithm, with variable-size simplex
 //   and modifications by Box for bounds and
 //   inequality constraints.
 //
-function this = neldermead_box ( this )
+function this = neldermead_constraints ( this )
   //
   // Order the vertices for the first time
   //
@@ -980,10 +1061,10 @@ function this = neldermead_box ( this )
     this = neldermead_log (this,sprintf("================================================================="));
     this = neldermead_log (this,sprintf("Iteration #%d (total = %d)",iter,totaliter));
     this = neldermead_log (this,sprintf("Function Eval #%d",funevals));
-    this = neldermead_log (this,sprintf("Xopt : [%s]",strcat(string(xlow)," ")));
+    this = neldermead_log (this,sprintf("Xopt : [%s]",_strvec(xlow)));
     this = neldermead_log (this,sprintf("Fopt : %e",flow));
     this = neldermead_log (this,sprintf("DeltaFv : %e",deltafv));
-    this = neldermead_log (this,sprintf("Center : [%s]",strcat(string(currentcenter)," ")));
+    this = neldermead_log (this,sprintf("Center : [%s]",_strvec(currentcenter)));
     this = neldermead_log (this,sprintf("Size : %e",ssize));
     str = optimsimplex_tostring ( simplex )
     for i = 1:nbve
@@ -997,7 +1078,7 @@ function this = neldermead_box ( this )
     if ( iter > 1 ) then
       [this , terminate , status] = neldermead_termination (this , ...
         fvinitial , oldfvmean , newfvmean , previouscenter , currentcenter , simplex );
-      if (terminate==1) then
+      if ( terminate ) then
         this = neldermead_log (this,sprintf("Terminate with status : %s",status));
         break
       end
@@ -1007,21 +1088,22 @@ function this = neldermead_box ( this )
     //
     this = neldermead_log (this,sprintf("Reflect"));
     xbar = optimsimplex_xbar ( simplex );
-    this = neldermead_log (this,sprintf("xbar=[%s]" , strcat(string(xbar)," ")));
+    this = neldermead_log (this,sprintf("xbar=[%s]" , _strvec(xbar)));
     //
     // Reflect the worst point with respect to center
     //
     xr = neldermead_interpolate ( xbar , xhigh , this.rho );
+    this = neldermead_log (this,sprintf("xr=[%s]" , _strvec(xr)));
     // Adjust point to satisfy bounds and nonlinear inequality constraints
     if ( hasbounds | nbnlc > 0 ) then
-      [ this , status , xr ] = _scaleinconstraints ( this , xr , xbar )
+      [ this , status , xr ] = _scaleinboundsandcons ( this , xr , xbar );
       if ( ~status ) then
         status = "impossibleconstr"
         break
       end
     end
     [ this , fr ] = neldermead_function ( this , xr );
-    this = neldermead_log (this,sprintf("xr=[%s], f(xr)=%f", strcat(string(xr)," ") , fr));
+    this = neldermead_log (this,sprintf("xr=[%s], f(xr)=%f", _strvec(xr) , fr));
     if ( fr >= flow & fr < fn ) then
       this = neldermead_log (this,sprintf("  > Perform reflection"));
       simplex = optimsimplex_setve ( simplex , ihigh , fr , xr )
@@ -1032,7 +1114,7 @@ function this = neldermead_box ( this )
       xe = neldermead_interpolate ( xbar , xhigh , this.rho*this.chi );
       // Adjust point to satisfy bounds and nonlinear inequality constraints
       if ( hasbounds | nbnlc > 0 ) then
-        [ this , status , xe ] = _scaleinconstraints ( this , xe , xbar )
+        [ this , status , xe ] = _scaleinboundsandcons ( this , xe , xbar );
         if ( ~status ) then
           status = "impossibleconstr"
           break
@@ -1055,7 +1137,7 @@ function this = neldermead_box ( this )
       xc = neldermead_interpolate ( xbar , xhigh , this.rho*this.gamma );
       // Adjust point to satisfy bounds and nonlinear inequality constraints
       if ( hasbounds | nbnlc > 0 ) then
-        [ this , status , xc ] = _scaleinconstraints ( this , xc , xbar )
+        [ this , status , xc ] = _scaleinboundsandcons ( this , xc , xbar );
         if ( ~status ) then
           status = "impossibleconstr"
           break
@@ -1080,7 +1162,7 @@ function this = neldermead_box ( this )
       xc = neldermead_interpolate ( xbar , xhigh , -this.gamma );
       // Adjust point to satisfy bounds and nonlinear inequality constraints
       if ( hasbounds | nbnlc > 0 ) then
-        [ this , status , xc ] = _scaleinconstraints ( this , xc , xbar )
+        [ this , status , xc ] = _scaleinboundsandcons ( this , xc , xbar );
         if ( ~status ) then
           status = "impossibleconstr"
           break
@@ -1110,3 +1192,254 @@ function this = neldermead_box ( this )
   this.optbase = optimbase_set ( this.optbase , "-status" , status );
   this.simplexopt = simplex;
 endfunction
+//
+// neldermead_box --
+//   The Nelder-Mead algorithm, with variable-size simplex
+//   and modifications by Box for bounds and
+//   inequality constraints.
+//
+function this = neldermead_box ( this )
+  //
+  // Order the vertices for the first time
+  //
+  simplex = this.simplex0;
+  n = optimbase_cget ( this.optbase , "-numberofvariables" );
+  fvinitial = optimbase_get ( this.optbase , "-fx0" );
+  // Sort function values and x points by increasing function value order
+  this = neldermead_log (this,"Step #1 : order");
+  simplex = optimsimplex_sort ( simplex );
+  currentcenter = optimsimplex_center ( simplex );
+  currentxopt = optimbase_cget ( this.optbase , "-x0" );
+  newfvmean = optimsimplex_fvmean ( simplex );
+  nbve = optimsimplex_getnbve ( simplex );
+  ihigh = nbve;
+  inext = ihigh - 1
+  ilow = 1
+  [ this.optbase , hasbounds ] = optimbase_hasbounds ( this.optbase );
+  nbnlc = optimbase_cget ( this.optbase , "-nbineqconst" )
+  rho = this.boxreflect;
+  //
+  // Initialize
+  //
+  terminate = 0;
+  iter = 0;
+  step = "init";
+  //
+  // Nelder-Mead Loop
+  //
+  while ( terminate == 0 )
+    this.optbase = optimbase_incriter ( this.optbase );
+    iter = iter + 1;
+    xlow = optimsimplex_getx ( simplex , ilow )
+    flow = optimsimplex_getfv ( simplex , ilow )
+    xhigh = optimsimplex_getx ( simplex , ihigh )
+    fhigh = optimsimplex_getfv ( simplex , ihigh )
+    xn = optimsimplex_getx ( simplex , inext )
+    fn = optimsimplex_getfv ( simplex , inext )
+    //
+    // Store history
+    //
+    xcoords = optimsimplex_getallx ( simplex )
+    this = neldermead_storehistory ( this , n , flow , xlow , xcoords );
+    deltafv = abs(optimsimplex_deltafvmax ( simplex ));
+    currentfopt = flow;
+    previousxopt = currentxopt;
+    currentxopt = xlow;
+    previouscenter = currentcenter;
+    currentcenter = optimsimplex_center ( simplex );
+    oldfvmean = newfvmean;
+    newfvmean = optimsimplex_fvmean ( simplex );
+    totaliter = optimbase_get ( this.optbase , "-iterations" );
+    funevals = optimbase_get ( this.optbase , "-funevals" );
+    ssize = optimsimplex_size ( simplex )
+    this = neldermead_log (this,sprintf("================================================================="));
+    this = neldermead_log (this,sprintf("Iteration #%d (total = %d)",iter,totaliter));
+    this = neldermead_log (this,sprintf("Function Eval #%d",funevals));
+    this = neldermead_log (this,sprintf("Xopt : [%s]",_strvec(xlow)));
+    this = neldermead_log (this,sprintf("Fopt : %e",flow));
+    this = neldermead_log (this,sprintf("DeltaFv : %e",deltafv));
+    this = neldermead_log (this,sprintf("Center : [%s]",_strvec(currentcenter)));
+    this = neldermead_log (this,sprintf("Size : %e",ssize));
+    str = optimsimplex_tostring ( simplex )
+    for i = 1:nbve
+      this = neldermead_log (this,str(i));
+    end
+    neldermead_outputcmd ( this, "iter" , simplex , step )
+
+    //
+    // Update termination flag
+    //
+    if ( iter > 1 ) then
+      [this , terminate , status] = neldermead_termination (this , ...
+        fvinitial , oldfvmean , newfvmean , previouscenter , currentcenter , simplex );
+      if ( terminate ) then
+        this = neldermead_log (this,sprintf("Terminate with status : %s",status));
+        break
+      end
+    end
+    //
+    // Compute xbar, center of better vertices
+    //
+    this = neldermead_log (this,sprintf("Reflect"));
+    xbar = optimsimplex_xbar ( simplex );
+    this = neldermead_log (this,sprintf("xbar=[%s]" , _strvec(xbar)));
+    //
+    // Search a point improving cost function
+    // and satisfying constraints.
+    //
+    [ this , scaled , xr , fr ] = _boxlinesearch ( this , n , xbar , xhigh , fhigh , rho );
+    if ( scaled == %f ) then
+      status = "impossibleimprovement"
+      break
+    end
+    this = neldermead_log (this,sprintf("xr=[%s], f(xr)=%f", strcat(string(xr)," ") , fr));
+    this = neldermead_log (this,sprintf("  > Perform Reflection"));
+    simplex = optimsimplex_setve ( simplex , ihigh , fr , xr )
+    step = "boxreflection";
+    //
+    // Sort simplex
+    //
+    this = neldermead_log (this,sprintf("Sort"));
+    simplex  = optimsimplex_sort ( simplex );
+  end
+  this.optbase = optimbase_set ( this.optbase , "-xopt" , xlow.' );
+  this.optbase = optimbase_set ( this.optbase , "-fopt" , flow );
+  this.optbase = optimbase_set ( this.optbase , "-status" , status );
+  this.simplexopt = simplex;
+endfunction
+
+  //
+  // _strvec --
+  //  Returns a string for the given vector.
+  //
+  function str = _strvec ( x )
+    str = strcat(string(x)," ")
+  endfunction
+  //
+  // _boxlinesearch --
+  //   For Box's method, perform a line search
+  //   from xbar, on the line (xhigh,xbar) and returns:
+  //   status : %t if the search is successful
+  //   xr : the reflected point
+  //   fr : the function value
+  //   The reflected point satisfies the following
+  //   constraints :
+  //   * fr < fhigh
+  //   * xr satisfies the bounds constraints
+  //   * xr satisfies the nonlinear positive inequality constraints
+  //   * xr satisfies the linear positive inequality constraints
+  //   The method is based on projection and
+  //   scaling toward the centroid.
+  //
+  // Arguments
+  //   n : number of variables
+  //   xbar : the centroid
+  //   xhigh : the high point
+  //   fhigh : function value at xhigh
+  //   rho : reflection factor
+  //
+  function [ this , status , xr , fr ] = _boxlinesearch ( this , n , xbar , xhigh , fhigh , rho )
+    this = neldermead_log (this,"_boxlinesearch");
+    this = neldermead_log (this, sprintf ("> xhigh=[%s], fhigh=%e",_strvec(xhigh),fhigh));
+    this = neldermead_log (this, sprintf ( "> xbar=[%s]" , _strvec(xbar) ) );
+    xr = neldermead_interpolate ( xbar , xhigh , rho );
+    this = neldermead_log (this, sprintf ( "> xr = [%s]" , _strvec ( xr ) ) );
+    status = %f
+    alphamin = this.guinalphamin
+    //
+    // Scale from xr toward xbar until fr < fhigh and update xr
+    //
+    xr0 = xr
+    alpha = 1.0
+    while ( alpha > alphamin )
+      [ this , fr ] = neldermead_function ( this , xr );
+      if ( fr < fhigh ) then
+        this = neldermead_log (this, sprintf ( "fr = %e improves %e : no need for scaling for f" , fr , fhigh ) );
+        status = %t;
+        break
+      end
+      alpha = alpha / 2.0;
+      this = neldermead_log (this, sprintf ( "Scaling for f with alpha=%e" , alpha ) );
+      xr = ( 1.0 - alpha ) * xbar + alpha * xr0;
+      this = neldermead_log (this, sprintf ( "> xr = %s" , _strvec ( xr ) ) );
+    end
+    // If the scaling for function improvement has failed,
+    // we return.
+    if ( ~status ) then
+      return;
+    end
+    // scaledc is set to %t if xr is updated during scaling into constraints 
+    // That implies that the function value is to update.
+    scaledc = %f
+    //
+    // Project xr into bounds, with an additionnal alpha inside the bounds.
+    // This algo is always succesful.
+    // Note:
+    //   If the alpha coefficient was not used, the
+    //   projectinbounds method could be used directly.
+    //
+    [ this.optbase , hasbounds ] = optimbase_hasbounds ( this.optbase );
+    if ( hasbounds ) then
+      boxboundsalpha = this.boxboundsalpha;
+      boundsmax = optimbase_cget ( this.optbase , "-boundsmax" );
+      boundsmin = optimbase_cget ( this.optbase , "-boundsmin" );
+      for ix = 1:n
+        xmin = boundsmin ( ix );
+        xmax = boundsmax ( ix );
+        xrix = xr ( ix );
+        if ( xrix > xmax ) then
+          this = neldermead_log (this, sprintf ( "Projecting index #%d = %e on max bound %e - %e" , ix , xrix , xmax , boxboundsalpha ) );
+          xr ( ix ) = xmax - boxboundsalpha;
+          if ( ~scaledc ) then
+            scaledc = %t
+          end
+        elseif ( xrix < xmin ) then
+          this = neldermead_log (this, sprintf ( "Projecting index #%e = %e on min bound %e - %e" , ix , xrix , xmin , boxboundsalpha ) );
+          xr ( ix ) = xmin + boxboundsalpha;
+          if ( ~scaledc ) then
+            scaledc = %t
+          end
+        end
+      end
+      this = neldermead_log (this, sprintf ( " > After projection into bounds xr = [%s]" , _strvec(xr)));
+    end
+    //
+    // Scale from xr to xbar into nonlinear inequality constraints
+    // and update xr. 
+    // Set status to 0 if the process fails.
+    //
+    nbnlc = optimbase_cget ( this.optbase , "-nbineqconst" );
+    if ( nbnlc == 0 ) then
+      status = %t
+    else
+      status = %f;
+      alpha = 1.0;
+      xr0 = xr;
+      while ( alpha > alphamin )
+        [ this.optbase , feasible ] = optimbase_isinnonlincons ( this.optbase , xr );
+        if ( feasible ) then
+          status = %t;
+          break
+        end
+        alpha = alpha / 2.0;
+        this = neldermead_log (this, sprintf ( "Scaling for nonlinear/linear inequality constraints with alpha=%e toward [%s]" , ...
+          alpha , _strvec(p0) ));
+        xr = ( 1.0 - alpha ) * xbar + alpha * xr0;
+        this = neldermead_log (this, sprintf ( "> xr = [%s]" , _strvec(xr) ));
+        if ( ~scaledc ) then
+          scaledc = %t;
+        end
+      end
+    end
+    // If scaling failed, returns immediately 
+    // (there is no need to update the function value).
+    if ( ~status ) then
+      return
+    end
+    if ( scaledc ) then
+      // Re-compute the function value at scaled point
+      [ this , fr ] = neldermead_function ( this , xr );
+    end
+    
+  endfunction
+
index 357666b..06c28f5 100644 (file)
@@ -18,7 +18,7 @@ function this = neldermead_updatesimp ( this )
   xopt = optimbase_get ( this.optbase , "-xopt" );
   select this.restartsimplexmethod
   case "oriented" then
-    [ simplex0 , this ] = optimsimplex_new ( "oriented" , this.simplexopt , neldermead_costf , this )
+    [ simplex0 , this ] = optimsimplex_new ( "oriented" , this.simplexopt , neldermead_costf , this );
   case "axes" then
     simplex0 = optimsimplex_destroy ( simplex0 )
     [ simplex0 , this ] = optimsimplex_new ( "axes" , ...
@@ -65,12 +65,11 @@ function this = neldermead_updatesimp ( this )
     nbve = optimsimplex_getnbve ( simplex0 )
     for ive = 1 : nbve
       // x is a row vector
-      x = optimsimplex_getx ( simplex0 , ive )
-      this = neldermead_log (this,sprintf("Scaling vertex #%d/%d at ["+...
-        strcat(string(x)," ")+"]... " , ...
-        ive , nbve ));
+      x = optimsimplex_getx ( simplex0 , ive );
+      this = neldermead_log (this,sprintf("Scaling vertex #%d/%d at [%s]... " , ...
+        ive , nbve , _strvec(x)));
       // Transpose x because xopt is a column vector : xp is now a column vector
-      [ this , status , xp ] = _scaleinconstraints ( this , x.' , xopt )
+      [ this , status , xp ] = _scaleinconstraints ( this , x.' , xopt );
       if ( ~status ) then
         errmsg = msprintf(gettext("Impossible to scale the vertex #%d/%d at [%s] into inequality constraints"), "neldermead_updatesimp", ...
           ive , nbve , strcat(string(x)," "));
index 29843ca..b223013 100644 (file)
@@ -57,6 +57,8 @@ function value = optimbase_cget (this,key)
     value = this.boundsmax;
   case "-nbineqconst" then
     value = this.nbineqconst;
+  case "-logfile" then
+    value = this.logfile;
   else
     errmsg = msprintf(gettext("%s: Unknown key %s") , "optimbase_cget" , key)
     error(errmsg)
index f915383..1e19516 100644 (file)
 // Arguments
 //   flag : %T or %F
 //
-function [ opt , isok ] = optimbase_checkx0 ( this )
-    optimbase_log ( this , sprintf ( "Checking initial guess..." ) )
+function [ this , isok ] = optimbase_checkx0 ( this )
+    this = optimbase_log ( this , sprintf ( "Checking initial guess..." ) )
     [ this , isfeasible ] = optimbase_isfeasible ( this , this.x0 )
     isok = ( isfeasible == 1 )
     if ( isok ) then
-      optimbase_log ( this , sprintf ( "... initial guess is feasible." ) )
+      this = optimbase_log ( this , sprintf ( "... initial guess is feasible." ) )
     else
-      optimbase_log ( this , sprintf ( "... initial guess is not feasible." ) )
+      this = optimbase_log ( this , sprintf ( "... initial guess is not feasible." ) )
     end
 endfunction
 
index b95cabf..eabe473 100644 (file)
@@ -81,6 +81,8 @@ function this = optimbase_configure (this,key,value)
     this.boundsmax = value;
   case "-nbineqconst" then
     this.nbineqconst = value;
+  case "-logfile" then
+    this.logfile = value;
   else
     errmsg = msprintf(gettext("%s: Unknown key %s"),"optimbase_configure",key)
     error(errmsg)
index 1626724..a17669f 100644 (file)
@@ -14,5 +14,8 @@
 function this = optimbase_destroy (this)
   this.historyfopt = [];
   this.historyxopt = [];
+  if ( this.logstartup ) then
+    this = optimbase_logshutdown ( this );
+  end
 endfunction
 
index 18b732c..a9001da 100644 (file)
@@ -52,23 +52,23 @@ function [ this , result ] = optimbase_function ( this , x , index )
       if ( this.nbineqconst == 0 ) then
         // There is one additionnal argument, but no nonlinear constraints,
         // therefore, there is no need for a index value.
-        result = this.fun ( x , this.costfargument );
+        [ result , this.costfargument ] = this.fun ( x , this.costfargument );
       else
         // Set the index, so that, if an additionnal cost function argument is provided,
         // it can be appended at the end.
         index = 1;
-        result = this.fun ( x , index , this.costfargument );
+        [ result , this.costfargument ] = this.fun ( x , index , this.costfargument );
       end
     else
       // There is one additionnal argument, and the caller provided the value of index.
-      result = this.fun ( x , index , this.costfargument );
+      [ result , this.costfargument ] = this.fun ( x , index , this.costfargument );
     end
   end
   this.funevals = this.funevals + 1;
   if this.verbose == 1 then
     msg = sprintf ( "Function Evaluation #%d is [%s] at [%s]" , ...
       this.funevals , strcat(string(result)," ") , strcat(string(x)," ") )
-    optimbase_log ( this , msg )
+    this = optimbase_log ( this , msg )
   end
 endfunction
 
index 1be98a6..37d8543 100644 (file)
@@ -30,13 +30,13 @@ function [ this , isfeasible ] = optimbase_isfeasible ( this , x )
           xix = x ( ix )
           if ( xix < xmin ) then
             isfeasible = 0
-            optimbase_log ( this , sprintf ( "Component #%d/%d of x is lower than min bound %e", ...
+            this = optimbase_log ( this , sprintf ( "Component #%d/%d of x is lower than min bound %e", ...
               ix , this.numberofvariables , xmin ) )
             break
           end
           if (xix > xmax) then
             isfeasible = 0
-            optimbase_log ( this , sprintf ( "Component #%d/%d of x is greater than max bound %e", ...
+            this = optimbase_log ( this , sprintf ( "Component #%d/%d of x is greater than max bound %e", ...
               ix , this.numberofvariables , xmax ) )
             break
           end
@@ -52,7 +52,7 @@ function [ this , isfeasible ] = optimbase_isfeasible ( this , x )
         index = 0
         for ic = 1 : this.nbineqconst
           if ( const ( ic ) < 0.0 ) then
-            optimbase_log ( this , sprintf ( "Inequality constraint #%d/%d is not satisfied for x", ...
+            this = optimbase_log ( this , sprintf ( "Inequality constraint #%d/%d is not satisfied for x", ...
               ic , this.nbineqconst ) )
             isfeasible = -1
             break
index 0769708..3d605fc 100644 (file)
@@ -25,13 +25,13 @@ function [ this , isfeasible ] = optimbase_isinbounds ( this , x )
           xix = x ( ix )
           if ( xix < xmin ) then
             isfeasible = %f
-            optimbase_log ( this , sprintf ( "Component #%d/%d of x is lower than min bound %e", ...
+            this = optimbase_log ( this , sprintf ( "Component #%d/%d of x is lower than min bound %e", ...
               ix , this.numberofvariables , xmin ) )
             break
           end
           if (xix > xmax) then
             isfeasible = %f
-            optimbase_log ( this , sprintf ( "Component #%d/%d of x is greater than max bound %e", ...
+            this = optimbase_log ( this , sprintf ( "Component #%d/%d of x is greater than max bound %e", ...
               ix , this.numberofvariables , xmax ) )
             break
           end
index f0bd793..442dac2 100644 (file)
@@ -22,7 +22,7 @@ function [ this , isfeasible ] = optimbase_isinnonlincons ( this , x )
         index = 0
         for ic = 1 : this.nbineqconst
           if ( const ( ic ) < 0.0 ) then
-            optimbase_log ( this , sprintf ( "Inequality constraint #%d/%d is not satisfied for x", ...
+            this = optimbase_log ( this , sprintf ( "Inequality constraint #%d/%d is not satisfied for x", ...
               ic , this.nbineqconst ) )
             isfeasible = %f
             break
index 728eca7..51e2ed7 100644 (file)
 //
 function this = optimbase_log (this,msg)
   if this.verbose == 1 then
-     mprintf("%s\n",msg);
+    if ( ~this.logstartup ) then
+      this = optimbase_logstartup ( this )
+    end
+    if ( this.logfile <> "" ) then
+      mfprintf ( this.logfilehandle , "%s\n" , msg );
+    else
+      mprintf("%s\n",msg);
+    end
   end
 endfunction
 
diff --git a/scilab/modules/optimization/macros/optimbase/optimbase_logshutdown.sci b/scilab/modules/optimization/macros/optimbase/optimbase_logshutdown.sci
new file mode 100644 (file)
index 0000000..0fd88fc
--- /dev/null
@@ -0,0 +1,26 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008-2009 - INRIA - Michael Baudin
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+
+//
+// optimbase_logshutdown --
+//   Shut down the logging.
+//   If the logging is already started up, generates an error.
+//   If the loggin is started up, if there is a log file, close it.
+//
+function this = optimbase_logshutdown ( this )
+  if ~this.logstartup then
+    error ( gettext ( "%s: Logging not started." , "optimbase_logstartup" ) )
+  else
+    this.logstartup = %f;
+    if ( this.logfile <> "" ) then
+      mclose( this.logfilehandle );
+    end
+  end
+endfunction
+
diff --git a/scilab/modules/optimization/macros/optimbase/optimbase_logstartup.sci b/scilab/modules/optimization/macros/optimbase/optimbase_logstartup.sci
new file mode 100644 (file)
index 0000000..523b88d
--- /dev/null
@@ -0,0 +1,39 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008-2009 - INRIA - Michael Baudin
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+
+//
+// optimbase_logstartup --
+//   Automatically startup logging.
+//   If the logging is already started up, generates an error.
+//   If the logging is not started up, if there is a log file configured,
+//   open that log file to append messages.
+//
+function this = optimbase_logstartup ( this )
+  if this.logstartup then
+    error ( gettext ( "%s: Logging already started." , "optimbase_logstartup" ) )
+  else
+    this.logstartup = %t;
+    if ( this.logfile <> "" ) then
+      if ( this.logfilehandle <> 0 ) then
+        error ( gettext ( "%s: Logging not started up, log file not empty but log file handle non zero." , "optimbase_logstartup" ) )
+      end
+      this.logfilehandle = mopen( this.logfile , "a" );
+      c = clock();
+      year = c(1);
+      month = c(2);
+      day = c(3);
+      hour = c(4);
+      minute = c(5);
+      seconds = c(6);
+      mfprintf ( this.logfilehandle , "Optimbase ready for logging at %d-%d-%d %d:%d:%d\n" , ...
+        year , month , day , hour , minute , seconds );
+    end
+  end
+endfunction
+
index 1743b9f..0f00a20 100644 (file)
 //   Creates a new Optimization object.
 //
 function newobj = optimbase_new ()
-  newobj = tlist(["T_OPTIMIZATION",...
-    "verbose",...
-    "x0","fx0",...
-    "xopt","fopt",...
-    "tolfunabsolute","tolfunrelative","tolfunmethod",...
-    "tolxabsolute","tolxrelative","tolxmethod" ,...
-    "funevals","maxfunevals",...
-    "iterations","maxiter",...
-    "fun","status", ...
-    "historyfopt","historyxopt",...
-    "verbosetermination",...
-    "outputcommand",...
-    "outputcommandarg","numberofvariables",...
-    "storehistory","costfargument",...
-    "boundsmin","boundsmax",...
-    "nbineqconst"]);
+  newobj = tlist(["T_OPTIMIZATION"
+    "verbose"
+    "x0"
+    "fx0"
+    "xopt"
+    "fopt"
+    "tolfunabsolute"
+    "tolfunrelative"
+    "tolfunmethod"
+    "tolxabsolute"
+    "tolxrelative"
+    "tolxmethod"
+    "funevals"
+    "maxfunevals"
+    "iterations"
+    "maxiter"
+    "fun"
+    "status"
+    "historyfopt"
+    "historyxopt"
+    "verbosetermination"
+    "outputcommand"
+    "outputcommandarg"
+    "numberofvariables"
+    "storehistory"
+    "costfargument"
+    "boundsmin"
+    "boundsmax"
+    "nbineqconst"
+    "logfile"
+    "logfilehandle"
+    "logstartup"
+    ]);
   // The number of variables to optimize
   newobj.numberofvariables = 0
   // The verbose option, controlling the amount of messages
@@ -94,4 +111,11 @@ function newobj = optimbase_new ()
   newobj.boundsmin = [];
   // The number of nonlinear inequality constraints
   newobj.nbineqconst = 0;
+  // The name of the log file
+  newobj.logfile = "";
+  // The handle for the log file
+  newobj.logfilehandle = 0;
+  // Set to %t when the logging is started up
+  newobj.logstartup = %f;
 endfunction
+
index 3598b41..a799fc2 100644 (file)
@@ -29,10 +29,10 @@ function [ this , p ] = optimbase_proj2bnds ( this ,  x )
       xmax = this.boundsmax ( ix )
       pix = p ( ix )
       if (pix > xmax) then
-        optimbase_log ( this , sprintf ( "Projecting p(%d) = %e on max bound %e" , ix , pix , xmax ))
+        this = optimbase_log ( this , sprintf ( "Projecting p(%d) = %e on max bound %e" , ix , pix , xmax ))
         p ( ix ) = xmax
       elseif ( pix < xmin) then
-        optimbase_log ( this , sprintf ( "Projecting p(%d) = %e on min bound %e" , ix , pix , xmin ))
+        this = optimbase_log ( this , sprintf ( "Projecting p(%d) = %e on min bound %e" , ix , pix , xmin ))
         p ( ix ) = xmin
       end
     end
index 0e65300..bd16904 100644 (file)
@@ -15,7 +15,7 @@
 function this = optimbase_stoplog ( this , msg )
   if this.verbose == 1 then
   if this.verbosetermination == 1 then
-     mprintf("%s\n",msg);
+    this = optimbase_log ( this , msg )
   end
   end
 endfunction
index 4ea4b6e..dd97063 100644 (file)
@@ -10,8 +10,8 @@
 
 //
 // optimbase_terminate --
-//   Returns 1 if the algorithm terminates.
-//   Returns 0 if the algorithm must continue.
+//   Returns %t if the algorithm terminates.
+//   Returns %f if the algorithm must continue.
 // Arguments, input
 //   this : the current object
 //   previousfopt : the previous value of the objective function
 //
 function [this , terminate , status] = optimbase_terminate (this , ...
   previousfopt , currentfopt , previousxopt , currentxopt )
-  terminate = 0;
+  terminate = %f;
   status = "continue";
   this = optimbase_stoplog (this,sprintf("  > Termination ?"));
   //
   // Criteria #1 : maximum number of iterations
   //
-  if (terminate == 0) then 
+  if ( ~terminate ) then 
     this = optimbase_stoplog (this,sprintf("  > iterations=%d >= maxiter=%d",this.iterations, this.maxiter));
     if this.iterations >= this.maxiter then
-      terminate = 1;
+      terminate = %t;
       status = "maxiter";
     end
   end
   //
   // Criteria #2 : maximum number of call to function
   //
-  if (terminate == 0) then 
+  if ( ~terminate ) then 
     this = optimbase_stoplog (this,sprintf("  > funevals=%d >= maxfunevals=%d",this.funevals, this.maxfunevals));
     if this.funevals >= this.maxfunevals then
-      terminate = 1;
+      terminate = %t;
       status = "maxfuneval";
     end
   end
@@ -64,13 +64,13 @@ function [this , terminate , status] = optimbase_terminate (this , ...
   //   and the optimum function value is strictly negative (e.g. f(x*)=-10),
   //   that criteria fails miserably.
   //
-  if (terminate == 0) then
+  if ( ~terminate ) then
     select this.tolfunmethod
     case "enabled" then
       this = optimbase_stoplog (this,sprintf("  > abs(currentfopt)=%e < tolfunrelative * abs(previousfopt) + tolfunabsolute=%e",...
         abs(currentfopt), this.tolfunrelative * abs(previousfopt) + this.tolfunabsolute));
       if abs(currentfopt) < this.tolfunrelative * abs(previousfopt) + this.tolfunabsolute then
-        terminate = 1;
+        terminate = %t;
         status = "tolf";
       end
     case "disabled" then
@@ -83,7 +83,7 @@ function [this , terminate , status] = optimbase_terminate (this , ...
   //
   // Criteria #4 : tolerance on x
   //
-  if (terminate == 0) then
+  if ( ~terminate ) then
     //
     // Note
     // What means a relative error on x ?
@@ -97,7 +97,7 @@ function [this , terminate , status] = optimbase_terminate (this , ...
       this = optimbase_stoplog (this,sprintf("  > e(x)=%e < tolxrelative = %e * %e + %e",...
         norm(currentxopt - previousxopt), this.tolxrelative , norm(currentxopt) , this.tolxabsolute ));
       if norm(currentxopt - previousxopt) < this.tolxrelative * norm(currentxopt) + this.tolxabsolute then
-        terminate = 1;
+        terminate = %t;
         status = "tolx";
       end
    case "disabled" then
index e2b2080..f6f4aed 100644 (file)
@@ -90,7 +90,7 @@ function result = optimtestcase ( x , index )
   end
 endfunction
 //
-// Test with Box algorithm and default axes initial simplex
+// Test with NM-constrained algorithm and default axes initial simplex
 //
 nm = neldermead_new ();
 nm = neldermead_configure(nm,"-numberofvariables",4);
@@ -100,7 +100,7 @@ nm = neldermead_configure(nm,"-maxiter",200);
 nm = neldermead_configure(nm,"-maxfunevals",400);
 nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-3);
 nm = neldermead_configure(nm,"-simplex0method","axes");
-nm = neldermead_configure(nm,"-method","box");
+nm = neldermead_configure(nm,"-method","nmconstraints");
 nm = neldermead_configure(nm,"-nbineqconst",3);
 //nm = neldermead_configure(nm,"-verbose",1);
 nm = neldermead_configure(nm,"-verbosetermination",1);
@@ -116,7 +116,7 @@ status = neldermead_get(nm,"-status");
 assert_equal ( status , "tolsize" );
 nm = neldermead_destroy(nm);
 //
-// Test with Box algorithm and restart
+// Test with NM-constrained algorithm and restart
 //
 nm = neldermead_new ();
 nm = neldermead_configure(nm,"-numberofvariables",4);
@@ -126,7 +126,7 @@ nm = neldermead_configure(nm,"-maxiter",200);
 nm = neldermead_configure(nm,"-maxfunevals",300);
 nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-1);
 nm = neldermead_configure(nm,"-simplex0method","axes");
-nm = neldermead_configure(nm,"-method","box");
+nm = neldermead_configure(nm,"-method","nmconstraints");
 nm = neldermead_configure(nm,"-nbineqconst",3);
 //nm = neldermead_configure(nm,"-verbose",1);
 nm = neldermead_configure(nm,"-verbosetermination",1);
@@ -143,7 +143,7 @@ status = neldermead_get(nm,"-status");
 assert_equal ( status , "maxfuneval" );
 nm = neldermead_destroy(nm);
 //
-// Test with Box algorithm and default axes initial simplex
+// Test with NM-constrained algorithm and default axes initial simplex
 // Add bounds and simplex initial length so that there is a need 
 // for variable projection.
 //
@@ -155,7 +155,7 @@ nm = neldermead_configure(nm,"-maxiter",200);
 nm = neldermead_configure(nm,"-maxfunevals",1000);
 nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-4);
 nm = neldermead_configure(nm,"-simplex0method","axes");
-nm = neldermead_configure(nm,"-method","box");
+nm = neldermead_configure(nm,"-method","nmconstraints");
 nm = neldermead_configure(nm,"-nbineqconst",3);
 //nm = neldermead_configure(nm,"-verbose",1);
 nm = neldermead_configure(nm,"-verbosetermination",1);
@@ -174,7 +174,7 @@ status = neldermead_get(nm,"-status");
 assert_equal ( status , "tolsize" );
 nm = neldermead_destroy(nm);
 //
-// Test with Box algorithm and default axes initial simplex
+// Test with NM-constrained algorithm and default axes initial simplex
 // Add bounds and simplex initial length so that there is a need 
 // for variable projection.
 // Here the initial simplex is computed with Box randomized bounds method
@@ -183,11 +183,11 @@ nm = neldermead_destroy(nm);
 // The convergence is not accurate in this case, whatever the 
 // value of the relative tolerance on simplex size.
 //
-    //
-    // Initialize the random number generator, so that the results are always the
-    // same.
-    //
-    rand("seed" , 0)
+//
+// Initialize the random number generator, so that the results are always the
+// same.
+//
+rand("seed" , 0)
 nm = neldermead_new ();
 nm = neldermead_configure(nm,"-numberofvariables",4);
 nm = neldermead_configure(nm,"-function",optimtestcase);
@@ -196,7 +196,7 @@ nm = neldermead_configure(nm,"-maxiter",300);
 nm = neldermead_configure(nm,"-maxfunevals",1000);
 nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-8);
 nm = neldermead_configure(nm,"-simplex0method","axes");
-nm = neldermead_configure(nm,"-method","box");
+nm = neldermead_configure(nm,"-method","nmconstraints");
 nm = neldermead_configure(nm,"-nbineqconst",3);
 //nm = neldermead_configure(nm,"-verbose",1);
 nm = neldermead_configure(nm,"-verbosetermination",1);
@@ -220,17 +220,17 @@ nbve = optimsimplex_getnbve ( simplexopt );
 assert_equal ( nbve , 8 );
 nm = neldermead_destroy(nm);
 //
-// Test with Box algorithm and default axes initial simplex
+// Test with NM-constrained algorithm and default axes initial simplex
 // Add bounds and simplex initial length so that there is a need 
 // for variable projection.
 // Here the initial simplex is computed with Box randomized bounds method
 // and user-defined number of points in the simplex, i.e. 6
 //
-    //
-    // Initialize the random number generator, so that the results are always the
-    // same.
-    //
-    rand("seed" , 0)
+//
+// Initialize the random number generator, so that the results are always the
+// same.
+//
+rand("seed" , 0)
 nm = neldermead_new ();
 nm = neldermead_configure(nm,"-numberofvariables",4);
 nm = neldermead_configure(nm,"-function",optimtestcase);
@@ -239,7 +239,7 @@ nm = neldermead_configure(nm,"-maxiter",300);
 nm = neldermead_configure(nm,"-maxfunevals",1000);
 nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-6);
 nm = neldermead_configure(nm,"-simplex0method","axes");
-nm = neldermead_configure(nm,"-method","box");
+nm = neldermead_configure(nm,"-method","nmconstraints");
 nm = neldermead_configure(nm,"-nbineqconst",3);
 //nm = neldermead_configure(nm,"-verbose",1);
 nm = neldermead_configure(nm,"-verbosetermination",1);
@@ -264,7 +264,7 @@ nbve = optimsimplex_getnbve ( simplexopt );
 assert_equal ( nbve , 6 );
 nm = neldermead_destroy(nm);
 //
-// Test with Box algorithm and default axes initial simplex
+// Test with NM-constrained algorithm and default axes initial simplex
 // Add bounds and simplex initial length so that there is a need 
 // for variable projection.
 // Here the initial simplex is user-defined.
@@ -279,7 +279,7 @@ nm = neldermead_configure(nm,"-maxiter",300);
 nm = neldermead_configure(nm,"-maxfunevals",1000);
 nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-3);
 nm = neldermead_configure(nm,"-simplex0method","axes");
-nm = neldermead_configure(nm,"-method","box");
+nm = neldermead_configure(nm,"-method","nmconstraints");
 nm = neldermead_configure(nm,"-nbineqconst",3);
 //nm = neldermead_configure(nm,"-verbose",1);
 nm = neldermead_configure(nm,"-verbosetermination",1);
@@ -311,3 +311,527 @@ simplexopt = neldermead_get ( nm , "-simplexopt" );
 nbve = optimsimplex_getnbve ( simplexopt );
 assert_equal ( nbve , 7 );
 nm = neldermead_destroy(nm);
+//
+// Test with NM-constrained algorithm and default axes initial simplex
+// Test that verbose mode works fine.
+//
+nm = neldermead_new ();
+nm = neldermead_configure(nm,"-numberofvariables",4);
+nm = neldermead_configure(nm,"-function",optimtestcase);
+nm = neldermead_configure(nm,"-x0",[0.0 0.0 0.0 0.0]');
+nm = neldermead_configure(nm,"-maxiter",5);
+nm = neldermead_configure(nm,"-maxfunevals",1000);
+nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-3);
+nm = neldermead_configure(nm,"-simplex0method","axes");
+nm = neldermead_configure(nm,"-method","nmconstraints");
+nm = neldermead_configure(nm,"-nbineqconst",3);
+nm = neldermead_configure(nm,"-verbose",1);
+nm = neldermead_configure(nm,"-verbosetermination",1);
+nm = neldermead_configure(nm,"-boundsmin",[-10.0 -10.0 -10.0 -10.0]);
+nm = neldermead_configure(nm,"-boundsmax",[10.0 10.0 10.0 10.0]);
+nm = neldermead_configure(nm,"-simplex0method","randbounds");
+nm = neldermead_configure(nm,"-coords0",coords);
+nm = neldermead_search(nm);
+Function Evaluation #1 is [0] at [0 0 0 0]
+
+Function Evaluation #2 is [0] at [0 0 0 0]
+
+Function Evaluation #3 is [8 10 5] at [0 0 0 0]
+
+Function Evaluation #4 is [0 8 10 5] at [0 0 0 0]
+
+Function Evaluation #1 is [0] at [0 0 0 0]
+
+Function Evaluation #2 is [238.57187] at [-3.8478185 8.6592324 -5.7079843 -3.7471601]
+
+Function Evaluation #3 is [32.861872] at [-2.767278 -4.1554667 1.3284976 -0.3470561]
+
+Function Evaluation #4 is [14.30847] at [-3.3565622 1.8701894 0.0306832 -1.2628248]
+
+Function Evaluation #5 is [213.76968] at [-4.6137504 2.6514897 -1.896092 8.3694157]
+
+Function Evaluation #6 is [265.53357] at [-9.1253313 -0.3629821 -4.720888 -1.7037926]
+
+Function Evaluation #7 is [71.673505] at [-4.387004 -7.4398831 5.5662572 -5.7619391]
+
+Function Evaluation #8 is [379.26278] at [-7.7572907 3.7137919 -6.9375666 3.9417012]
+
+Scaling initial simplex into nonlinear inequality constraints...
+
+Scaling vertex #2/8 at [-3.8478185 8.6592324 -5.7079843 -3.7471601]... 
+
+Function Evaluation #9 is [-113.94243 -223.0288 -119.5671] at [-3.8478185 8.6592324 -5.7079843 -3.7471601]
+
+Inequality constraint #1/3 is not satisfied for x
+
+Scaling inequality constraint with alpha = 1.000000e+000
+
+Function Evaluation #10 is [-18.868639 -50.155945 -22.989847] at [-1.9239093 4.3296162 -2.8539921 -1.87358]
+
+Inequality constraint #1/3 is not satisfied for x
+
+Scaling inequality constraint with alpha = 5.000000e-001
+
+Function Evaluation #11 is [3.0913246 -5.9883585 -0.4214980] at [-0.9619546 2.1648081 -1.4269961 -0.9367900]
+
+Inequality constraint #2/3 is not satisfied for x
+
+Scaling inequality constraint with alpha = 2.500000e-001
+
+Function Evaluation #12 is [7.6770734 5.5282242 4.4326073] at [-0.4809773 1.0824041 -0.7134980 -0.4683950]
+
+ > After scaling into inequality constraints p = [-0.4809773 1.0824041 -0.7134980 -0.4683950]
+
+Function Evaluation #13 is [11.33805] at [-0.4809773 1.0824041 -0.7134980 -0.4683950]
+
+Scaling vertex #3/8 at [-2.767278 -4.1554667 1.3284976 -0.3470561]... 
+
+Function Evaluation #14 is [-21.874827 -37.313771 -28.316431] at [-2.767278 -4.1554667 1.3284976 -0.3470561]
+
+Inequality constraint #1/3 is not satisfied for x
+
+Scaling inequality constraint with alpha = 1.000000e+000
+
+Function Evaluation #15 is [-0.2346425 -2.6070262 -3.0710995] at [-1.383639 -2.0777334 0.6642488 -0.1735280]
+
+Inequality constraint #1/3 is not satisfied for x
+
+Scaling inequality constraint with alpha = 5.000000e-001
+
+Function Evaluation #16 is [5.5583716 6.4589517 3.1112293] at [-0.6918195 -1.0388667 0.3321244 -0.0867640]
+
+ > After scaling into inequality constraints p = [-0.6918195 -1.0388667 0.3321244 -0.0867640]
+
+Function Evaluation #17 is [2.8574697] at [-0.6918195 -1.0388667 0.3321244 -0.0867640]
+
+Scaling vertex #4/8 at [-3.3565622 1.8701894 0.0306832 -1.2628248]... 
+
+Function Evaluation #18 is [-4.4265425 -16.071508 -13.71108] at [-3.3565622 1.8701894 0.0306832 -1.2628248]
+
+Inequality constraint #1/3 is not satisfied for x
+
+Scaling inequality constraint with alpha = 1.000000e+000
+
+Function Evaluation #19 is [5.8766753 2.3272762 2.1523522] at [-1.6782811 0.9350947 0.0153416 -0.6314124]
+
+ > After scaling into inequality constraints p = [-1.6782811 0.9350947 0.0153416 -0.6314124]
+
+Function Evaluation #20 is [3.0640533] at [-1.6782811 0.9350947 0.0153416 -0.6314124]
+
+Scaling vertex #5/8 at [-4.6137504 2.6514897 -1.896092 8.3694157]... 
+
+Function Evaluation #21 is [-76.428626 -165.28122 -27.950542] at [-4.6137504 2.6514897 -1.896092 8.3694157]
+
+Inequality constraint #1/3 is not satisfied for x
+
+Scaling inequality constraint with alpha = 1.000000e+000
+
+Function Evaluation #22 is [-8.7244695 -32.88139 1.8244661] at [-2.3068752 1.3257449 -0.9480460 4.1847078]
+
+Inequality constraint #1/3 is not satisfied for x
+
+Scaling inequality constraint with alpha = 5.000000e-001
+
+Function Evaluation #23 is [6.0102261 -0.2508893 6.7371673] at [-1.1534376 0.6628724 -0.4740230 2.0923539]
+
+Inequality constraint #2/3 is not satisfied for x
+
+Scaling inequality constraint with alpha = 2.500000e-001
+
+Function Evaluation #24 is [8.5982283 7.6720068 6.6998172] at [-0.5767188 0.3314362 -0.2370115 1.046177]
+
+ > After scaling into inequality constraints p = [-0.5767188 0.3314362 -0.2370115 1.046177]
+
+Function Evaluation #25 is [15.176183] at [-0.5767188 0.3314362 -0.2370115 1.046177]
+
+Scaling vertex #6/8 at [-9.1253313 -0.3629821 -4.720888 -1.7037926]... 
+
+Function Evaluation #26 is [-88.813676 -112.45691 -167.778] at [-9.1253313 -0.3629821 -4.720888 -1.7037926]
+
+Inequality constraint #1/3 is not satisfied for x
+
+Scaling inequality constraint with alpha = 1.000000e+000
+
+Function Evaluation #27 is [-13.258558 -23.321508 -34.148527] at [-4.5626657 -0.1814911 -2.360444 -0.8518963]
+
+Inequality constraint #1/3 is not satisfied for x
+
+Scaling inequality constraint with alpha = 5.000000e-001
+
+Function Evaluation #28 is [4.1577911 0.3159824 -2.7641457] at [-2.2813328 -0.0907455 -1.180222 -0.4259481]
+
+Inequality constraint #3/3 is not satisfied for x
+
+Scaling inequality constraint with alpha = 2.500000e-001
+
+Function Evaluation #29 is [7.7756631 6.9021754 4.0704566] at [-1.1406664 -0.0453728 -0.590111 -0.2129741]
+
+ > After scaling into inequality constraints p = [-1.1406664 -0.0453728 -0.590111 -0.2129741]
+
+Function Evaluation #30 is [18.876707] at [-1.1406664 -0.0453728 -0.590111 -0.2129741]
+
+Scaling vertex #7/8 at [-4.387004 -7.4398831 5.5662572 -5.7619391]... 
+
+Function Evaluation #31 is [-145.1619 -227.48157 -124.2545] at [-4.387004 -7.4398831 5.5662572 -5.7619391]
+
+Inequality constraint #1/3 is not satisfied for x
+
+Scaling inequality constraint with alpha = 1.000000e+000
+
+Function Evaluation #32 is [-33.885744 -51.907628 -28.420579] at [-2.193502 -3.7199415 2.7831286 -2.8809696]
+
+Inequality constraint #1/3 is not satisfied for x
+
+Scaling inequality constraint with alpha = 5.000000e-001
+
+Function Evaluation #33 is [-4.2690704 -6.745525 -3.9086215] at [-1.096751 -1.8599708 1.3915643 -1.4404848]
+
+Inequality constraint #1/3 is not satisfied for x
+
+Scaling inequality constraint with alpha = 2.500000e-001
+
+Function Evaluation #34 is [4.0339152 5.1793098 2.4961062] at [-0.5483755 -0.9299854 0.6957821 -0.7202424]
+
+ > After scaling into inequality constraints p = [-0.5483755 -0.9299854 0.6957821 -0.7202424]
+
+Function Evaluation #35 is [-9.6087543] at [-0.5483755 -0.9299854 0.6957821 -0.7202424]
+
+Scaling vertex #8/8 at [-7.7572907 3.7137919 -6.9375666 3.9417012]... 
+
+Function Evaluation #36 is [-107.2843 -160.7795 -154.10312] at [-7.7572907 3.7137919 -6.9375666 3.9417012]
+
+Inequality constraint #1/3 is not satisfied for x
+
+Scaling inequality constraint with alpha = 1.000000e+000
+
+Function Evaluation #37 is [-15.233487 -33.648771 -28.983262] at [-3.8786453 1.856896 -3.4687833 1.9708506]
+
+Inequality constraint #1/3 is not satisfied for x
+
+Scaling inequality constraint with alpha = 5.000000e-001
+
+Function Evaluation #38 is [4.9854221 -1.3891415 -0.5995563] at [-1.9393227 0.9284480 -1.7343917 0.9854253]
+
+Inequality constraint #2/3 is not satisfied for x
+
+Scaling inequality constraint with alpha = 2.500000e-001
+
+Function Evaluation #39 is [8.6432524 6.9142403 5.0482406] at [-0.9696613 0.4642240 -0.8671958 0.4927127]
+
+ > After scaling into inequality constraints p = [-0.9696613 0.4642240 -0.8671958 0.4927127]
+
+Function Evaluation #40 is [27.089858] at [-0.9696613 0.4642240 -0.8671958 0.4927127]
+
+Step #1 : order
+
+=================================================================
+
+Iteration #1 (total = 1)
+
+Function Eval #40
+
+Xopt : [-0.5483755 -0.9299854 0.6957821 -0.7202424]
+
+Fopt : -9.608754e+000
+
+DeltaFv : 3.669861e+001
+
+Center : [-0.7608125 0.0998668 -0.1705710 -0.0726123]
+
+Size : 2.470577e+000
+
+Vertex #1/8 : fv=-9.608754e+000, x=-5.483755e-001 -9.299854e-001 6.957821e-001 -7.202424e-001
+
+Vertex #2/8 : fv=0.000000e+000, x=0.000000e+000 0.000000e+000 0.000000e+000 0.000000e+000
+
+Vertex #3/8 : fv=2.857470e+000, x=-6.918195e-001 -1.038867e+000 3.321244e-001 -8.676401e-002
+
+Vertex #4/8 : fv=3.064053e+000, x=-1.678281e+000 9.350947e-001 1.534160e-002 -6.314124e-001
+
+Vertex #5/8 : fv=1.133805e+001, x=-4.809773e-001 1.082404e+000 -7.134980e-001 -4.683950e-001
+
+Vertex #6/8 : fv=1.517618e+001, x=-5.767188e-001 3.314362e-001 -2.370115e-001 1.046177e+000
+
+Vertex #7/8 : fv=1.887671e+001, x=-1.140666e+000 -4.537277e-002 -5.901110e-001 -2.129741e-001
+
+Vertex #8/8 : fv=2.708986e+001, x=-9.696613e-001 4.642240e-001 -8.671958e-001 4.927127e-001
+
+Reflect
+
+xbar=[-0.7309769 0.0478157 -0.0710532 -0.1533730]
+
+ > After projection into bounds p = [-0.4922926 -0.3685925 0.7250894 -0.7994586]
+
+Function Evaluation #41 is [5.0560508 6.3901531 3.6702149] at [-0.4922926 -0.3685925 0.7250894 -0.7994586]
+
+ > After scaling into inequality constraints p = [-0.4922926 -0.3685925 0.7250894 -0.7994586]
+
+Function Evaluation #42 is [-14.449807] at [-0.4922926 -0.3685925 0.7250894 -0.7994586]
+
+xr=[-0.4922926 -0.3685925 0.7250894 -0.7994586], f(xr)=-14.449807
+
+Expand
+
+ > After projection into bounds p = [-0.2536082 -0.7850008 1.5212321 -1.4455443]
+
+Function Evaluation #43 is [-0.5824575 0.5107345 0.2176639] at [-0.2536082 -0.7850008 1.5212321 -1.4455443]
+
+Inequality constraint #1/3 is not satisfied for x
+
+Scaling inequality constraint with alpha = 1.000000e+000
+
+Function Evaluation #44 is [5.0560508 6.3901531 3.6702149] at [-0.4922926 -0.3685925 0.7250894 -0.7994586]
+
+ > After scaling into inequality constraints p = [-0.4922926 -0.3685925 0.7250894 -0.7994586]
+
+Function Evaluation #45 is [-14.449807] at [-0.4922926 -0.3685925 0.7250894 -0.7994586]
+
+xe=[-0.4922926 -0.3685925 0.7250894 -0.7994586], f(xe)=-14.449807
+
+  > Perform reflection
+
+Sort
+
+=================================================================
+
+Iteration #2 (total = 2)
+
+Function Eval #45
+
+Xopt : [-0.4922926 -0.3685925 0.7250894 -0.7994586]
+
+Fopt : -1.444981e+001
+
+DeltaFv : 3.332651e+001
+
+Center : [-0.7011414 -0.0042353 0.0284646 -0.2341337]
+
+Size : 2.197539e+000
+
+Vertex #1/8 : fv=-1.444981e+001, x=-4.922926e-001 -3.685925e-001 7.250894e-001 -7.994586e-001
+
+Vertex #2/8 : fv=-9.608754e+000, x=-5.483755e-001 -9.299854e-001 6.957821e-001 -7.202424e-001
+
+Vertex #3/8 : fv=0.000000e+000, x=0.000000e+000 0.000000e+000 0.000000e+000 0.000000e+000
+
+Vertex #4/8 : fv=2.857470e+000, x=-6.918195e-001 -1.038867e+000 3.321244e-001 -8.676401e-002
+
+Vertex #5/8 : fv=3.064053e+000, x=-1.678281e+000 9.350947e-001 1.534160e-002 -6.314124e-001
+
+Vertex #6/8 : fv=1.133805e+001, x=-4.809773e-001 1.082404e+000 -7.134980e-001 -4.683950e-001
+
+Vertex #7/8 : fv=1.517618e+001, x=-5.767188e-001 3.314362e-001 -2.370115e-001 1.046177e+000
+
+Vertex #8/8 : fv=1.887671e+001, x=-1.140666e+000 -4.537277e-002 -5.901110e-001 -2.129741e-001
+
+  > Termination ?
+
+  > iterations=2 >= maxiter=5
+
+  > funevals=45 >= maxfunevals=1000
+
+  > e(x)=2.830234e-001 < tolxrelative = 2.220446e-016 * 7.397608e-001 + 0.000000e+000
+
+  > simplex size=2.197539e+000 < 0.000000e+000 + 1.000000e-003 * 2.022362e+000
+
+Reflect
+
+xbar=[-0.6383521 0.0016415 0.1168326 -0.2371565]
+
+ > After projection into bounds p = [-0.1360378 0.0486557 0.8237762 -0.2613389]
+
+Function Evaluation #46 is [6.3317996 8.764179 4.3414053] at [-0.1360378 0.0486557 0.8237762 -0.2613389]
+
+ > After scaling into inequality constraints p = [-0.1360378 0.0486557 0.8237762 -0.2613389]
+
+Function Evaluation #47 is [-17.245376] at [-0.1360378 0.0486557 0.8237762 -0.2613389]
+
+xr=[-0.1360378 0.0486557 0.8237762 -0.2613389], f(xr)=-17.245376
+
+Expand
+
+ > After projection into bounds p = [0.3662765 0.0956700 1.5307197 -0.2855214]
+
+Function Evaluation #48 is [3.3452158 7.4221434 1.457023] at [0.3662765 0.0956700 1.5307197 -0.2855214]
+
+ > After scaling into inequality constraints p = [0.3662765 0.0956700 1.5307197 -0.2855214]
+
+Function Evaluation #49 is [-31.542457] at [0.3662765 0.0956700 1.5307197 -0.2855214]
+
+xe=[0.3662765 0.0956700 1.5307197 -0.2855214], f(xe)=-31.542457
+
+  > Perform Expansion
+
+Sort
+
+=================================================================
+
+Iteration #3 (total = 3)
+
+Function Eval #49
+
+Xopt : [0.3662765 0.0956700 1.5307197 -0.2855214]
+
+Fopt : -3.154246e+001
+
+DeltaFv : 4.671864e+001
+
+Center : [-0.5127735 0.0133950 0.2935685 -0.2432021]
+
+Size : 2.702011e+000
+
+Vertex #1/8 : fv=-3.154246e+001, x=3.662765e-001 9.566998e-002 1.530720e+000 -2.855214e-001
+
+Vertex #2/8 : fv=-1.444981e+001, x=-4.922926e-001 -3.685925e-001 7.250894e-001 -7.994586e-001
+
+Vertex #3/8 : fv=-9.608754e+000, x=-5.483755e-001 -9.299854e-001 6.957821e-001 -7.202424e-001
+
+Vertex #4/8 : fv=0.000000e+000, x=0.000000e+000 0.000000e+000 0.000000e+000 0.000000e+000
+
+Vertex #5/8 : fv=2.857470e+000, x=-6.918195e-001 -1.038867e+000 3.321244e-001 -8.676401e-002
+
+Vertex #6/8 : fv=3.064053e+000, x=-1.678281e+000 9.350947e-001 1.534160e-002 -6.314124e-001
+
+Vertex #7/8 : fv=1.133805e+001, x=-4.809773e-001 1.082404e+000 -7.134980e-001 -4.683950e-001
+
+Vertex #8/8 : fv=1.517618e+001, x=-5.767188e-001 3.314362e-001 -2.370115e-001 1.046177e+000
+
+  > Termination ?
+
+  > iterations=3 >= maxiter=5
+
+  > funevals=49 >= maxfunevals=1000
+
+  > e(x)=3.258152e-001 < tolxrelative = 2.220446e-016 * 6.390977e-001 + 0.000000e+000
+
+  > simplex size=2.702011e+000 < 0.000000e+000 + 1.000000e-003 * 2.022362e+000
+
+Reflect
+
+xbar=[-0.5036385 -0.0320394 0.3693656 -0.4273991]
+
+ > After projection into bounds p = [-0.4305582 -0.3955150 0.9757427 -1.9009752]
+
+Function Evaluation #50 is [0.2507322 -1.0092652 2.0853594] at [-0.4305582 -0.3955150 0.9757427 -1.9009752]
+
+Inequality constraint #2/3 is not satisfied for x
+
+Scaling inequality constraint with alpha = 1.000000e+000
+
+Function Evaluation #51 is [4.3450374 4.8961397 3.6218408] at [-0.4670983 -0.2137772 0.6725542 -1.1641872]
+
+ > After scaling into inequality constraints p = [-0.4670983 -0.2137772 0.6725542 -1.1641872]
+
+Function Evaluation #52 is [-16.344698] at [-0.4670983 -0.2137772 0.6725542 -1.1641872]
+
+xr=[-0.4670983 -0.2137772 0.6725542 -1.1641872], f(xr)=-16.344698
+
+  > Perform reflection
+
+Sort
+
+=================================================================
+
+Iteration #4 (total = 4)
+
+Function Eval #52
+
+Xopt : [0.3662765 0.0956700 1.5307197 -0.2855214]
+
+Fopt : -3.154246e+001
+
+DeltaFv : 4.288051e+001
+
+Center : [-0.4990710 -0.0547566 0.4072642 -0.5194976]
+
+Size : 2.702011e+000
+
+Vertex #1/8 : fv=-3.154246e+001, x=3.662765e-001 9.566998e-002 1.530720e+000 -2.855214e-001
+
+Vertex #2/8 : fv=-1.634470e+001, x=-4.670983e-001 -2.137772e-001 6.725542e-001 -1.164187e+000
+
+Vertex #3/8 : fv=-1.444981e+001, x=-4.922926e-001 -3.685925e-001 7.250894e-001 -7.994586e-001
+
+Vertex #4/8 : fv=-9.608754e+000, x=-5.483755e-001 -9.299854e-001 6.957821e-001 -7.202424e-001
+
+Vertex #5/8 : fv=0.000000e+000, x=0.000000e+000 0.000000e+000 0.000000e+000 0.000000e+000
+
+Vertex #6/8 : fv=2.857470e+000, x=-6.918195e-001 -1.038867e+000 3.321244e-001 -8.676401e-002
+
+Vertex #7/8 : fv=3.064053e+000, x=-1.678281e+000 9.350947e-001 1.534160e-002 -6.314124e-001
+
+Vertex #8/8 : fv=1.133805e+001, x=-4.809773e-001 1.082404e+000 -7.134980e-001 -4.683950e-001
+
+  > Termination ?
+
+  > iterations=4 >= maxiter=5
+
+  > funevals=52 >= maxfunevals=1000
+
+  > e(x)=3.067545e-001 < tolxrelative = 2.220446e-016 * 8.293443e-001 + 0.000000e+000
+
+  > simplex size=2.702011e+000 < 0.000000e+000 + 1.000000e-003 * 2.022362e+000
+
+Reflect
+
+xbar=[-0.5016558 -0.2172082 0.5673731 -0.5267980]
+
+ > After projection into bounds p = [-0.5223342 -1.5168204 1.8482442 -0.5852010]
+
+Function Evaluation #53 is [-1.7599751 -0.0827833 -2.3197696] at [-0.5223342 -1.5168204 1.8482442 -0.5852010]
+
+Inequality constraint #1/3 is not satisfied for x
+
+Scaling inequality constraint with alpha = 1.000000e+000
+
+Function Evaluation #54 is [3.0993829 5.0893666 1.8661831] at [-0.5119950 -0.8670143 1.2078086 -0.5559995]
+
+ > After scaling into inequality constraints p = [-0.5119950 -0.8670143 1.2078086 -0.5559995]
+
+Function Evaluation #55 is [-18.12034] at [-0.5119950 -0.8670143 1.2078086 -0.5559995]
+
+xr=[-0.5119950 -0.8670143 1.2078086 -0.5559995], f(xr)=-18.120340
+
+  > Perform reflection
+
+Sort
+
+=================================================================
+
+Iteration #5 (total = 5)
+
+Function Eval #55
+
+Xopt : [0.3662765 0.0956700 1.5307197 -0.2855214]
+
+Fopt : -3.154246e+001
+
+DeltaFv : 3.460651e+001
+
+Center : [-0.5029482 -0.2984339 0.6474275 -0.5304482]
+
+Size : 2.702011e+000
+
+Vertex #1/8 : fv=-3.154246e+001, x=3.662765e-001 9.566998e-002 1.530720e+000 -2.855214e-001
+
+Vertex #2/8 : fv=-1.812034e+001, x=-5.119950e-001 -8.670143e-001 1.207809e+000 -5.559995e-001
+
+Vertex #3/8 : fv=-1.634470e+001, x=-4.670983e-001 -2.137772e-001 6.725542e-001 -1.164187e+000
+
+Vertex #4/8 : fv=-1.444981e+001, x=-4.922926e-001 -3.685925e-001 7.250894e-001 -7.994586e-001
+
+Vertex #5/8 : fv=-9.608754e+000, x=-5.483755e-001 -9.299854e-001 6.957821e-001 -7.202424e-001
+
+Vertex #6/8 : fv=0.000000e+000, x=0.000000e+000 0.000000e+000 0.000000e+000 0.000000e+000
+
+Vertex #7/8 : fv=2.857470e+000, x=-6.918195e-001 -1.038867e+000 3.321244e-001 -8.676401e-002
+
+Vertex #8/8 : fv=3.064053e+000, x=-1.678281e+000 9.350947e-001 1.534160e-002 -6.314124e-001
+
+  > Termination ?
+
+  > iterations=5 >= maxiter=5
+
+Terminate with status : maxiter
+
+nm = neldermead_destroy(nm);
index d2c3253..a0b9282 100644 (file)
@@ -47,6 +47,20 @@ function flag = assert_equal ( computed , expected )
 endfunction
 
 //
+//  Reference:
+//
+//    An extension of the simplex method to constrained
+//    nonlinear optimization
+//    M.B. Subrahmanyam
+//    Journal of optimization theory and applications
+//    Vol. 62, August 1989
+//
+//    Gould F.J.
+//    Nonlinear Tolerance Programming
+//    Numerical methods for Nonlinear optimization
+//    Edited by F.A. Lootsma, pp 349-366, 1972
+
+//
 // optimtestcase --
 //   Non linear inequality constraints are positive.
 //    
@@ -95,7 +109,7 @@ function result = optimtestcase ( x , index )
 endfunction
 
 //
-// Test with Box algorithm and default axes initial simplex
+// Test with NM-constrained algorithm and default axes initial simplex
 //
 nm = neldermead_new ();
 nm = neldermead_configure(nm,"-numberofvariables",4);
@@ -105,7 +119,7 @@ nm = neldermead_configure(nm,"-maxiter",200);
 nm = neldermead_configure(nm,"-maxfunevals",400);
 nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-3);
 nm = neldermead_configure(nm,"-simplex0method","axes");
-nm = neldermead_configure(nm,"-method","box");
+nm = neldermead_configure(nm,"-method","nmconstraints");
 nm = neldermead_configure(nm,"-nbineqconst",3);
 //nm = neldermead_configure(nm,"-verbose",1);
 nm = neldermead_configure(nm,"-verbosetermination",1);
@@ -122,7 +136,7 @@ assert_equal ( status , "tolsize" );
 nm = neldermead_destroy(nm);
 
 //
-// Test with Box algorithm and restart
+// Test with NM-constrained algorithm and restart
 //
 nm = neldermead_new ();
 nm = neldermead_configure(nm,"-numberofvariables",4);
@@ -132,7 +146,7 @@ nm = neldermead_configure(nm,"-maxiter",200);
 nm = neldermead_configure(nm,"-maxfunevals",300);
 nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-1);
 nm = neldermead_configure(nm,"-simplex0method","axes");
-nm = neldermead_configure(nm,"-method","box");
+nm = neldermead_configure(nm,"-method","nmconstraints");
 nm = neldermead_configure(nm,"-nbineqconst",3);
 //nm = neldermead_configure(nm,"-verbose",1);
 nm = neldermead_configure(nm,"-verbosetermination",1);
@@ -150,7 +164,7 @@ assert_equal ( status , "maxfuneval" );
 nm = neldermead_destroy(nm);
 
 //
-// Test with Box algorithm and default axes initial simplex
+// Test with NM-constrained algorithm and default axes initial simplex
 // Add bounds and simplex initial length so that there is a need 
 // for variable projection.
 //
@@ -162,7 +176,7 @@ nm = neldermead_configure(nm,"-maxiter",200);
 nm = neldermead_configure(nm,"-maxfunevals",1000);
 nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-4);
 nm = neldermead_configure(nm,"-simplex0method","axes");
-nm = neldermead_configure(nm,"-method","box");
+nm = neldermead_configure(nm,"-method","nmconstraints");
 nm = neldermead_configure(nm,"-nbineqconst",3);
 //nm = neldermead_configure(nm,"-verbose",1);
 nm = neldermead_configure(nm,"-verbosetermination",1);
@@ -181,7 +195,7 @@ status = neldermead_get(nm,"-status");
 assert_equal ( status , "tolsize" );
 nm = neldermead_destroy(nm);
 //
-// Test with Box algorithm and default axes initial simplex
+// Test with NM-constrained algorithm and default axes initial simplex
 // Add bounds and simplex initial length so that there is a need 
 // for variable projection.
 // Here the initial simplex is computed with Box randomized bounds method
@@ -190,11 +204,11 @@ nm = neldermead_destroy(nm);
 // The convergence is not accurate in this case, whatever the 
 // value of the relative tolerance on simplex size.
 //
-    //
-    // Initialize the random number generator, so that the results are always the
-    // same.
-    //
-    rand("seed" , 0)
+//
+// Initialize the random number generator, so that the results are always the
+// same.
+//
+rand("seed" , 0)
 nm = neldermead_new ();
 nm = neldermead_configure(nm,"-numberofvariables",4);
 nm = neldermead_configure(nm,"-function",optimtestcase);
@@ -203,7 +217,7 @@ nm = neldermead_configure(nm,"-maxiter",300);
 nm = neldermead_configure(nm,"-maxfunevals",1000);
 nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-8);
 nm = neldermead_configure(nm,"-simplex0method","axes");
-nm = neldermead_configure(nm,"-method","box");
+nm = neldermead_configure(nm,"-method","nmconstraints");
 nm = neldermead_configure(nm,"-nbineqconst",3);
 //nm = neldermead_configure(nm,"-verbose",1);
 nm = neldermead_configure(nm,"-verbosetermination",1);
@@ -229,17 +243,17 @@ nm = neldermead_destroy(nm);
 
 
 //
-// Test with Box algorithm and default axes initial simplex
+// Test with NM-constrained algorithm and default axes initial simplex
 // Add bounds and simplex initial length so that there is a need 
 // for variable projection.
 // Here the initial simplex is computed with Box randomized bounds method
 // and user-defined number of points in the simplex, i.e. 6
 //
-    //
-    // Initialize the random number generator, so that the results are always the
-    // same.
-    //
-    rand("seed" , 0)
+//
+// Initialize the random number generator, so that the results are always the
+// same.
+//
+rand("seed" , 0)
 nm = neldermead_new ();
 nm = neldermead_configure(nm,"-numberofvariables",4);
 nm = neldermead_configure(nm,"-function",optimtestcase);
@@ -248,7 +262,7 @@ nm = neldermead_configure(nm,"-maxiter",300);
 nm = neldermead_configure(nm,"-maxfunevals",1000);
 nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-6);
 nm = neldermead_configure(nm,"-simplex0method","axes");
-nm = neldermead_configure(nm,"-method","box");
+nm = neldermead_configure(nm,"-method","nmconstraints");
 nm = neldermead_configure(nm,"-nbineqconst",3);
 //nm = neldermead_configure(nm,"-verbose",1);
 nm = neldermead_configure(nm,"-verbosetermination",1);
@@ -273,7 +287,7 @@ nbve = optimsimplex_getnbve ( simplexopt );
 assert_equal ( nbve , 6 );
 nm = neldermead_destroy(nm);
 //
-// Test with Box algorithm and default axes initial simplex
+// Test with NM-constrained algorithm and default axes initial simplex
 // Add bounds and simplex initial length so that there is a need 
 // for variable projection.
 // Here the initial simplex is user-defined.
@@ -288,7 +302,7 @@ nm = neldermead_configure(nm,"-maxiter",300);
 nm = neldermead_configure(nm,"-maxfunevals",1000);
 nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-3);
 nm = neldermead_configure(nm,"-simplex0method","axes");
-nm = neldermead_configure(nm,"-method","box");
+nm = neldermead_configure(nm,"-method","nmconstraints");
 nm = neldermead_configure(nm,"-nbineqconst",3);
 //nm = neldermead_configure(nm,"-verbose",1);
 nm = neldermead_configure(nm,"-verbosetermination",1);
@@ -321,3 +335,26 @@ nbve = optimsimplex_getnbve ( simplexopt );
 assert_equal ( nbve , 7 );
 nm = neldermead_destroy(nm);
 
+//
+// Test with NM-constrained algorithm and default axes initial simplex
+// Test that verbose mode works fine.
+//
+nm = neldermead_new ();
+nm = neldermead_configure(nm,"-numberofvariables",4);
+nm = neldermead_configure(nm,"-function",optimtestcase);
+nm = neldermead_configure(nm,"-x0",[0.0 0.0 0.0 0.0]');
+nm = neldermead_configure(nm,"-maxiter",5);
+nm = neldermead_configure(nm,"-maxfunevals",1000);
+nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-3);
+nm = neldermead_configure(nm,"-simplex0method","axes");
+nm = neldermead_configure(nm,"-method","nmconstraints");
+nm = neldermead_configure(nm,"-nbineqconst",3);
+nm = neldermead_configure(nm,"-verbose",1);
+nm = neldermead_configure(nm,"-verbosetermination",1);
+nm = neldermead_configure(nm,"-boundsmin",[-10.0 -10.0 -10.0 -10.0]);
+nm = neldermead_configure(nm,"-boundsmax",[10.0 10.0 10.0 10.0]);
+nm = neldermead_configure(nm,"-simplex0method","randbounds");
+nm = neldermead_configure(nm,"-coords0",coords);
+nm = neldermead_search(nm);
+nm = neldermead_destroy(nm);
+
index 6cfba0c..52d789b 100644 (file)
@@ -100,7 +100,7 @@ nm = neldermead_configure(nm,"-maxiter",200);
 nm = neldermead_configure(nm,"-maxfunevals",1000);
 nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-4);
 nm = neldermead_configure(nm,"-simplex0method","axes");
-nm = neldermead_configure(nm,"-method","box");
+nm = neldermead_configure(nm,"-method","nmconstraints");
 nm = neldermead_configure(nm,"-nbineqconst",3);
 //nm = neldermead_configure(nm,"-verbose",1);
 nm = neldermead_configure(nm,"-verbosetermination",1);
index 4ebbc1d..3afc222 100644 (file)
@@ -106,7 +106,7 @@ nm = neldermead_configure(nm,"-maxiter",200);
 nm = neldermead_configure(nm,"-maxfunevals",1000);
 nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-4);
 nm = neldermead_configure(nm,"-simplex0method","axes");
-nm = neldermead_configure(nm,"-method","box");
+nm = neldermead_configure(nm,"-method","nmconstraints");
 nm = neldermead_configure(nm,"-nbineqconst",3);
 //nm = neldermead_configure(nm,"-verbose",1);
 nm = neldermead_configure(nm,"-verbosetermination",1);
@@ -157,3 +157,21 @@ assert_equal ( computed , expected );
 //
 nm = neldermead_destroy(nm);
 
+
+//
+// Test search with verbose to log file
+//
+nm = neldermead_new ();
+nm = neldermead_configure(nm,"-numberofvariables",4);
+nm = neldermead_configure(nm,"-function",optimtestcase);
+nm = neldermead_configure(nm,"-x0",[0.0 0.0 0.0 0.0]');
+nm = neldermead_configure(nm,"-maxiter",10);
+nm = neldermead_configure(nm,"-verbose",1);
+nm = neldermead_configure(nm,"-logfile" , "search.txt" );
+nm = neldermead_configure(nm,"-verbosetermination",1);
+nm = neldermead_search(nm);
+nm = neldermead_destroy(nm);
+computed = deletefile("search.txt");
+assert_equal ( computed , %t );
+
+
index b4b5d86..efffa2e 100644 (file)
@@ -51,7 +51,7 @@ endfunction
 // So the actual name "mydata" does not matter
 // and whatever variable name can be used.
 //
-function y = rosenbrock ( x , mydata )
+function [ y , mydata ] = rosenbrock ( x , mydata )
   a = mydata.a
   y = 100*(x(2)-x(1)^2)^2 + ( a - x(1))^2;
 endfunction
@@ -66,21 +66,21 @@ nm = neldermead_new ();
 nm = neldermead_configure(nm,"-numberofvariables",2);
 nm = neldermead_configure(nm,"-function",rosenbrock);
 nm = neldermead_configure(nm,"-x0",[-1.2 1.0]');
-nm = neldermead_configure(nm,"-maxiter",200);
-nm = neldermead_configure(nm,"-maxfunevals",300);
+nm = neldermead_configure(nm,"-maxiter",400);
+nm = neldermead_configure(nm,"-maxfunevals",800);
 nm = neldermead_configure(nm,"-tolfunrelative",10*%eps);
 nm = neldermead_configure(nm,"-tolxrelative",10*%eps);
 nm = neldermead_configure(nm,"-simplex0method","axes");
 nm = neldermead_configure(nm,"-simplex0length",1.0);
 nm = neldermead_configure(nm,"-method","variable");
-nm = neldermead_configure(nm,"-verbose",0);
+nm = neldermead_configure(nm,"-verbose",1);
 nm = neldermead_configure(nm,"-verbosetermination",0);
 nm = neldermead_configure(nm,"-storehistory",1);
 nm = neldermead_configure(nm,"-costfargument",mystuff);
 nm = neldermead_search(nm);
 // Check optimum point
 xopt = neldermead_get(nm,"-xopt");
-assert_close ( xopt , [10.294738520749856292014 105.94430388598999570604]', 1e-6 );
+assert_close ( xopt , [12.0 144.0]', 1e-6 );
 nm = neldermead_destroy(nm);
 
 
index b3faabc..e669a75 100644 (file)
@@ -288,3 +288,25 @@ assert_close ( f , [32745.827    261254.17    96991.969    197008.03    130368.4
 assert_close ( f , [-2351243.5    32745.827    261254.17    96991.969    197008.03    130368.43    146831.57] , 1.e-7 );
 opt = optimbase_destroy(opt);
 
+//
+// Test with an additional argument
+// which is modified by the call.
+//
+function [ y , mydata ] = rosenbrock3 ( x , mydata )
+  mydata.n = mydata.n + 1
+  y = 100*(x(2)-x(1)^2)^2 + ( 1.0 - x(1))^2;
+endfunction
+
+mystuff = tlist(["T_MYSTUFF","n"]);
+mystuff.n = 0;
+
+opt = optimbase_new ();
+opt = optimbase_configure(opt,"-numberofvariables",2);
+opt = optimbase_configure(opt,"-function",rosenbrock3);
+opt = optimbase_configure(opt,"-costfargument",mystuff);
+[ opt , f ] = optimbase_function ( opt , [0.0 0.0] );
+assert_close ( f , 1.0 , %eps );
+mystuff = optimbase_cget ( opt , "-costfargument" );
+assert_equal ( mystuff.n , 1 );
+opt = optimbase_destroy(opt);
+
diff --git a/scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_log.tst b/scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_log.tst
new file mode 100644 (file)
index 0000000..0a4d996
--- /dev/null
@@ -0,0 +1,99 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008-2009 - INRIA - Michael Baudin
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+
+
+//
+// assert_close --
+//   Returns 1 if the two real matrices computed and expected are close,
+//   i.e. if the relative distance between computed and expected is lesser than epsilon.
+// Arguments
+//   computed, expected : the two matrices to compare
+//   epsilon : a small number
+//
+function flag = assert_close ( computed, expected, epsilon )
+  if expected==0.0 then
+    shift = norm(computed-expected);
+  else
+    shift = norm(computed-expected)/norm(expected);
+  end
+  if shift < epsilon then
+    flag = 1;
+  else
+    flag = 0;
+  end
+  if flag <> 1 then pause,end
+endfunction
+//
+// assert_equal --
+//   Returns 1 if the two real matrices computed and expected are equal.
+// Arguments
+//   computed, expected : the two matrices to compare
+//   epsilon : a small number
+//
+function flag = assert_equal ( computed , expected )
+  if computed==expected then
+    flag = 1;
+  else
+    flag = 0;
+  end
+  if flag <> 1 then pause,end
+endfunction
+function y = rosenbrock (x)
+  y = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;
+endfunction
+
+//
+// Test basic verbose
+//
+opt = optimbase_new ();
+computed = optimbase_cget(opt,"-verbose");
+assert_equal ( computed , 0 );
+opt = optimbase_configure(opt,"-verbose",1);
+computed = optimbase_cget(opt,"-verbose");
+assert_equal ( computed , 1 );
+opt = optimbase_log ( opt , "My interesting message" );
+opt = optimbase_configure(opt,"-verbose",0);
+// Log a message relative to the stopping rule
+computed = optimbase_cget(opt,"-verbosetermination");
+assert_equal ( computed , 0 );
+opt = optimbase_configure(opt,"-verbosetermination",1);
+computed = optimbase_cget(opt,"-verbosetermination");
+assert_equal ( computed , 1 );
+opt = optimbase_stoplog ( opt , "My interesting stop message" );
+opt = optimbase_configure(opt,"-verbosetermination",0);
+opt = optimbase_destroy(opt);
+
+//
+// Test logging into a file
+//
+opt = optimbase_new ();
+opt = optimbase_configure(opt,"-verbose",1);
+computed = optimbase_cget(opt,"-logfile");
+assert_equal ( computed , "" );
+opt = optimbase_configure(opt,"-logfile","mylogfile.txt");
+computed = optimbase_cget(opt,"-logfile");
+assert_equal ( computed , "mylogfile.txt" );
+opt = optimbase_log ( opt , "My interesting message" );
+opt = optimbase_configure(opt,"-verbose",0);
+opt = optimbase_log ( opt , "My NOT interesting message" );
+opt = optimbase_destroy(opt);
+computed = fileinfo ( "mylogfile.txt" );
+assert_equal ( ( computed <> [] ) , %t );
+// Check content
+fd = mopen( "mylogfile.txt" , "r" );
+computed = mgetl ( fd , 1 ); // The date is not checked
+computed = mgetl ( fd , 1 );
+assert_equal ( computed , "My interesting message" );
+computed = meof ( fd );
+assert_equal ( computed , 0 );
+mclose ( fd );
+computed = deletefile ( "mylogfile.txt" );
+assert_equal ( computed , %t );
+
+