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