scilab/modules/optimization/demos/neldermead/nmplot_mckinnon.sce

   1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   2 // Copyright (C) 2008-2009 - INRIA - Michael Baudin
   3 //
   4 // This file must be used under the terms of the CeCILL.
   5 // This source file is licensed as described in the file COPYING, which
   6 // you should have received as part of this distribution.  The terms
   7 // are also available at
   8 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
   9
  10
  11 //dirname = get_absolute_file_path('nmplot_mckinnon.sce');
  12 //% MCKINNON computes the McKinnon function.
  13 //
  14 //  Discussion:
  15 //
  16 //    This function has a global minimizer:
  17 //
  18 //      X* = ( 0.0, -0.5 ), F(X*) = -0.25
  19 //
  20 //    There are three parameters, TAU, THETA and PHI.
  21 //
  22 //    1 < TAU, then F is strictly convex.
  23 //             and F has continuous first derivatives.
  24 //    2 < TAU, then F has continuous second derivatives.
  25 //    3 < TAU, then F has continuous third derivatives.
  26 //
  27 //    However, this function can cause the Nelder-Mead optimization
  28 //    algorithm to "converge" to a point which is not the minimizer
  29 //    of the function F.
  30 //
  31 //    Sample parameter values which cause problems for Nelder-Mead
  32 //    include:
  33 //
  34 //      TAU = 1, THETA = 15, PHI =  10;
  35 //      TAU = 2, THETA =  6, PHI =  60;
  36 //      TAU = 3, THETA =  6, PHI = 400;
  37 //
  38 //    To get the bad behavior, we also assume the initial simplex has the form
  39 //
  40 //      X1 = (0,0),
  41 //      X2 = (1,1),
  42 //      X3 = (A,B),
  43 //
  44 //    where
  45 //
  46 //      A = (1+sqrt(33))/8 =  0.84307...
  47 //      B = (1-sqrt(33))/8 = -0.59307...
  48 //
  49 //  Licensing:
  50 //
  51 //    This code is distributed under the GNU LGPL license.
  52 //
  53 //  Modified:
  54 //
  55 //    09 February 2008
  56 //
  57 //  Author:
  58 //
  59 //    John Burkardt
  60 //
  61 //  Reference:
  62 //
  63 //    Ken McKinnon,
  64 //    Convergence of the Nelder-Mead simplex method to a nonstationary point,
  65 //    SIAM Journal on Optimization,
  66 //    Volume 9, Number 1, 1998, pages 148-158.
  67 //
  68 //  Parameters:
  69 //
  70 //    Input, real X(2), the argument of the function.
  71 //
  72 //    Output, real F, the value of the function at X.
  73 //
  74 // Copyright (C) 2009 - INRIA - Michael Baudin, Scilab port
  75
  76 function f = mckinnon3 ( x )
  77
  78   if ( length ( x ) ~= 2 )
  79     error ( 'Error: function expects a two dimensional input\n' );
  80   end
  81
  82   tau = 3.0;
  83   theta = 6.0;
  84   phi = 400.0;
  85
  86   if ( x(1) <= 0.0 )
  87     f = theta * phi * abs ( x(1) ).^tau + x(2) * ( 1.0 + x(2) );
  88   else
  89     f = theta       *       x(1).^tau   + x(2) * ( 1.0 + x(2) );
  90   end
  91 endfunction
  92
  93
  94 lambda1 = (1.0 + sqrt(33.0))/8.0
  95 lambda2 = (1.0 - sqrt(33.0))/8.0
  96 coords0 = [
  97 1.0 0.0 lambda1
  98 1.0 0.0 lambda2
  99 ]
 100
 101
 102 nm = nmplot_new ();
 103 nm = nmplot_configure(nm,"-numberofvariables",2);
 104 nm = nmplot_configure(nm,"-function",mckinnon3);
 105 nm = nmplot_configure(nm,"-x0",[1.0 1.0]');
 106 nm = nmplot_configure(nm,"-maxiter",200);
 107 nm = nmplot_configure(nm,"-maxfunevals",300);
 108 nm = nmplot_configure(nm,"-tolfunrelative",10*%eps);
 109 nm = nmplot_configure(nm,"-tolxrelative",10*%eps);
 110 nm = nmplot_configure(nm,"-simplex0method","given");
 111 nm = nmplot_configure(nm,"-coords0",coords0);
 112 nm = nmplot_configure(nm,"-simplex0length",1.0);
 113 nm = nmplot_configure(nm,"-method","variable");
 114 //nm = nmplot_configure(nm,"-verbose",1);
 115 nm = nmplot_configure(nm,"-verbosetermination",1);
 116 //
 117 // Setup output files
 118 //
 119 nm = nmplot_configure(nm,"-simplexfn","mckinnon.history.simplex.txt");
 120 nm = nmplot_configure(nm,"-fbarfn","mckinnon.history.fbar.txt");
 121 nm = nmplot_configure(nm,"-foptfn","mckinnon.history.fopt.txt");
 122 nm = nmplot_configure(nm,"-sigmafn","mckinnon.history.sigma.txt");
 123 //
 124 // Perform optimization
 125 //
 126 nm = nmplot_search(nm);
 127 //
 128 // Plot
 129 //
 130 [nm , xdata , ydata , zdata ] = nmplot_contour ( nm , xmin = -0.2 , xmax = 1.2 , ymin = -2.0 , ymax = 2.0 , nx = 50 , ny = 50 );
 131 f = scf();
 132 xset("fpf"," ")
 133 contour ( xdata , ydata , zdata , 40 )
 134 nmplot_simplexhistory ( nm );
 135 xs2png(0,"mckinnon.history.simplex.png");
 136 f = scf();
 137 nmplot_historyplot ( nm , "mckinnon.history.fbar.txt" , ...
 138   mytitle = "Function Value Average" , myxlabel = "Iterations" );
 139 xs2png(1,"mckinnon.history.fbar.png");
 140 f = scf();
 141 nmplot_historyplot ( nm , "mckinnon.history.fopt.txt" , ...
 142   mytitle = "Minimum Function Value" , myxlabel = "Iterations" );
 143 xs2png(2,"mckinnon.history.fopt.png");
 144 f = scf();
 145 nmplot_historyplot ( nm , "mckinnon.history.sigma.txt" , ...
 146   mytitle = "Maximum Oriented length" , myxlabel = "Iterations" );
 147 xs2png(3,"mckinnon.history.sigma.png");
 148 deletefile("mckinnon.history.simplex.txt");
 149 deletefile("mckinnon.history.fbar.txt");
 150 deletefile("mckinnon.history.fopt.txt");
 151 deletefile("mckinnon.history.sigma.txt");
 152 nm = nmplot_destroy(nm);
 153
 154