b23b4ac64fa0c0aa2536083d8cedbfdeb348b83e
[scilab.git] / 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 mprintf("Defining McKinnon function...\n");
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 , index ] = mckinnon3 ( x , index )
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 1.0
98 0.0 0.0
99 lambda1 lambda2
100 ];
101
102
103 x0 = [1.0 1.0]';
104 mprintf("x0=%s\n",strcat(string(x0)," "));
105 mprintf("Creating object...\n");
106 nm = nmplot_new ();
107 nm = nmplot_configure(nm,"-numberofvariables",2);
108 nm = nmplot_configure(nm,"-function",mckinnon3);
109 nm = nmplot_configure(nm,"-x0",x0);
110 nm = nmplot_configure(nm,"-maxiter",200);
111 nm = nmplot_configure(nm,"-maxfunevals",300);
112 nm = nmplot_configure(nm,"-tolfunrelative",10*%eps);
113 nm = nmplot_configure(nm,"-tolxrelative",10*%eps);
114 nm = nmplot_configure(nm,"-simplex0method","given");
115 nm = nmplot_configure(nm,"-coords0",coords0);
116 nm = nmplot_configure(nm,"-simplex0length",1.0);
117 nm = nmplot_configure(nm,"-method","variable");
118 //nm = nmplot_configure(nm,"-verbose",1);
119 nm = nmplot_configure(nm,"-verbosetermination",1);
120 //
121 // Setup output files
122 //
123 nm = nmplot_configure(nm,"-simplexfn","mckinnon.history.simplex.txt");
124 nm = nmplot_configure(nm,"-fbarfn","mckinnon.history.fbar.txt");
125 nm = nmplot_configure(nm,"-foptfn","mckinnon.history.fopt.txt");
126 nm = nmplot_configure(nm,"-sigmafn","mckinnon.history.sigma.txt");
127 //
128 // Perform optimization
129 //
130 nm = nmplot_search(nm);
131 nmplot_display(nm);
132 //
133 // Plot
134 //
135 mprintf("Plot contour...\n");
136 [nm , xdata , ydata , zdata ] = nmplot_contour ( nm , xmin = -0.2 , xmax = 1.2 , ymin = -2.0 , ymax = 2.0 , nx = 50 , ny = 50 );
137 f = scf(100001);
138 xset("fpf"," ")
139 drawlater();
140 contour ( xdata , ydata , zdata , [-0.2 0.0 1.0 2.0 5.0 10.0 20.0] )
141 nmplot_simplexhistory ( nm );
142 drawnow();
143 f = scf(100002);
144 nmplot_historyplot ( nm , "mckinnon.history.fbar.txt" , ...
145   mytitle = "Function Value Average" , myxlabel = "Iterations" );
146 f = scf(100003);
147 nmplot_historyplot ( nm , "mckinnon.history.fopt.txt" , ...
148   mytitle = "Minimum Function Value" , myxlabel = "Iterations" );
149 f = scf(100004);
150 nmplot_historyplot ( nm , "mckinnon.history.sigma.txt" , ...
151   mytitle = "Maximum Oriented length" , myxlabel = "Iterations" );
152 deletefile("mckinnon.history.simplex.txt");
153 deletefile("mckinnon.history.fbar.txt");
154 deletefile("mckinnon.history.fopt.txt");
155 deletefile("mckinnon.history.sigma.txt");
156 nm = nmplot_destroy(nm);
157
158 //
159 // Load this script into the editor
160 //
161 filename = 'nmplot_mckinnon.sce';
162 dname = get_absolute_file_path(filename);
163 editor ( dname + filename );
164