Deleted vectorized computation feature. Deleted neldermead_contour. Fixed the demos.
Michaƫl Baudin [Fri, 25 Sep 2009 09:05:44 +0000 (11:05 +0200)]
26 files changed:
scilab/modules/optimization/demos/neldermead/mckinnon.sci [deleted file]
scilab/modules/optimization/demos/neldermead/neldermead.dem.sce [deleted file]
scilab/modules/optimization/demos/neldermead/neldermead_boxproblemA.sce [moved from scilab/modules/optimization/demos/neldermead/nmplot_boxproblemA.sce with 50% similarity]
scilab/modules/optimization/demos/neldermead/neldermead_rosenbrock.sce
scilab/modules/optimization/demos/neldermead/nmplot_han1.sce
scilab/modules/optimization/demos/neldermead/nmplot_han2.sce
scilab/modules/optimization/demos/neldermead/nmplot_mckinnon.sce
scilab/modules/optimization/demos/neldermead/nmplot_mckinnon2.sce
scilab/modules/optimization/demos/neldermead/nmplot_quadratic.fixed.sce
scilab/modules/optimization/demos/neldermead/nmplot_quadratic.fixed2.sce [new file with mode: 0644]
scilab/modules/optimization/demos/neldermead/nmplot_rosenbrock.fixed.sce
scilab/modules/optimization/demos/neldermead/nmplot_rosenbrock.sce
scilab/modules/optimization/demos/optimization.dem.gateway.sce
scilab/modules/optimization/help/en_US/neldermead/neldermead.xml
scilab/modules/optimization/help/en_US/neldermead/nmplot.xml
scilab/modules/optimization/help/en_US/optimsimplex/optimsimplex.xml
scilab/modules/optimization/macros/neldermead/neldermead_cget.sci
scilab/modules/optimization/macros/neldermead/neldermead_configure.sci
scilab/modules/optimization/macros/neldermead/neldermead_contour.sci [deleted file]
scilab/modules/optimization/macros/neldermead/neldermead_costf.sci
scilab/modules/optimization/macros/neldermead/neldermead_display.sci
scilab/modules/optimization/macros/neldermead/neldermead_new.sci
scilab/modules/optimization/macros/neldermead/neldermead_search.sci
scilab/modules/optimization/macros/neldermead/nmplot_contour.sci
scilab/modules/optimization/macros/optimsimplex/optimsimplex_print.sci
scilab/modules/optimization/macros/optimsimplex/optimsimplex_tostring.sci

diff --git a/scilab/modules/optimization/demos/neldermead/mckinnon.sci b/scilab/modules/optimization/demos/neldermead/mckinnon.sci
deleted file mode 100644 (file)
index 082cbc2..0000000
+++ /dev/null
@@ -1,112 +0,0 @@
-//% MCKINNON computes the McKinnon function.
-//
-//  Discussion:
-//
-//    This function has a global minimizer:
-//
-//      X* = ( 0.0, -0.5 ), F(X*) = -0.25
-//
-//    There are three parameters, TAU, THETA and PHI.
-//
-//    1 < TAU, then F is strictly convex.
-//             and F has continuous first derivatives.
-//    2 < TAU, then F has continuous second derivatives.
-//    3 < TAU, then F has continuous third derivatives.
-//
-//    However, this function can cause the Nelder-Mead optimization
-//    algorithm to "converge" to a point which is not the minimizer
-//    of the function F.
-//
-//    Sample parameter values which cause problems for Nelder-Mead 
-//    include:
-//
-//      TAU = 1, THETA = 15, PHI =  10;
-//      TAU = 2, THETA =  6, PHI =  60;
-//      TAU = 3, THETA =  6, PHI = 400;
-//
-//    To get the bad behavior, we also assume the initial simplex has the form
-//
-//      X1 = (0,0),
-//      X2 = (1,1),
-//      X3 = (A,B), 
-//
-//    where 
-//
-//      A = (1+sqrt(33))/8 =  0.84307...
-//      B = (1-sqrt(33))/8 = -0.59307...
-//
-//  Licensing:
-//
-//    This code is distributed under the GNU LGPL license.
-//
-//  Modified:
-//
-//    09 February 2008
-//
-//  Author:
-//
-//    John Burkardt
-//
-//  Reference:
-//
-//    Ken McKinnon,
-//    Convergence of the Nelder-Mead simplex method to a nonstationary point,
-//    SIAM Journal on Optimization,
-//    Volume 9, Number 1, 1998, pages 148-158.
-//
-//  Parameters:
-//
-//    Input, real X(2), the argument of the function.
-//
-//    Output, real F, the value of the function at X.
-//
-// Copyright (C) 2009 - INRIA - Michael Baudin, Scilab port
-function f = mckinnon1 ( x )
-
-  if ( length ( x ) ~= 2 )
-    error ( 'Error: function expects a two dimensional input\n' );
-  end
-
-  tau = 1.0;
-  theta = 15.0;
-  phi = 10.0;
-
-  if ( x(1) <= 0.0 )
-    f = theta * phi * abs ( x(1) ).^tau + x(2) * ( 1.0 + x(2) );
-  else
-    f = theta       *       x(1).^tau   + x(2) * ( 1.0 + x(2) );
-  end
-endfunction
-function f = mckinnon2 ( x )
-
-  if ( length ( x ) ~= 2 )
-    error ( 'Error: function expects a two dimensional input\n' );
-  end
-
-  tau = 2.0;
-  theta = 6.0;
-  phi = 60.0;
-
-  if ( x(1) <= 0.0 )
-    f = theta * phi * abs ( x(1) ).^tau + x(2) * ( 1.0 + x(2) );
-  else
-    f = theta       *       x(1).^tau   + x(2) * ( 1.0 + x(2) );
-  end
-endfunction
-function f = mckinnon3 ( x )
-
-  if ( length ( x ) ~= 2 )
-    error ( 'Error: function expects a two dimensional input\n' );
-  end
-
-  tau = 3.0;
-  theta = 6.0;
-  phi = 400.0;
-
-  if ( x(1) <= 0.0 )
-    f = theta * phi * abs ( x(1) ).^tau + x(2) * ( 1.0 + x(2) );
-  else
-    f = theta       *       x(1).^tau   + x(2) * ( 1.0 + x(2) );
-  end
-endfunction
-
diff --git a/scilab/modules/optimization/demos/neldermead/neldermead.dem.sce b/scilab/modules/optimization/demos/neldermead/neldermead.dem.sce
deleted file mode 100644 (file)
index 14c83d1..0000000
+++ /dev/null
@@ -1,12 +0,0 @@
-//
-// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) ????-2008 - INRIA
-//
-// This file is distributed under the same license as the Scilab package.
-//
-
-mode(-1);
-
-path=get_absolute_file_path('datafit.dem.sce');
-exec(path+'neldermead_rosenbrock.sce');
-demo_datafit();
@@ -13,6 +13,9 @@
 //   numerical experiment presented in Box's paper.
 //
 
+mprintf("Illustrates Box''' algorithm on Box problem A.\n");
+mprintf("Defining Box Problem A function...\n");
+
 //
 // boxproblemA --
 //   Computes the Box problem A cost function and 
@@ -34,7 +37,7 @@
 //    constraints
 //  The inequality constraints are expected to be positive.
 //
-function f = boxproblemA ( x , index , data )
+function result = boxproblemA ( x , index , data )
   if (~isdef('index','local')) then
     index = 1
   end
@@ -56,7 +59,7 @@ function f = boxproblemA ( x , index , data )
   if ( index==1 | index==3 ) then
     f = (data.a2 * y1 + data.a3 * y2 + data.a4 * y3 + data.a5 * y4 ...
              + 7840 * data.a6 - 100000 * data.a0 ...
-             - 50800 * $b * data.a7 + data.k31 + data.k32 * x(2) + data.k33 * x(3) ...
+             - 50800 * b * data.a7 + data.k31 + data.k32 * x(2) + data.k33 * x(3) ...
              + data.k34 * x(4) + data.k35 * x(5)) * x(1) ...
              - 24345 + data.a1 * x6
     f = -f
@@ -82,77 +85,101 @@ function f = boxproblemA ( x , index , data )
   end
 endfunction
 
