Michaël Baudin [Fri, 23 Oct 2009 14:37:16 +0000 (16:37 +0200)]

index b4acff3..513f699 100644 (file)
@@ -14,7 +14,7 @@
//

mprintf("Illustrates Box'' algorithm on Box problem B.\n");
-mprintf("Defining Box Problem A function...\n");
+mprintf("Defining Box Problem B function...\n");

//
//  Reference:
@@ -57,10 +57,10 @@ function [ f , c , index ] = boxproblemB ( x , index )
f = []
c = []
x3 = x(1) + sqrt(3.0) * x(2)
-  if ( index==1 | index==3 ) then
+  if ( index==2 | index==6 ) then
f = -(9.0 - (x(1) - 3.0) ^ 2) * x(2) ^ 3 / 27.0 / sqrt(3.0)
end
-  if ( index==2 | index==3 ) then
+  if ( index==5 | index==6 ) then
c1 = x(1) / sqrt(3.0) - x(2)
c2 = x3
c3 = 6.0 - x3
@@ -77,18 +77,18 @@ rand("seed" , 0)

x0 = [1.0 0.5].';
// Compute f(x0) : should be close to -0.0133645895646
-fx0 = boxproblemB ( x0 );
+fx0 = boxproblemB ( x0 , 2 );
mprintf("Computed fx0 = %e (expected = %e)\n",fx0 , -0.0133645895646 );
-result = boxproblemB ( x0 , 2 );
+[fx0 , cx0 , index ] = boxproblemB ( x0 , 6 );
mprintf("Computed Constraints(x0) = [%e %e %e]\n", ...
-  result(1), result(2), result(3) );
+  cx0(1), cx0(2), cx0(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 );
+fopt = boxproblemB ( xopt , 2 );
mprintf("Computed fopt = %e (expected = %e)\n", fopt , -1.0 );

