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