-boxparams = struct()
-boxparams.a0 = 9
-boxparams.a1 = 15
-boxparams.a2 = 50
-boxparams.a3 = 9.583
-boxparams.a4 = 20
-boxparams.a5 = 15
-boxparams.a6 = 6
-boxparams.a7 = 0.75
-boxparams.k1 = -145421.402
-boxparams.k2 = 2931.1506
-boxparams.k3 = -40.427932
-boxparams.k4 = 5106.192
-boxparams.k5 = 15711.36
-boxparams.k6 = -161622.577
-boxparams.k7 = 4176.15328
-boxparams.k8 = 2.8260078
-boxparams.k9 = 9200.476
-boxparams.k10 = 13160.295
-boxparams.k11 = -21686.9194
-boxparams.k12 = 123.56928
-boxparams.k13 = -21.1188894
-boxparams.k14 = 706.834
-boxparams.k15 = 2898.573
-boxparams.k16 = 28298.388
-boxparams.k17 = 60.81096
-boxparams.k18 = 31.242116
-boxparams.k19 = 329.574
-boxparams.k20 = -2882.082
-boxparams.k21 = 74095.3845
-boxparams.k22 = -306.262544
-boxparams.k23 = 16.243649
-boxparams.k24 = -3094.252
-boxparams.k25 = -5566.2628
-boxparams.k26 = -26237.0
-boxparams.k27 = 99.0
-boxparams.k28 = -0.42
-boxparams.k29 = 1300.0
-boxparams.k30 = 2100.0
-boxparams.k31 = 925548.252
-boxparams.k32 = -61968.8432
-boxparams.k33 = 23.3088196
-boxparams.k34 = -27097.648
-boxparams.k35 = -50843.766
-
+boxparams = struct();
+boxparams.a0 = 9;
+boxparams.a1 = 15;
+boxparams.a2 = 50;
+boxparams.a3 = 9.583;
+boxparams.a4 = 20;
+boxparams.a5 = 15;
+boxparams.a6 = 6;
+boxparams.a7 = 0.75;
+boxparams.k1 = -145421.402;
+boxparams.k2 = 2931.1506;
+boxparams.k3 = -40.427932;
+boxparams.k4 = 5106.192;
+boxparams.k5 = 15711.36;
+boxparams.k6 = -161622.577;
+boxparams.k7 = 4176.15328;
+boxparams.k8 = 2.8260078;
+boxparams.k9 = 9200.476;
+boxparams.k10 = 13160.295;
+boxparams.k11 = -21686.9194;
+boxparams.k12 = 123.56928;
+boxparams.k13 = -21.1188894;
+boxparams.k14 = 706.834;
+boxparams.k15 = 2898.573;
+boxparams.k16 = 28298.388;
+boxparams.k17 = 60.81096;
+boxparams.k18 = 31.242116;
+boxparams.k19 = 329.574;
+boxparams.k20 = -2882.082;
+boxparams.k21 = 74095.3845;
+boxparams.k22 = -306.262544;
+boxparams.k23 = 16.243649;
+boxparams.k24 = -3094.252;
+boxparams.k25 = -5566.2628;
+boxparams.k26 = -26237.0;
+boxparams.k27 = 99.0;
+boxparams.k28 = -0.42;
+boxparams.k29 = 1300.0;
+boxparams.k30 = 2100.0;
+boxparams.k31 = 925548.252;
+boxparams.k32 = -61968.8432;
+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)
 
+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 );
+xopt = [4.53743 2.4 60.0 9.3 7.0].';
+// Compute f(xopt) : should be -5280334.0
+fopt = boxproblemA ( xopt , data = boxparams );
 
-nm = nmplot_new ();
-nm = nmplot_configure(nm,"-numberofvariables",2);
-nm = nmplot_configure(nm,"-function",boxproblemA);
-nm = nmplot_configure(nm,"-costfargument",boxparams);
-nm = nmplot_configure(nm,"-x0",[1.0 1.0]');
-nm = nmplot_configure(nm,"-maxiter",200);
-nm = nmplot_configure(nm,"-maxfunevals",300);
-nm = nmplot_configure(nm,"-simplex0method","given");
-nm = nmplot_configure(nm,"-tolsimplexizerelative",1.e-3);
-nm = nmplot_configure(nm,"-coords0",coords0);
-nm = nmplot_configure(nm,"-simplex0length",1.0);
-nm = nmplot_configure(nm,"-method","box");
-//nm = nmplot_configure(nm,"-verbose",1);
-nm = nmplot_configure(nm,"-verbosetermination",1);
+nm = neldermead_new ();
+nm = neldermead_configure(nm,"-numberofvariables",5);
+nm = neldermead_configure(nm,"-function",boxproblemA);
+nm = neldermead_configure(nm,"-costfargument",boxparams);
+nm = neldermead_configure(nm,"-x0",x0);
+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,"-verbosetermination",1);
 // Configure like Box
 nm = neldermead_configure(nm,"-boundsmin",[0.0 1.2 20.0 9.0 6.0]);
 nm = neldermead_configure(nm,"-boundsmax",[5.0 2.5 60.0 9.3 7.0]);
 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 = nmplot_search(nm);
-nm = nmplot_destroy(nm);
+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);
 
 
index 657715d..31976f4 100644 (file)
@@ -32,7 +32,18 @@ mprintf("Searching for minimum...\n");
 nm = neldermead_search(nm);
 neldermead_display(nm);
 mprintf("Plot contour...\n");
-[nm , xdata , ydata , zdata ] = neldermead_contour ( nm , xmin = -2.0 , xmax = 2.0 , ymin = -2.0 , ymax = 2.0 , nx = 100 , ny = 100 );
+xmin = -2.0 ; xmax = 2.0 ; ymin = -2.0 ; ymax = 2.0 ; nx = 100 ; ny = 100;
+stepx = (xmax - xmin)/nx
+xdata = xmin:stepx:xmax;
+stepy = (ymax - ymin)/ny
+ydata = ymin:stepy:ymax;
+for ix = 1:length(xdata)
+  for iy = 1:length(ydata)
+    experiment = [xdata(ix) ydata(iy)];
+    [ nm , fiexp ] = neldermead_function ( nm , experiment );
+    zdata ( ix , iy ) = fiexp;
+  end
+end
 wnum = 100001;
 my_handle             = scf(wnum);
 contour ( xdata , ydata , zdata , [1 10 100 500 1000 2000] )
index 6d6a82c..299ab13 100644 (file)
@@ -7,6 +7,7 @@
 // are also available at
 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 
+mprintf("Illustrates the 1st counter example given by Han et al.\n");
 
 //
 // han1 --
@@ -68,8 +69,7 @@ nmplot_display(nm);
 mprintf("Plotting contour...\n");
 [nm , xdata , ydata , zdata ] = nmplot_contour ( nm , xmin = -0.2 , xmax = 1.2 , ymin = -2.0 , ymax = 2.0 , nx = 50 , ny = 50 );
 //[nm , xdata , ydata , zdata ] = nmplot_contour ( nm , xmin = -0.2 , xmax = 1.2 , ymin = -1.2 , ymax = 1.2 , nx = 50 , ny = 50 );
-wnum = 100001;
-f = scf(wnum);
+f = scf(100001);
 xset("fpf"," ")
 contour ( xdata , ydata , zdata , [-5 -4 -2 -1 0 1 1.5] )
 nmplot_simplexhistory ( nm );
index e7fa80d..85f79c4 100644 (file)
@@ -7,6 +7,7 @@
 // are also available at
 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 
+mprintf("Illustrates the 2nd counter example given by Han et al.\n");
 
 //
 // han2 --
@@ -31,11 +32,13 @@ endfunction
 
 
 coords0 = [
-0.0 0.0 1.0
-0.5 -0.5 0.0
+    0.    0.5  
+    0.   -0.5  
+    1.    0.   
 ]
 
 
+mprintf("Creating nmplot object...\n");
 nm = nmplot_new ();
 nm = nmplot_configure(nm,"-numberofvariables",2);
 nm = nmplot_configure(nm,"-function",han2);
@@ -57,16 +60,17 @@ nm = nmplot_configure(nm,"-simplexfn","han2-history-simplex.txt");
 //
 // Perform optimization
 //
+mprintf("Searching for minimum...\n");
 nm = nmplot_search(nm);