@@ -117,7 +117,7 @@ nm = neldermead_configure(nm,"-boxboundsalpha" , 0.0001 );
// The index must be provided, because the additionnal argument "data"
// comes after.
//
-[ nm , result ] = neldermead_function ( nm , x0 , 1 );
+[ nm , f ] = neldermead_function ( nm , x0 );
//
// Perform optimization
//
new file mode 100644 (file)
index 0000000..6863644
--- /dev/null
@@ -0,0 +1,85 @@
+// 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 the sensitivity to dimension of the Nelder-Mead algorithm\n");
+
+function [ f , index ] = quadracticn ( x , index )
+  f = sum(x.^2);
+endfunction
+//
+// solvepb --
+//   Find the solution for the given number of dimensions
+//
+function [nbfevals , niter , rho] = solvepb ( n )
+  rand("seed",0)
+  coords (1,1:n) = zeros(1,n);
+  for i = 2:n+1
+    coords (i,1:n) = 2.0 * rand(1,n) - 1.0;
+  end
+  //
+  // Perform optimization
+  //
+  si0 = neldermead_get ( nm , "-simplex0" );
+  sigma0 = optimsimplex_size ( si0 , "sigmaplus" );
+  siopt = neldermead_get ( nm , "-simplexopt" );
+  sigmaopt = optimsimplex_size ( siopt , "sigmaplus" );
+  niter = neldermead_get ( nm , "-iterations" );
+  rho = (sigmaopt/sigma0)^(1.0/niter);
+  nbfevals = neldermead_get ( nm , "-funevals" );
+  mprintf ( "%d %d %d %f\n", n , nbfevals , niter , rho );
+endfunction
+
+
+for n = 1:20
+  [nbfevals niter rho] = solvepb ( n );
+  array_rho(n) = rho;
+  array_nbfevals(n) = nbfevals;
+  array_niter(n) = niter;
+end
+// Plot rate of convergence
+hh = scf();
+plot(1:20,array_rho)
+hh.children.x_label.text = "Number of parameters"
+hh.children.y_label.text = "Rate of convergence"
+hh.children.children.children.mark_mode = "on";
+hh.children.children.children.mark_style = 9;
+hh.children.children.children.mark_size = 10;
+
+// Plot number of function evaluations
+hh = scf();
+plot(1:20,array_nbfevals)
+hh.children.x_label.text = "Number of parameters"
+hh.children.y_label.text = "Number of function evaluations"
+hh.children.children.children.mark_mode = "on";
+hh.children.children.children.mark_style = 9;
+hh.children.children.children.mark_size = 10;
+//
+// Load this script into the editor
+//
+dname = get_absolute_file_path(filename);
+editor ( dname + filename );
+
index ef715a4..4b82f2c 100644 (file)
@@ -116,8 +116,8 @@ nm = nmplot_configure(nm,"-simplex0length",1.0);
nm = nmplot_configure(nm,"-method","variable");
//nm = nmplot_configure(nm,"-verbose",0);
nm = nmplot_configure(nm,"-verbosetermination",0);
-nm = nmplot_configure(nm,"-kelleystagnationflag",1);
-nm = nmplot_configure(nm,"-restartflag",1);
+nm = nmplot_configure(nm,"-kelleystagnationflag",%t);
+nm = nmplot_configure(nm,"-restartflag",%t);
nm = nmplot_configure(nm,"-restartdetection","kelley");
//
// Setup output files
new file mode 100644 (file)
index 0000000..276f220
--- /dev/null
@@ -0,0 +1,89 @@
+// 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 variable-shape Nelder-Mead algorithm performs well on a quadratic test case.\n");
+function [ y , index ] = quadratic ( x , index )
+  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,"-x0",[2.0 2.0]');
+nm = nmplot_configure(nm,"-maxiter",100);
+nm = nmplot_configure(nm,"-maxfunevals",300);
+nm = nmplot_configure(nm,"-tolxmethod",%f);
+nm = nmplot_configure(nm,"-tolsimplexizerelative",1.e-8);
+nm = nmplot_configure(nm,"-simplex0method","spendley");
+nm = nmplot_configure(nm,"-method","variable");
+//nm = nmplot_configure(nm,"-verbose",1);
+//nm = nmplot_configure(nm,"-verbosetermination",1);
+//
+// Setup output files
+//
+nm = nmplot_configure(nm,"-simplexfn","history.simplex.txt");
+nm = nmplot_configure(nm,"-fbarfn","history.fbar.txt");
+nm = nmplot_configure(nm,"-foptfn","history.fopt.txt");
+nm = nmplot_configure(nm,"-sigmafn","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 , "history.fbar.txt" , ...
+  mytitle = "Function Value Average" , myxlabel = "Iterations" );
+mprintf("Plotting history of fopt...\n");
+f = scf();
+nmplot_historyplot ( nm , "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 , "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 = -2.0 , xmax = 4.0 , ymin = -2.0 , ymax = 4.0 , nx = 50 , ny = 50 );
+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("history.simplex.txt");
+deletefile("history.fbar.txt");
+deletefile("history.fopt.txt");
+deletefile("history.sigma.txt");
+nm = nmplot_destroy(nm);
+
+//
+// Load this script into the editor
+//
+dname = get_absolute_file_path(filename);
+editor ( dname + filename );
+
new file mode 100644 (file)
index 0000000..4c829ca
--- /dev/null
@@ -0,0 +1,45 @@
+// 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
+
+
+a = 100.0;
+function [ y , index ] = quadratic ( x , index )
+  y = a * 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,"-x0",[10.0 10.0]');
+nm = nmplot_configure(nm,"-maxiter",400);
+nm = nmplot_configure(nm,"-maxfunevals",400);
+nm = nmplot_configure(nm,"-tolxmethod",%f);
+nm = nmplot_configure(nm,"-tolsimplexizerelative",1.e-8);
+nm = nmplot_configure(nm,"-simplex0method","spendley");
+nm = nmplot_configure(nm,"-method","variable");
+//nm = nmplot_configure(nm,"-verbose",1);
+//nm = nmplot_configure(nm,"-verbosetermination",0);
+//
+// Perform optimization
+//
+mprintf("Searching for minimum...\n");
+nm = nmplot_search(nm);
+nmplot_display(nm);
+nm = nmplot_destroy(nm);
+
+//
+// Load this script into the editor
+//