Renamed mtlb_fminsearch into fminsearch. Added demos.
[scilab.git] / scilab / modules / optimization / demos / neldermead / nmplot_mckinnon2.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
13 //% MCKINNON computes the McKinnon function.
14 //
15 //  Discussion:
16 //
17 //    This function has a global minimizer:
18 //
19 //      X* = ( 0.0, -0.5 ), F(X*) = -0.25
20 //
21 //    There are three parameters, TAU, THETA and PHI.
22 //
23 //    1 < TAU, then F is strictly convex.
24 //             and F has continuous first derivatives.
25 //    2 < TAU, then F has continuous second derivatives.
26 //    3 < TAU, then F has continuous third derivatives.
27 //
28 //    However, this function can cause the Nelder-Mead optimization
29 //    algorithm to "converge" to a point which is not the minimizer
30 //    of the function F.
31 //
32 //    Sample parameter values which cause problems for Nelder-Mead 
33 //    include:
34 //
35 //      TAU = 1, THETA = 15, PHI =  10;
36 //      TAU = 2, THETA =  6, PHI =  60;
37 //      TAU = 3, THETA =  6, PHI = 400;
38 //
39 //    To get the bad behavior, we also assume the initial simplex has the form
40 //
41 //      X1 = (0,0),
42 //      X2 = (1,1),
43 //      X3 = (A,B), 
44 //
45 //    where 
46 //
47 //      A = (1+sqrt(33))/8 =  0.84307...
48 //      B = (1-sqrt(33))/8 = -0.59307...
49 //
50 //  Licensing:
51 //
52 //    This code is distributed under the GNU LGPL license.
53 //
54 //  Modified:
55 //
56 //    09 February 2008
57 //
58 //  Author:
59 //
60 //    John Burkardt
61 //
62 //  Reference:
63 //
64 //    Ken McKinnon,
65 //    Convergence of the Nelder-Mead simplex method to a nonstationary point,
66 //    SIAM Journal on Optimization,
67 //    Volume 9, Number 1, 1998, pages 148-158.
68 //
69 //  Parameters:
70 //
71 //    Input, real X(2), the argument of the function.
72 //
73 //    Output, real F, the value of the function at X.
74 //
75 // Copyright (C) 2009 - INRIA - Michael Baudin, Scilab port
76
77 function f = mckinnon3 ( x )
78
79   if ( length ( x ) ~= 2 )
80     error ( 'Error: function expects a two dimensional input\n' );
81   end
82
83   tau = 3.0;
84   theta = 6.0;
85   phi = 400.0;
86
87   if ( x(1) <= 0.0 )
88     f = theta * phi * abs ( x(1) ).^tau + x(2) * ( 1.0 + x(2) );
89   else
90     f = theta       *       x(1).^tau   + x(2) * ( 1.0 + x(2) );
91   end
92 endfunction
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,"-tolsimplexizerelative",1.e-6);
109 nm = nmplot_configure(nm,"-simplex0method","given");
110 nm = nmplot_configure(nm,"-coords0",coords0);
111 nm = nmplot_configure(nm,"-simplex0length",1.0);
112 nm = nmplot_configure(nm,"-method","variable");
113 nm = nmplot_configure(nm,"-verbose",0);
114 nm = nmplot_configure(nm,"-verbosetermination",0);
115 nm = nmplot_configure(nm,"-kelleystagnationflag",1);
116 nm = nmplot_configure(nm,"-restartflag",1);
117 nm = nmplot_configure(nm,"-restartdetection","kelley");
118 //
119 // Setup output files
120 //
121 nm = nmplot_configure(nm,"-simplexfn","mckinnon.history.restart.simplex.txt");
122 nm = nmplot_configure(nm,"-fbarfn","mckinnon.history.restart.fbar.txt");
123 nm = nmplot_configure(nm,"-foptfn","mckinnon.history.restart.fopt.txt");
124 nm = nmplot_configure(nm,"-sigmafn","mckinnon.history.restart.sigma.txt");
125 //
126 // Perform optimization
127 //
128 nm = nmplot_search(nm);
129 //
130 // Plot
131 //
132 [nm , xdata , ydata , zdata ] = nmplot_contour ( nm , xmin = -0.2 , xmax = 2.0 , ymin = -2.0 , ymax = 2.0 , nx = 100 , ny = 100 );
133 f = scf();
134 xset("fpf"," ")
135 contour ( xdata , ydata , zdata , 40 )
136 nmplot_simplexhistory ( nm );
137 xs2png(0,"mckinnon.history.restart.simplex.png");
138 f = scf();
139 nmplot_historyplot ( nm , "mckinnon.history.restart.fbar.txt" , ...
140   mytitle = "Function Value Average" , myxlabel = "Iterations" );
141 xs2png(1,"mckinnon.history.restart.fbar.png");
142 f = scf();
143 nmplot_historyplot ( nm , "mckinnon.history.restart.fopt.txt" , ...
144   mytitle = "Minimum Function Value" , myxlabel = "Iterations" );
145 xs2png(2,"mckinnon.history.restart.fopt.png");
146 f = scf();
147 nmplot_historyplot ( nm , "mckinnon.history.restart.sigma.txt" , ...
148   mytitle = "Maximum Oriented length" , myxlabel = "Iterations" );
149 xs2png(3,"mckinnon.history.restart.sigma.png");
150 deletefile("mckinnon.history.restart.simplex.txt");
151 deletefile("mckinnon.history.restart.fbar.txt");
152 deletefile("mckinnon.history.restart.fopt.txt");
153 deletefile("mckinnon.history.restart.sigma.txt");
154 nm = nmplot_destroy(nm);
155
156