+nmplot_display(nm);
 //
 // Plot
 //
 [nm , xdata , ydata , zdata ] = nmplot_contour ( nm , xmin = -0.2 , xmax = 1.2 , ymin = -1.5 , ymax = 1.5 , nx = 50 , ny = 50 );
-f = scf();
+f = scf(100001);
 xset("fpf"," ")
-contour ( xdata , ydata , zdata , 40 )
+contour ( xdata , ydata , zdata , [0.1 0.2 0.5 1.0 1.5 1.9] )
 nmplot_simplexhistory ( nm );
-xs2png(0,"han2-history-simplex.png");
 deletefile("han2-history-simplex.txt");
 nm = nmplot_destroy(nm);
 
index fc1af7b..901afad 100644 (file)
@@ -8,7 +8,7 @@
 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 
 
-//dirname = get_absolute_file_path('nmplot_mckinnon.sce');
+mprintf("Defining McKinnon function...\n");
 //% MCKINNON computes the McKinnon function.
 //
 //  Discussion:
@@ -91,18 +91,22 @@ function f = mckinnon3 ( x )
 endfunction
 
 
-lambda1 = (1.0 + sqrt(33.0))/8.0
-lambda2 = (1.0 - sqrt(33.0))/8.0
+lambda1 = (1.0 + sqrt(33.0))/8.0;
+lambda2 = (1.0 - sqrt(33.0))/8.0;
 coords0 = [
-1.0 0.0 lambda1
-1.0 0.0 lambda2
-]
+1.0 1.0
+0.0 0.0
+lambda1 lambda2
+];
 
 
+x0 = [1.0 1.0]';
+mprintf("x0=%s\n",strcat(string(x0)," "));
+mprintf("Creating object...\n");
 nm = nmplot_new ();
 nm = nmplot_configure(nm,"-numberofvariables",2);
 nm = nmplot_configure(nm,"-function",mckinnon3);
-nm = nmplot_configure(nm,"-x0",[1.0 1.0]');
+nm = nmplot_configure(nm,"-x0",x0);
 nm = nmplot_configure(nm,"-maxiter",200);
 nm = nmplot_configure(nm,"-maxfunevals",300);
 nm = nmplot_configure(nm,"-tolfunrelative",10*%eps);
@@ -128,18 +132,19 @@ nmplot_display(nm);
 //
 // Plot
 //
+mprintf("Plot contour...\n");
 [nm , xdata , ydata , zdata ] = nmplot_contour ( nm , xmin = -0.2 , xmax = 1.2 , ymin = -2.0 , ymax = 2.0 , nx = 50 , ny = 50 );
-f = scf();
+f = scf(100001);
 xset("fpf"," ")
-contour ( xdata , ydata , zdata , 40 )
+contour ( xdata , ydata , zdata , [-0.2 0.0 1.0 2.0 5.0 10.0 20.0] )
 nmplot_simplexhistory ( nm );
-f = scf();
+f = scf(100002);
 nmplot_historyplot ( nm , "mckinnon.history.fbar.txt" , ...
   mytitle = "Function Value Average" , myxlabel = "Iterations" );
-f = scf();
+f = scf(100003);
 nmplot_historyplot ( nm , "mckinnon.history.fopt.txt" , ...
   mytitle = "Minimum Function Value" , myxlabel = "Iterations" );
-f = scf();
+f = scf(100004);
 nmplot_historyplot ( nm , "mckinnon.history.sigma.txt" , ...
   mytitle = "Maximum Oriented length" , myxlabel = "Iterations" );
 deletefile("mckinnon.history.simplex.txt");
index 9238927..c99c74f 100644 (file)
@@ -8,7 +8,7 @@
 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 
 
-//dirname = get_absolute_file_path('nmplot_mckinnon.sce');
+mprintf("Defining McKinnon function...\n");
 
 //% MCKINNON computes the McKinnon function.
 //
@@ -91,18 +91,22 @@ function f = mckinnon3 ( x )
   end
 endfunction
 
-lambda1 = (1.0 + sqrt(33.0))/8.0
-lambda2 = (1.0 - sqrt(33.0))/8.0
+lambda1 = (1.0 + sqrt(33.0))/8.0;
+lambda2 = (1.0 - sqrt(33.0))/8.0;
 coords0 = [
-1.0 0.0 lambda1
-1.0 0.0 lambda2
-]
+1.0  1.0
+0.0  0.0
+lambda1 lambda2
+];
 
 
+x0 = [1.0 1.0]';
+mprintf("x0=%s\n",strcat(string(x0)," "));
+mprintf("Creating object...\n");
 nm = nmplot_new ();
 nm = nmplot_configure(nm,"-numberofvariables",2);
 nm = nmplot_configure(nm,"-function",mckinnon3);
-nm = nmplot_configure(nm,"-x0",[1.0 1.0]');
+nm = nmplot_configure(nm,"-x0",x0);
 nm = nmplot_configure(nm,"-maxiter",200);
 nm = nmplot_configure(nm,"-maxfunevals",300);
 nm = nmplot_configure(nm,"-tolsimplexizerelative",1.e-6);
@@ -130,10 +134,11 @@ nmplot_display(nm);
 //
 // Plot
 //
+mprintf("Plot contour...\n");
 [nm , xdata , ydata , zdata ] = nmplot_contour ( nm , xmin = -0.2 , xmax = 2.0 , ymin = -2.0 , ymax = 2.0 , nx = 100 , ny = 100 );
 f = scf();
 xset("fpf"," ")
-contour ( xdata , ydata , zdata , 40 )
+contour ( xdata , ydata , zdata , [-0.2 0.0 1.0 2.0 5.0 10.0 20.0] )
 nmplot_simplexhistory ( nm );
 f = scf();
 nmplot_historyplot ( nm , "mckinnon.history.restart.fbar.txt" , ...
index c7893bf..71bb2f3 100644 (file)
@@ -8,6 +8,7 @@
 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 
 
+mprintf("Illustrates that the fixed-shape Spendley et al. algorithm performs well on a quadratic test case.\n");
 mprintf("Defining quadratic function...\n");
 function y = quadratic (x)
   y = x(1)^2 + x(2)^2 - x(1) * x(2);
@@ -20,12 +21,11 @@ nm = nmplot_configure(nm,"-function",quadratic);
 nm = nmplot_configure(nm,"-x0",[2.0 2.0]');
 nm = nmplot_configure(nm,"-maxiter",100);
 nm = nmplot_configure(nm,"-maxfunevals",300);
-//nm = nmplot_configure(nm,"-tolfunrelative",10*%eps);
-nm = nmplot_configure(nm,"-tolxrelative",1.e-8);
+nm = nmplot_configure(nm,"-tolxmethod","disabled");
+nm = nmplot_configure(nm,"-tolsimplexizerelative",1.e-8);
 nm = nmplot_configure(nm,"-simplex0method","spendley");
-//nm = nmplot_configure(nm,"-simplex0length",1.0);
 nm = nmplot_configure(nm,"-method","fixed");
-//nm = nmplot_configure(nm,"-verbose",1);
+nm = nmplot_configure(nm,"-verbose",1);
 nm = nmplot_configure(nm,"-verbosetermination",0);
 //
 // Setup output files
@@ -40,12 +40,7 @@ nm = nmplot_configure(nm,"-sigmafn","rosenbrock.fixed.history.sigma.txt");
 mprintf("Searching for minimum...\n");
 nm = nmplot_search(nm);
 nmplot_display(nm);
-// Plot the contours of the cost function and the simplex history
-mprintf("Plotting contour...\n");
-[nm , xdata , ydata , zdata ] = nmplot_contour ( nm , xmin = -2.0 , xmax = 4.0 , ymin = -2.0 , ymax = 4.0 , nx = 100 , ny = 100 );
-f = scf();
-contour ( xdata , ydata , zdata , [0.1 1.0 2.0 5.0 10.0 15.0 20.0] )
-nmplot_simplexhistory ( nm );
+// Plot various histories
 mprintf("Plotting history of fbar...\n");
 f = scf();
 nmplot_historyplot ( nm , "rosenbrock.fixed.history.fbar.txt" , ...
@@ -53,11 +48,32 @@ nmplot_historyplot ( nm , "rosenbrock.fixed.history.fbar.txt" , ...
 mprintf("Plotting history of fopt...\n");
 f = scf();
 nmplot_historyplot ( nm , "rosenbrock.fixed.history.fopt.txt" , ...
-  mytitle = "Minimum Function Value" , myxlabel = "Iterations" );
+  mytitle = "Logarithm Minimum Function Value" , myxlabel = "Iterations" );
+f.children.log_flags = "nln";
+newticks = tlist(["ticks","locations","labels"]);
+newticks.labels = ["1.e-20" "1.e-10" "1.e-1"];
+newticks.locations = [1.e-20 1.e-10 1.e-1];
+f.children.y_ticks = newticks;
+f.children.children(1).children.mark_mode = "on";
+f.children.children(1).children.mark_style = 9;
 mprintf("Plotting history of sigma...\n");
 f = scf();
 nmplot_historyplot ( nm , "rosenbrock.fixed.history.sigma.txt" , ...
-  mytitle = "Maximum Oriented length" , myxlabel = "Iterations" );
+  mytitle = "Logarithm Maximum Oriented length" , myxlabel = "Iterations" );
+f.children.log_flags = "nln";
+f.children.y_ticks = newticks;
+f.children.children(1).children.mark_mode = "on";
+f.children.children(1).children.mark_style = 9;
+// Plot the contours of the cost function and the simplex history
+mprintf("Plotting contour...\n");
+nm = nmplot_configure(nm,"-verbose",0);
+[nm , xdata , ydata , zdata ] = nmplot_contour ( nm , xmin = -2.0 , xmax = 4.0 , ymin = -2.0 , ymax = 4.0 , nx = 100 , ny = 100 );
+f = scf();
+drawlater();
+contour ( xdata , ydata , zdata , [0.1 1.0 2.0 5.0 10.0 15.0 20.0] )
+nmplot_simplexhistory ( nm );
+drawnow();
+// Clean-up
 deletefile("rosenbrock.fixed.history.simplex.txt");
 deletefile("rosenbrock.fixed.history.fbar.txt");
 deletefile("rosenbrock.fixed.history.fopt.txt");
diff --git a/scilab/modules/optimization/demos/neldermead/nmplot_quadratic.fixed2.sce b/scilab/modules/optimization/demos/neldermead/nmplot_quadratic.fixed2.sce
new file mode 100644 (file)
index 0000000..a76bc8b
--- /dev/null
@@ -0,0 +1,83 @@
+// 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
+
+
+mprintf("Illustrates that the fixed-shape Spendley et al. algorithm performs badly on a badly quadratic test case.\n");
+mprintf("Defining quadratic function...\n");
+function y = quadratic (x)
+  y = 100 * x(1)^2 + x(2)^2;
+endfunction
+
+mprintf("Creating nmplot object...\n");
+nm = nmplot_new ();
+nm = nmplot_configure(nm,"-numberofvariables",2);
+nm = nmplot_configure(nm,"-function",quadratic);
+nm = nmplot_configure(nm,"-x0",[10.0 10.0]');
+nm = nmplot_configure(nm,"-maxiter",400);
+nm = nmplot_configure(nm,"-maxfunevals",400);
+nm = nmplot_configure(nm,"-tolxmethod","disabled");
+nm = nmplot_configure(nm,"-tolsimplexizerelative",1.e-8);
+nm = nmplot_configure(nm,"-simplex0method","spendley");
+nm = nmplot_configure(nm,"-method","fixed");
+nm = nmplot_configure(nm,"-verbose",1);
+nm = nmplot_configure(nm,"-verbosetermination",0);
+//
+// Setup output files
+//
+nm = nmplot_configure(nm,"-simplexfn","rosenbrock.fixed.history.simplex.txt");
+nm = nmplot_configure(nm,"-fbarfn","rosenbrock.fixed.history.fbar.txt");
+nm = nmplot_configure(nm,"-foptfn","rosenbrock.fixed.history.fopt.txt");
+nm = nmplot_configure(nm,"-sigmafn","rosenbrock.fixed.history.sigma.txt");
+//
+// Perform optimization
+//
+mprintf("Searching for minimum...\n");
+nm = nmplot_search(nm);
+nmplot_display(nm);
+// Plot various histories
+mprintf("Plotting history of fbar...\n");
+f = scf();
+nmplot_historyplot ( nm , "rosenbrock.fixed.history.fbar.txt" , ...
+  mytitle = "Function Value Average" , myxlabel = "Iterations" );
+mprintf("Plotting history of fopt...\n");
+f = scf();
+nmplot_historyplot ( nm , "rosenbrock.fixed.history.fopt.txt" , ...
+  mytitle = "Logarithm Minimum Function Value" , myxlabel = "Iterations" );
+f.children.log_flags = "nln";
+newticks = tlist(["ticks","locations","labels"]);
+newticks.labels = ["1.e-20" "1.e-10" "1.e-1"];
+newticks.locations = [1.e-20 1.e-10 1.e-1];
+f.children.y_ticks = newticks;
+f.children.children(1).children.mark_mode = "on";
+f.children.children(1).children.mark_style = 9;
+mprintf("Plotting history of sigma...\n");
+f = scf();
+nmplot_historyplot ( nm , "rosenbrock.fixed.history.sigma.txt" , ...
+  mytitle = "Logarithm Maximum Oriented length" , myxlabel = "Iterations" );
+f.children.log_flags = "nln";
+f.children.y_ticks = newticks;
+f.children.children(1).children.mark_mode = "on";
+f.children.children(1).children.mark_style = 9;
+// Plot the contours of the cost function and the simplex history
+mprintf("Plotting contour...\n");
+nm = nmplot_configure(nm,"-verbose",0);
+[nm , xdata , ydata , zdata ] = nmplot_contour ( nm , xmin = -5.0 , xmax = 12.0 , ymin = -2.0 , ymax = 12.0 , nx = 100 , ny = 100 );
+f = scf();
+drawlater();
+contour ( xdata , ydata , zdata , [10.0 50 100 1000 2000 5000 10000 20000] )
+nmplot_simplexhistory ( nm );
+drawnow();
+// Clean-up
+deletefile("rosenbrock.fixed.history.simplex.txt");
+deletefile("rosenbrock.fixed.history.fbar.txt");
+deletefile("rosenbrock.fixed.history.fopt.txt");
+deletefile("rosenbrock.fixed.history.sigma.txt");
+nm = nmplot_destroy(nm);
+
+
index 6d55904..407ae02 100644 (file)
@@ -8,6 +8,7 @@
 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 
 
+mprintf("Illustrates that the fixed-shape Spendley et al. algorithm does NOT perform well on Rosenbrock test case.\n");
 mprintf("Defining Rosenbrock function...\n");
 function y = rosenbrock (x)
   y = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;
@@ -39,11 +40,12 @@ nm = nmplot_configure(nm,"-sigmafn","rosenbrock.fixed.history.sigma.txt");
 //
 mprintf("Searching for minimum...\n");
 nm = nmplot_search(nm);
+nmplot_display(nm);
 // Plot the contours of the cost function and the simplex history
 mprintf("Plotting contour...\n");
 [nm , xdata , ydata , zdata ] = nmplot_contour ( nm , xmin = -2.0 , xmax = 2.0 , ymin = -2.0 , ymax = 2.0 , nx = 100 , ny = 100 );
 f = scf();
-contour ( xdata , ydata , zdata , 20 )
+contour ( xdata , ydata , zdata , [1 10 100 500 1000 2000] )
 nmplot_simplexhistory ( nm );
 mprintf("Plotting history of fbar...\n");
 f = scf();
index 040d5b8..88a0430 100644 (file)
@@ -8,6 +8,7 @@
 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 
 
+mprintf("Illustrates the Nelder-Mead algorithm on Rosenbrock test case.\n");
 mprintf("Defining Rosenbrock function...\n");
 function y = rosenbrock (x)
   y = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;
index 9bb824c..98abe49 100644 (file)
@@ -10,12 +10,13 @@ subdemolist = [
 "Non linear data fitting"         "datafit/datafit.dem.sce"       
 "fminsearch"                      "neldermead/fminsearch.sce"       
 "neldermead/Rosenbrock Variable"  "neldermead/neldermead_rosenbrock.sce"       
-"nmplot/BoxA"               "neldermead/nmplot_boxproblemA.sce"       
+"neldermead/Box A"                "neldermead/neldermead_boxproblemA.sce"       
 "nmplot/Han #1"             "neldermead/nmplot_han1.sce"       
 "nmplot/Han #2"             "neldermead/nmplot_han2.sce"       
 "nmplot/McKinnon #1"        "neldermead/nmplot_mckinnon.sce"       
 "nmplot/McKinnon #2"        "neldermead/nmplot_mckinnon2.sce"       
-"nmplot/Quadratic Fixed"   "neldermead/nmplot_quadratic.fixed.sce"       
+"nmplot/Quadratic Fixed #1"   "neldermead/nmplot_quadratic.fixed.sce"       
+"nmplot/Quadratic Fixed #2"   "neldermead/nmplot_quadratic.fixed2.sce"       
 "nmplot/Rosenbrock Fixed"   "neldermead/nmplot_rosenbrock.fixed.sce"       
 "nmplot/Rosenbrock"         "neldermead/nmplot_rosenbrock.sce"       
 ];
index 3640ef4..5ae7e0f 100644 (file)
@@ -32,7 +32,6 @@ this = neldermead_display ( this )
 value = neldermead_get ( this , key )
 this = neldermead_search ( this )
 this = neldermead_restart ( this )
-[ this , xdata , ydata , zdata ] = neldermead_contour ( this , xmin , xmax , ymin , ymax , nx , ny )
 </synopsis>
   </refsynopsisdiv>
 
@@ -1656,6 +1655,71 @@ alpha = alpha0 * sigma0 / nsg
   </refsection>
 
   <refsection>
+    <title>Spendley et al. implementation notes</title>
+
+    <para>The original paper may be implemented with several variations, which
+    might lead to different results. This section defines what algorithmic
+    choices have been used.</para>
+
+    <para>The paper states the following rules.</para>
+
+    <itemizedlist>
+      <listitem>
+        <para>"Rule 1. Ascertain the lowest reading y, of yi ... yk+1 Complete
+        a new simplex Sp by excluding the point Vp corresponding to y, and
+        replacing it by V* defined as above."</para>
+      </listitem>
+
+      <listitem>
+        <para>"Rule 2. If a result has occurred in (k + 1) successive
+        simplexes, and is not then eliminated by application of Rule 1, do not
+        move in the direction indicated by Rule 1, or at all, but discard the
+        result and replace it by a new observation at the same point."</para>
+      </listitem>
+
+      <listitem>
+        <para>"Rule 3. If y is the lowest reading in So , and if the next
+        observation made, y* , is the lowest reading in the new simplex S , do
+        not apply Rule 1 and return to So from Sp . Move out of S, by
+        rejecting the second lowest reading (which is also the second lowest
+        reading in So)."</para>
+      </listitem>
+    </itemizedlist>
+
+    <para>We implement the following "rules" of the Spendley et al.
+    method.</para>
+
+    <itemizedlist>
+      <listitem>
+        <para>Rule 1 is strictly applied, but the reflection is done by
+        reflection the high point, since we minimize a function instead of
+        maximizing it, like Spendley.</para>
+      </listitem>
+
+      <listitem>
+        <para>Rule 2 is NOT implemented, as we expect that the function
+        evaluation is not subject to errors.</para>
+      </listitem>
+
+      <listitem>
+        <para>Rule 3 is applied, ie reflection with respect to next to high
+        point. </para>
+      </listitem>
+    </itemizedlist>
+
+    <para>The original paper does not mention any shrink step. When the
+    original algorithm cannot improve the function value with reflection
+    steps, the basic algorithm stops. In order to make the current
+    implementation of practical value, a shrink step is included, with
+    shrinkage factor sigma. This perfectly fits into to the spirit of the
+    original paper. Notice that the shrink step make the rule #3 (reflection
+    with respect to next-to-worst vertex) unnecessary. Indeed, the minimum
+    required steps are the reflection and shrinkage. Never the less, the rule
+    #3 has been kept in order to make the algorithm as close as it can be to
+    the original.</para>
+  </refsection>
+
+  <refsection>
     <title>Example #1</title>
 
     <para>In the following example, we solve the Rosenbrock test case. We
@@ -1702,7 +1766,18 @@ nm = neldermead_configure(nm,"-verbosetermination",0);
 mprintf("Searching for minimum...\n");
 nm = neldermead_search(nm);
 mprintf("Plot contour...\n");
-[nm , xdata , ydata , zdata ] = neldermead_contour ( nm , xmin = -2.0 , xmax = 2.0 , ymin = -2.0 , ymax = 2.0 , nx = 100 , ny = 100 );
+xmin = -2.0 ; xmax = 2.0 ; ymin = -2.0 ; ymax = 2.0 ; nx = 100 ; ny = 100;
+stepx = (xmax - xmin)/nx
+xdata = xmin:stepx:xmax;
+stepy = (ymax - ymin)/ny
+ydata = ymin:stepy:ymax;
+for ix = 1:length(xdata)
+  for iy = 1:length(ydata)
+    experiment = [xdata(ix) ydata(iy)];
+    [ nm , fiexp ] = neldermead_function ( nm , experiment );
+    zdata ( ix , iy ) = fiexp;
+  end
+end
 wnum = 100001;
 my_handle             = scf(wnum);
 contour ( xdata , ydata , zdata , [1 10 100 500 1000 2000] )
@@ -1758,28 +1833,20 @@ Sort
   <refsection>
     <title>TODO</title>
 
-    <variablelist>
-      <varlistentry>
-        <term>Box-Guin algorithm</term>
-
-        <listitem>
-          <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 is an update of Box's algorithm.</para>
-        </listitem>
-      </varlistentry>
+    <itemizedlist>
+      <listitem>
+        <para>add the Box-Guin algoritm as a 4th method.</para>
 
-      <varlistentry>
-        <term>add the optimization of the Rosenbrock test case, with the
-        interactive plot thanks to the -outputcommand option</term>
-      </varlistentry>
+        <para>This algorithm solves an constrained optimization problem with a
+        variable sized simplex made of an arbitrary k number of vertices. This
+        is an update of Box's algorithm.</para>
+      </listitem>
 
-      <varlistentry>
-        <term>add a sample output with the verbose option.</term>
-      </varlistentry>
-    </variablelist>
+      <listitem>
+        <para>add the optimization of the Rosenbrock test case, with the
+        interactive plot thanks to the -outputcommand option</para>
+      </listitem>
+    </itemizedlist>
   </refsection>
 
   <refsection>
@@ -1815,7 +1882,9 @@ Sort
   <refsection>
     <title>Authors</title>
 
-    <para>Michael Baudin, 2008-2009</para>
+    <para>Michael Baudin - INRIA - 2008-2009</para>
+
+    <para>Michael Baudin - Digiteo - 2009</para>
   </refsection>
 
   <refsection>
index afc0bb2..0ba83cd 100644 (file)
@@ -46,7 +46,7 @@ this = nmplot_simplexhistory ( this , colorforeground , markforeground , marksty
 
     <para>The goal of this class is to provide a neldermead component with
     plotting features. It enables to make fast plots of the algorithm progress
-    through the iterations. </para>
+    through the iterations.</para>
 
     <para>It is a specialized neldermead class, with a specific output
     command. This output function allows to store the history of several datas
@@ -1328,17 +1328,99 @@ this = nmplot_simplexhistory ( this , colorforeground , markforeground , marksty
   </refsection>
 
   <refsection>
-    <title>TODO</title>
+    <title>Example</title>
+
+    <para>In the following example, we use the fixed shape Spendley et al.
+    simplex algorithm and find the minimum of a quadratic function. We begin
+    by defining a quadratic function associated with 2 input variables. We
+    then define an nmplot object and configure the object so that the "fixed"
+    shape simplex algorithm is used with the regular initial simplex
+    associated with the "spendley" key. Four files are configured, which will
+    contain the history of the simplex, the history of fbar, fopt and sigma
+    through the iterations. The search is performed by the nmplot_search
+    function, which writes the 4 data files during the iterations. The
+    nmplot_contour function is called in order to compute the arrays xdata,
+    ydata and zdata which are required as input to the contour function. The
+    nmplot_simplexhistory then uses the history of the simplex, as stored in
+    the rosenbrock.fixed.history.simplex.txt data file, and plot the various
+    simplices on the contour plot. The nmplot_historyplot is used with the
+    files rosenbrock.fixed.history.fbar.txt, rosenbrock.fixed.history.fopt.txt
+    and rosenbrock.fixed.history.sigma.txt, which produces 3 plots of the
+    history of the optimization algorithm through the iterations.</para>
+
+    <programlisting role="example"> 
+mprintf("Defining quadratic function...\n");
+function y = quadratic (x)
+  y = x(1)^2 + x(2)^2 - x(1) * x(2);
+endfunction
+
+mprintf("Creating nmplot object...\n");
+nm = nmplot_new ();
+nm = nmplot_configure(nm,"-numberofvariables",2);
+nm = nmplot_configure(nm,"-function",quadratic);
+nm = nmplot_configure(nm,"-x0",[2.0 2.0]');
+nm = nmplot_configure(nm,"-maxiter",100);
+nm = nmplot_configure(nm,"-maxfunevals",300);
+nm = nmplot_configure(nm,"-tolxrelative",1.e-8);
+nm = nmplot_configure(nm,"-simplex0method","spendley");
+nm = nmplot_configure(nm,"-method","fixed");
+//
+// Setup output files
+//
+nm = nmplot_configure(nm,"-simplexfn","rosenbrock.fixed.history.simplex.txt");
+nm = nmplot_configure(nm,"-fbarfn","rosenbrock.fixed.history.fbar.txt");
+nm = nmplot_configure(nm,"-foptfn","rosenbrock.fixed.history.fopt.txt");
+nm = nmplot_configure(nm,"-sigmafn","rosenbrock.fixed.history.sigma.txt");
+//
+// Perform optimization
+//
+mprintf("Searching for minimum...\n");
+nm = nmplot_search(nm);
+nmplot_display(nm);
+// Plot the contours of the cost function and the simplex history
+mprintf("Plotting contour...\n");
+[nm , xdata , ydata , zdata ] = nmplot_contour ( nm , xmin = -2.0 , xmax = 4.0 , ymin = -2.0 , ymax = 4.0 , nx = 100 , ny = 100 );
+f = scf();
+contour ( xdata , ydata , zdata , [0.1 1.0 2.0 5.0 10.0 15.0 20.0] )
+nmplot_simplexhistory ( nm );
+mprintf("Plotting history of fbar...\n");
+f = scf();
+nmplot_historyplot ( nm , "rosenbrock.fixed.history.fbar.txt" , ...
+  mytitle = "Function Value Average" , myxlabel = "Iterations" );
+mprintf("Plotting history of fopt...\n");
+f = scf();
+nmplot_historyplot ( nm , "rosenbrock.fixed.history.fopt.txt" , ...
+  mytitle = "Minimum Function Value" , myxlabel = "Iterations" );
+mprintf("Plotting history of sigma...\n");
+f = scf();
+nmplot_historyplot ( nm , "rosenbrock.fixed.history.sigma.txt" , ...
+  mytitle = "Maximum Oriented length" , myxlabel = "Iterations" );
+deletefile("rosenbrock.fixed.history.simplex.txt");
+deletefile("rosenbrock.fixed.history.fbar.txt");
+deletefile("rosenbrock.fixed.history.fopt.txt");
+deletefile("rosenbrock.fixed.history.sigma.txt");
+nm = nmplot_destroy(nm);
+
+
+ </programlisting>
   </refsection>
 
   <refsection>
-    <title>Bibliography</title>
+    <title>TODO</title>
+
+    <itemizedlist>
+      <listitem>
+        <para>add an example</para>
+      </listitem>
+    </itemizedlist>
   </refsection>
 
   <refsection>
     <title>Authors</title>
 
-    <para>Michael Baudin, 2008-2009</para>
+    <para>Michael Baudin - INRIA - 2008-2009</para>
+
+    <para>Michael Baudin - Digiteo - 2009</para>
   </refsection>
 
   <refsection>
index 4fce9e5..e7df8c2 100644 (file)
@@ -190,7 +190,7 @@ cen = optimsimplex_xbar ( this , iexcl )</synopsis>
     <para>In the following functions, simplices and vertices are, depending on
     the functions either input or output arguments. The following general
     principle have been used to manage the storing of the coordinates of the
-    points. </para>
+    points.</para>
 
     <itemizedlist>
       <listitem>
@@ -1720,8 +1720,9 @@ function [ y , data ] = myfunction ( x , data )
 
     <programlisting role="example"> 
 coords = [
-0.0 1.0 0.0
-0.0 0.0 1.0
+    0.    0.  
+    1.    0.  
+    0.    1.  
 ];
 s1 = optimsimplex_new ( coords );
 computed = optimsimplex_getallx ( s1 );
@@ -1737,22 +1738,26 @@ s1 = optimsimplex_destroy(s1);
     <para>In the following example, one creates a simplex with in the 2D
     domain [-5 5]^2, with [-1.2 1.0] as the first vertex. One uses the
     randomized bounds method to generate a simplex with 5 vertices. The
-    function takes an additionnal argument myobj, which is counts the number
+    function takes an additionnal argument mystuff, which is counts the number
     of times the function is called. After the creation of the simplex, the
-    value of mydude.nb is 5, which is the expected result because there is one
-    function call by vertex.</para>
+    value of mystuff.nb is 5, which is the expected result because there is
+    one function call by vertex.</para>
 
     <programlisting role="example"> 
-function [ y , myobj ] = mycostf ( x , myobj )
+function y = rosenbrock (x)
+  y = 100*(x(2)-x(1)^2)^2+(1-x(1))^2;
+endfunction
+function [ y , mystuff ] = mycostf ( x , mystuff )
   y = rosenbrock(x);
-  myobj.nb = myobj.nb + 1
+  mystuff.nb = mystuff.nb + 1
 endfunction
 
-mydude = tlist(["T_MYSTUFF","nb"]);
-mydude.nb = 0;
+mystuff = tlist(["T_MYSTUFF","nb"]);
+mystuff.nb = 0;
 s1 = optimsimplex_new ();
-[ s1 , mydude ] = optimsimplex_randbounds ( s1 , x0 = [-1.2 1.0], fun = mycostf, ...
-  boundsmin = [-5.0 -5.0] , boundsmax = [5.0 5.0], nbve=5 , data = mydude );
+[ s1 , mystuff ] = optimsimplex_randbounds ( s1 , x0 = [-1.2 1.0], fun = mycostf, ...
+  boundsmin = [-5.0 -5.0] , boundsmax = [5.0 5.0], nbve=5 , data = mystuff );
+mprintf("Function evaluations: %d\n",mystuff.nb)
 s1 = optimsimplex_destroy ( s1 );
  </programlisting>
   </refsection>
@@ -1801,7 +1806,9 @@ s1 = optimsimplex_destroy ( s1 );
   <refsection>
     <title>Authors</title>
 
-    <para>Michael Baudin, 2008-2009</para>
+    <para>Michael Baudin - INRIA - 2008-2009</para>
+
+    <para>Michael Baudin - Digiteo - 2009</para>
   </refsection>
 
   <refsection>
index 0349b8b..ae01ffe 100644 (file)
@@ -75,8 +75,6 @@ function value = neldermead_cget (this,key)
     value = this.nbineqloops;
   case "-ineqscaling" then
     value = this.ineqscaling;
-  case "-vectorizefunction" then
-    value = this.vectorizefunction;
   else
     // Delegate to the optimization object
     value = optimbase_cget ( this.optbase , key );
index 5fde63b..f84c1fd 100644 (file)
@@ -119,17 +119,6 @@ function this = neldermead_configure (this,key,value)
     this.nbineqloops = value;
   case "-ineqscaling" then
     this.ineqscaling = value;
-  case "-vectorizefunction" then
-    this.vectorizefunction = value;
-    select value
-    case 0 then
-      this.vectorizefunction = 0;
-    case 1 then
-      this.vectorizefunction = 1;
-    else
-      errmsg = msprintf(gettext("%s: Unknown value %s for -vectorizefunction option"),"neldermead_configure", value);
-      error(errmsg);
-    end
   else
     // Delegate to the optimization object
     this.optbase = optimbase_configure ( this.optbase , key , value );
diff --git a/scilab/modules/optimization/macros/neldermead/neldermead_contour.sci b/scilab/modules/optimization/macros/neldermead/neldermead_contour.sci
deleted file mode 100644 (file)
index be1f795..0000000
+++ /dev/null
@@ -1,59 +0,0 @@
-// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 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
-
-//
-// neldermead_contour --
-//   Computes data necessary to create a contour plot of the cost function.
-//   The cost function must be a function with two parameters.
-// TODO : move to optimbase_compcontour
-// Arguments
-//   xmin , xmax , ymin , ymax : the bounds for the contour plot
-//   nx , ny : the number of points in the directions x, y
-//   xdata , ydata , zdata : vectors of data, as required by the contour command
-//
-function [ this , xdata , ydata , zdata ] = neldermead_contour ( this , xmin , xmax , ymin , ymax , nx , ny )
-  // Check that there are only 2 parameters
-  n = neldermead_cget ( this , "-numberofvariables" );
-  if n <> 2 then
-    errmsg = msprintf(gettext("%s: Unexpected number of variables %d. Cannot draw contour plot for functions which do not have two parameters."),"nmplot_contour",n)
-    error(errmsg)
-  end
-  stepx = (xmax - xmin)/nx
-  xdata = xmin:stepx:xmax;
-  stepy = (ymax - ymin)/ny
-  ydata = ymin:stepy:ymax;
-  // 1. Compute the matrix of experiments
-  nexp = length(xdata) * length(ydata)
-  iexp = 1;
-  for ix = 1:length(xdata)
-    for iy = 1:length(ydata)
-      x(iexp,1:2) = [xdata(ix) ydata(iy)];
-      iexp = iexp + 1;
-    end
-  end
-  // 2. Perform the experiments
-  // Vectorize the call to the cost function if possible
-  vectorizefunction = neldermead_cget ( this , "-vectorizefunction" );
-  if vectorizefunction == 0 then
-    for iexp = 1:nexp
-      [ this , f(iexp) ] = neldermead_function ( this , x(iexp,1:2) );
-    end
-  else
-    [ this , f ] = neldermead_function ( this , x );
-  end
-  // 3. Store the experiments results in a matrix, suitable for input to the contour function
-  iexp = 1;
-  for ix = 1:length(xdata)
-    for iy = 1:length(ydata)
-      zdata ( ix , iy ) = f(iexp);
-      iexp = iexp + 1;
-    end
-  end
-endfunction
-
index db518da..a823ca6 100644 (file)
 //   the requirements of simplex methods.
 //
 function [ f , this ] = neldermead_costf ( x , this )
-  [this.optbase , f] = optimbase_function ( this.optbase , x );
+  nbnlc = optimbase_cget ( this.optbase , "-nbineqconst" )
+  if ( nbnlc == 0 ) then
+    [this.optbase , f] = optimbase_function ( this.optbase , x );
+  else
+    // Set the index, so that, if an additionnal cost function argument is provided,
+    // it can be appended at the end.
+    index = 1;
+    [this.optbase , f] = optimbase_function ( this.optbase , x , index );
+  end
 endfunction
 
index 66beb65..c1bde27 100644 (file)
@@ -20,6 +20,11 @@ function this = neldermead_display ( this )
   for k =1:size(str,1)
     mprintf("%s\n",str(k));
   end
+  mprintf("Simplex Optimum : \n");
+  str = optimsimplex_tostring ( this.simplexopt );
+  for k =1:size(str,1)
+    mprintf("%s\n",str(k));
+  end
   mprintf("Simplex0 Method : %s\n", string(this.simplex0method));
   mprintf("Simplex0 Length : %s\n", string(this.simplex0length));
   mprintf("Termination Method on simplex size : %s\n", string(this.tolsimplexizemethod));
index dd933aa..023c417 100644 (file)
@@ -28,8 +28,7 @@ function newobj = neldermead_new ()
     "kelleynormalizationflag","kelleystagnationalpha0", ...
     "kelleyalpha","restartnb","restartflag","restartdetection" , ...
     "startupflag" , ...
-    "boxnbpoints" , "boxnbpointseff" , "nbineqloops" , "ineqscaling" , ...
-    "vectorizefunction" ]);
+    "boxnbpoints" , "boxnbpointseff" , "nbineqloops" , "ineqscaling" ]);
   newobj.optbase = optimbase_new();
   // Possible values "variable", "fixed".
   newobj.method = "variable";
@@ -114,7 +113,5 @@ function newobj = neldermead_new ()
   // The scaling coefficient in nonlinear inequality constraints
   // in Box method, in (0,1) range
   newobj.ineqscaling = 0.5
-  // Set to 1 if the cost function is vectorized
-  newobj.vectorizefunction = 0
 endfunction
 
index 63b39f3..a9d55ce 100644 (file)
@@ -247,6 +247,29 @@ endfunction
 //
 // neldermead_fixed --
 //   The simplex algorithm with fixed size simplex.
+// Implementation note:
+//   We implement the following "rules" of the Spendley et al.
+//   method.
+//   Rule 1 is strictly applied, but the reflection is done
+//   by reflection the high point, since we minimize a function
+//   instead of maximizing it, like Spendley.
+//   The Rule 2 is NOT implemented, as we expect that the
+//   function evaluation is not subject to errors.
+//   The Rule 3 is applied, ie reflection with respect
+//   to next to high point.
+//   A shrink step is included, with shrinkage factor sigma.
+//
+//   "Rule 1. Ascertain the lowest reading y, of yi ... Yk+1
+//   Complete a new simplex Sp by excluding the point Vp corresponding to
+//   y, and replacing it by V* defined as above."
+//   "Rule 2. If a result has occurred in (k + 1) successive simplexes, and is not
+//   then eliminated by application of Rule 1, do not move in the direction
+//   indicated by Rule 1, or at all, but discard the result and replace it by a new
+//   observation at the same point."
+//   "Rule 3. If y is the lowest reading in So , and if the next observation made,
+//   y* , is the lowest reading in the new simplex S , do not apply Rule 1 and
+//   return to So from Sp . Move out of S, by rejecting the second lowest reading
+//   (which is also the second lowest reading in So)."
 //
 function this = neldermead_fixed (this)
   //
@@ -266,6 +289,12 @@ function this = neldermead_fixed (this)
   newfvmean = optimsimplex_fvmean ( simplex );
   currentxopt = optimbase_cget ( this.optbase , "-x0" );
   //
+  // Set indices for "clarity"
+  //
+  ilow = 1
+  ihigh = n + 1
+  inext = n
+  //
   // Initialize
   //
   terminate = 0;
@@ -276,10 +305,10 @@ function this = neldermead_fixed (this)
   while (terminate == 0)
     this.optbase = optimbase_incriter ( this.optbase );
     iter = iter + 1;
-    xlow = optimsimplex_getx ( simplex , 1 )
-    flow = optimsimplex_getfv ( simplex , 1 )
-    xhigh = optimsimplex_getx ( simplex , n+1 )
-    fhigh = optimsimplex_getfv ( simplex , n+1 )
+    xlow = optimsimplex_getx ( simplex , ilow )
+    flow = optimsimplex_getfv ( simplex , ilow )
+    xhigh = optimsimplex_getx ( simplex , ihigh )
+    fhigh = optimsimplex_getfv ( simplex , ihigh )
     //
     // Store history
     //
@@ -332,18 +361,31 @@ function this = neldermead_fixed (this)
     // Reflect the worst point with respect to center
     //
     xr = neldermead_interpolate ( xbar , xhigh , this.rho );
-    [ this ,fr] = neldermead_function ( this ,xr);
+    [ this , fr ] = neldermead_function ( this , xr );
     this = neldermead_log (this,sprintf("xr="+strcat(string(xr)," ")+", f(xr)=%f",fr));
     //
     // Replace worst point by xr if it is better
     //
     if ( fr < fhigh ) then
       this = neldermead_log (this,sprintf("  > Perform reflect"));
-      simplex = optimsimplex_setve ( simplex , n+1 , fr , xr )
+      simplex = optimsimplex_setve ( simplex , ihigh , fr , xr )
     else
-      //  Shrink
-      this = neldermead_log (this,sprintf("  > Perform Shrink"));
-      [ simplex , this ] = optimsimplex_shrink ( simplex , neldermead_costf , this.sigma , this )
+      // Reflect / xnext
+      xnext = optimsimplex_getx ( simplex , inext );
+      fnext = optimsimplex_getfv ( simplex , inext );
+      xbar2 = optimsimplex_xbar ( simplex , inext );
+      this = neldermead_log (this,sprintf("xbar2="+strcat(string(xbar2)," ")+""));
+      xr2 = neldermead_interpolate ( xbar2 , xnext , this.rho );
+      [ this , fr2 ] = neldermead_function ( this ,xr2 );
+      this = neldermead_log (this,sprintf("xr2="+strcat(string(xr2)," ")+", f(xr2)=%f",fr2));
+      if ( fr2 < fnext ) then
+        this = neldermead_log (this,sprintf("  > Perform reflect / next"));
+        simplex = optimsimplex_setve ( simplex , inext , fr2 , xr2 )
+      else
+        //  Shrink
+        this = neldermead_log (this,sprintf("  > Perform Shrink"));
+        [ simplex , this ] = optimsimplex_shrink ( simplex , neldermead_costf , this.sigma , this )
+      end
     end
     //
     // Sort simplex
@@ -670,7 +712,10 @@ function this = neldermead_startup (this)
         error(errmsg);
       end
       if ( or ( x <> xp ) ) then
-        [ this , fv ] = neldermead_function ( this , xp );
+        // Set the index, so that, if an additionnal cost function argument is provided,
+        // it can be appended at the end.
+        index = 1;
+        [ this , fv ] = neldermead_function ( this , xp , index );
         // Transpose xp, which is a column vector
         simplex0 = optimsimplex_setve ( simplex0 , ive , fv , xp.' );
       end
@@ -896,7 +941,14 @@ function this = neldermead_box ( this )
         break
       end
     end
-    [ this , fr ] = neldermead_function ( this , xr );
+    if ( nbnlc > 0 ) then
+      // Set the index, so that, if an additionnal cost function argument is provided,
+      // it can be appended at the end.
+      index = 1;
+      [ this , fr ] = neldermead_function ( this , xr , index );
+    else
+      [ this , fr ] = neldermead_function ( this , xr );
+    end
     this = neldermead_log (this,sprintf("xr=[%s], f(xr)=%f", strcat(string(xr)," ") , fr));
     if ( fr >= flow & fr < fn ) then
       this = neldermead_log (this,sprintf("  > Perform reflection"));
@@ -913,7 +965,14 @@ function this = neldermead_box ( this )
           break
         end
       end
-      [ this ,fe] = neldermead_function ( this ,xe);
+      if ( nbnlc > 0 ) then
+        // Set the index, so that, if an additionnal cost function argument is provided,
+        // it can be appended at the end.
+        index = 1;
+        [ this , fe ] = neldermead_function ( this , xe , index );
+      else
+        [ this , fe ] = neldermead_function ( this , xe );
+      end
       this = neldermead_log (this,sprintf("xe=[%s], f(xe)=%f", strcat(string(xe)," ") , fe ));
       if (fe < fr) then
         this = neldermead_log (this,sprintf("  > Perform Expansion"));
@@ -934,7 +993,14 @@ function this = neldermead_box ( this )
           break
         end
       end
-      [ this ,fc] = neldermead_function ( this ,xc);
+      if ( nbnlc > 0 ) then
+        // Set the index, so that, if an additionnal cost function argument is provided,
+        // it can be appended at the end.
+        index = 1;
+        [ this , fc ] = neldermead_function ( this , xc , index );
+      else
+        [ this , fc ] = neldermead_function ( this , xc );
+      end
       this = neldermead_log (this,sprintf("xc=[%s], f(xc)=%f", strcat(string(xc)," ") , fc));
       if ( fc <= fr ) then
         this = neldermead_log (this,sprintf("  > Perform Outside Contraction"));
@@ -957,7 +1023,14 @@ function this = neldermead_box ( this )
           break
         end
       end
-      [ this ,fc] = neldermead_function ( this ,xc);
+      if ( nbnlc > 0 ) then
+        // Set the index, so that, if an additionnal cost function argument is provided,
+        // it can be appended at the end.
+        index = 1;
+        [ this , fc ] = neldermead_function ( this , xc , index );
+      else
+        [ this , fc ] = neldermead_function ( this , xc );
+      end
       this = neldermead_log (this,sprintf("xc=[%s], f(xc)=%f", strcat(string(xc)," ") , fc));
       if ( fc < fhigh ) then
         this = neldermead_log (this,sprintf("  > Perform Inside Contraction"));
index 98aab9b..34f2af5 100644 (file)
@@ -20,37 +20,18 @@ function [ this , xdata , ydata , zdata ] = nmplot_contour ( this , xmin , xmax
   // Check that there are only 2 parameters
   n = neldermead_cget ( this.nmbase , "-numberofvariables" );
   if n <> 2 then
-    errmsg = msprintf(gettext("%s: Unexpected number of variables %d. Cannot draw contour plot for functions which do not have 2 parameters."),"nmplot_contour",n)
+    errmsg = msprintf(gettext("%s: Unexpected number of variables %d. Cannot draw contour plot for functions which do not have two parameters."),"nmplot_contour",n)
     error(errmsg)
   end
   stepx = (xmax - xmin)/nx
   xdata = xmin:stepx:xmax;
   stepy = (ymax - ymin)/ny
   ydata = ymin:stepy:ymax;
-  nexp = length(xdata) * length(ydata)
-  iexp = 1;
   for ix = 1:length(xdata)
     for iy = 1:length(ydata)
-      x(iexp,1:2) = [xdata(ix) ydata(iy)];
-      iexp = iexp + 1;
-    end
-  end
-  // 2. Perform the experiments
-  // Vectorize the call to the cost function if possible
-  vectorizefunction = neldermead_cget ( this.nmbase , "-vectorizefunction" );
-  if vectorizefunction == 0 then
-    for iexp = 1:nexp
-      [ this.nmbase , f(iexp) ] = neldermead_function ( this.nmbase , x(iexp,1:2) );
-    end
-  else
-    [ this.nmbase , f ] = neldermead_function ( this.nmbase , x );
-  end
-  // 3. Store the experiments results in a matrix, suitable for input to the contour function
-  iexp = 1;
-  for ix = 1:length(xdata)
-    for iy = 1:length(ydata)
-      zdata ( ix , iy ) = f(iexp);
-      iexp = iexp + 1;
+      experiment = [xdata(ix) ydata(iy)];
+      [ this.nmbase , fiexp ] = neldermead_function ( this.nmbase , experiment );
+      zdata ( ix , iy ) = fiexp;
     end
   end
 endfunction
index 703d348..7fca7bf 100644 (file)
@@ -24,16 +24,12 @@ function optimsimplex_print ( this )
   elseif this.fv == [] then
     mprintf("Empty simplex (zero function values)\n");
   else
-  mprintf("Dimension : %d\n" , this.n );
-  mprintf("Number of vertices : %d\n" , this.nbve );
-  for k = 1:this.nbve
-    // Compute a string for x
-    ss = sprintf("%f", this.x(k,1));
-    for i = 2:this.n
-      ss = ss + " " + sprintf("%f", this.x(k,i));
+    mprintf("Dimension : %d\n" , this.n );
+    mprintf("Number of vertices : %d\n" , this.nbve );
+    str = optimsimplex_tostring ( this );
+    for k = 1:this.nbve
+      mprintf("%s\n" , str(k) );
     end
-    mprintf("Vertex #%d/%d : fv=%f, x=%s\n" , k , this.n+1 , this.fv(k), ss );
-  end
   end
 endfunction
 
index 44ac919..cf67dab 100644 (file)
@@ -27,11 +27,11 @@ function str = optimsimplex_tostring ( this )
   str = []
   for k = 1:this.nbve
     // Compute a string for x
-    ss = sprintf("%f", this.x(k,1));
+    ss = sprintf("%e", this.x(k,1));
     for i = 2:this.n
-      ss = ss + " " + sprintf("%f", this.x(k,i));
+      ss = ss + " " + sprintf("%e", this.x(k,i));
     end
-    str(k) = sprintf("Vertex #%d/%d : fv=%f, x=%s\n" , k , this.n+1 , this.fv(k), ss );
+    str(k) = sprintf("Vertex #%d/%d : fv=%e, x=%s\n" , k , this.nbve , this.fv(k), ss );
   end
   end
 endfunction