remove reference to comp and remove variable type 11
[scilab.git] / scilab / modules / optimization / macros / neldermead / fminsearch.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 2008-2009 - INRIA - Michael Baudin
3 // Copyright (C) 2009-2011 - DIGITEO - Michael Baudin
4 //
5 // Copyright (C) 2012 - 2016 - Scilab Enterprises
6 //
7 // This file is hereby licensed under the terms of the GNU GPL v2.0,
8 // pursuant to article 5.3.4 of the CeCILL v.2.1.
9 // This file was originally licensed under the terms of the CeCILL v2.1,
10 // and continues to be available under such terms.
11 // For more information, see the COPYING file which you should have received
12 // along with this program.
13
14 //
15 // fminsearch --
16 //   Emulate the fminsearch command of Matlab.
17 //   Search the minimum with Nelder-Mead algorithm.
18 //   [x,fval,exitflag,output] = fminsearch(fun,x0,options)
19 // Arguments, input
20 //   fun : the function to minimize
21 //   x0 : a row vector with dimension n where n is the number of parameters
22 //        to optimize.
23 //        Initial guess for optimization algorithm.
24 //  options : an optional struct, as provided by optimset
25 //
26 function [x,fval,exitflag,output] = fminsearch ( varargin )
27     //
28     // The output function called back by fminsearch
29     //
30     // Arguments
31     //  state : the current state of the algorithm
32     //    "init", "iter", "done"
33     //  data : the data at the current state
34     //    This is a tlist with the following entries:
35     //    * x : the optimal vector of parameters
36     //    * fval : the minimum function value
37     //    * simplex : the simplex, as a simplex object
38     //    * iteration : the number of iterations performed
39     //    * funccount : the number of function evaluations
40     //    * step : the type of step in the previous iteration
41     //  fmsdata : this is a tlist which contains specific data of the
42     //    fminsearch algorithm
43     //    * Display : what to display
44     //    * OutputFcn : the array of output functions
45     //
46     function stop = fminsearch_outputfun ( state , data , fmsdata )
47         //
48         // Compute procedure
49         //
50         select data.step
51         case "init" then
52             if ( data.iteration == 0 ) then
53                 procedure = "";
54             else
55                 procedure = "initial simplex";
56             end
57         case "done" then
58             procedure = ""
59         case "reflection" then
60             procedure = "reflect"
61         case "expansion" then
62             procedure = "expand"
63         case "insidecontraction" then
64             procedure = "contract inside"
65         case "outsidecontraction" then
66             procedure = "contract outside"
67         case "shrink" then
68             procedure = "shrink"
69         else
70             errmsg = msprintf(gettext("%s: Unknown step %s"), "fminsearch", data.step)
71             error(errmsg)
72         end
73         //
74         // Display a message
75         //
76         if ( fmsdata.Display == "iter" ) then
77             if ( data.step <> "done" ) then
78                 mprintf ( "%6s        %5s     %12s         %-20s\n", ...
79                 string(data.iteration) , string(data.funccount) , string(data.fval) , procedure )
80             else
81                 mprintf ( "\n" )
82             end
83         end
84         //
85         // Process output functions
86         //
87         stop = %f
88         optimValues = struct(...
89         "funccount" ,data.funccount , ...
90         "fval" ,data.fval , ...
91         "iteration" , data.iteration , ...
92         "procedure" , procedure ...
93         );
94         if ( fmsdata.OutputFcn <> [] ) then
95             if type ( fmsdata.OutputFcn ) == 13 then
96                 // The output function is a macro
97                 stop = fmsdata.OutputFcn ( data.x , optimValues , state );
98                 //
99                 // Backward-compatibility: define the stop variable
100                 //
101                 if ( exists("stop")==0 ) then
102                     fms_warnheaderobsolete ( "outputfun(x,optimValues , state )" , "stop=outputfun(x,optimValues , state )", "5.4.1" )
103                     stop = %f
104                 end
105             elseif ( type ( fmsdata.OutputFcn ) == 15 ) then
106                 // The output function is a list of macros
107                 for i = 1:length(fmsdata.OutputFcn)
108                     stop = fmsdata.OutputFcn(i) ( data.x , optimValues , state );
109                 end
110                 //
111                 // Backward-compatibility: define the stop variable
112                 //
113                 if ( exists("stop")==0 ) then
114                     fms_warnheaderobsolete ( "outputfun(x,optimValues , state )" , "stop=outputfun(x,optimValues , state )", "5.4.1" )
115                     stop = %f
116                 end
117             else
118                 // The user did something wrong...
119                 errmsg = msprintf(gettext("%s: The value of the ''OutputFcn'' option is neither a function nor a list."), "fminsearch")
120                 error(errmsg)
121             end
122         end
123         // Process plot functions
124         if ( fmsdata.PlotFcns <> [] ) then
125             if type ( fmsdata.PlotFcns ) == 13 then
126                 // The output function is a macro
127                 fmsdata.PlotFcns ( data.x , optimValues , state );
128             elseif ( type ( fmsdata.PlotFcns ) == 15 ) then
129                 // The output function is a list of macros
130                 for i = 1:length(fmsdata.PlotFcns)
131                     fmsdata.PlotFcns(i) ( data.x , optimValues , state );
132                 end
133             else
134                 // The user did something wrong...
135                 errmsg = msprintf(gettext("%s: The value of the ''PlotFcns'' option is neither a function nor a list."), "fminsearch")
136                 error(errmsg)
137             end
138         end
139     endfunction
140     //
141     // fminsearch_function --
142     //   Calls the cost function and make it match
143     //   neldermead requirements.
144     //
145     function [ f , index , fmsfundata ] = fminsearch_function ( x , index , fmsfundata )
146         funtype = typeof(fmsfundata.Fun)
147         if ( funtype == "function" ) then
148             __fminsearch_f__ = fmsfundata.Fun
149             __fminsearch_args__ = list()
150         else
151             __fminsearch_f__ = fmsfundata.Fun(1)
152             __fminsearch_args__ = list(fmsfundata.Fun(2:$))
153         end
154         f = __fminsearch_f__ ( x , __fminsearch_args__(:))
155     endfunction
156
157     function fms_warnheaderobsolete ( oldheader , newheader , removedVersion )
158         warnMessage = msprintf(_("Syntax %s is obsolete."),oldheader)
159         warnMessage = [warnMessage, msprintf(_("Please use %s instead."),newheader)]
160         warnMessage = [warnMessage, msprintf(_("This feature will be permanently removed in Scilab %s"), removedVersion)]
161         warning(warnMessage);
162     endfunction
163
164     function errMessage = fms_errheaderobsolete (oldheader, newheader)
165         errMessage = msprintf(_("Calling sequence %s is obsolete."),oldheader)
166         errMessage = [errMessage, msprintf(_("Please use %s instead."),newheader)]
167     endfunction
168
169
170     function assert_typecallable ( var , varname , ivar )
171         // Check that var is a function or a list
172         if ( and ( type ( var ) <> [13 15] ) ) then
173             errmsg = msprintf(gettext("%s: Expected function or list for variable %s at input #%d, but got %s instead."),"assert_typecallable", varname , ivar , typeof(var) );
174             error(errmsg);
175         end
176         if ( type ( var ) == 15 ) then
177             // Check that var(1) is a function
178             if type ( var(1) ) <> 13 then
179                 errmsg = msprintf(gettext("%s: Expected function for variable %s(1) at input #%d, but got %s instead."),"assert_typecallable", varname , ivar , typeof(var) );
180                 error(errmsg);
181             end
182         end
183     endfunction
184     // Generates an error if the given variable is not of type real
185     function assert_typereal ( var , varname , ivar )
186         if ( type ( var ) <> 1 ) then
187             errmsg = msprintf(gettext("%s: Expected real variable for variable %s at input #%d, but got %s instead."),"assert_typereal", varname , ivar , typeof(var) );
188             error(errmsg);
189         end
190     endfunction
191
192     [lhs,rhs]=argn();
193     if rhs<>2 & rhs<>3 then
194         errmsg = msprintf(gettext("%s: Wrong number of input arguments: %d or %d expected.\n"), "fminsearch", 2,3);
195         error(errmsg)
196     end
197     fun = varargin(1);
198     x0 = varargin(2);
199     // Get x0 and change it into a column vector
200     x0t = size(x0,"*");
201     x0 = matrix(x0,x0t,1);
202     defaultoptions = optimset ("fminsearch");
203     msg="";
204     if rhs==2 then
205         // No options on the command line
206         // Set default values
207         options = defaultoptions;
208     elseif rhs==3 then
209         // One options struc on the command line : use it !
210         options = varargin(3);
211     end
212     // Compute options from the options struct
213     numberofvariables = size(x0,"*");
214     MaxFunEvals = optimget ( options , "MaxFunEvals" , defaultoptions.MaxFunEvals );
215     MaxIter = optimget ( options , "MaxIter" , defaultoptions.MaxIter );
216     TolFun = optimget ( options , "TolFun" , defaultoptions.TolFun );
217     TolX = optimget ( options , "TolX" , defaultoptions.TolX );
218     Display = optimget ( options , "Display" , defaultoptions.Display );
219     OutputFcn = optimget ( options , "OutputFcn" , defaultoptions.OutputFcn );
220     PlotFcns = optimget ( options , "PlotFcns" , defaultoptions.PlotFcns );
221     // If the MaxIter option is a string, we make the assumption that it is the default 200 value.
222     // If not, this is the actual value.
223     if ( type ( MaxIter ) == 10 ) then
224         if ( MaxIter == "200*numberofvariables" ) then
225             MaxIter = 200 * numberofvariables;
226         else
227             errmsg = msprintf(gettext("%s: Unexpected maximum number of iterations %s."), "fminsearch", MaxIter );
228             error(errmsg)
229         end
230     end
231     // If the MaxFunEvals option is a string, this is the default 200 value
232     // If not, this is the actual value.
233     if ( type ( MaxFunEvals ) == 10 ) then
234         if ( MaxFunEvals == "200*numberofvariables" ) then
235             MaxFunEvals = 200 * numberofvariables;
236         else
237             errmsg = msprintf(gettext("%s: Unexpected maximum number of function evaluations %s."), "fminsearch", MaxFunEvals );
238             error(errmsg)
239         end
240     end
241     if ( Display == "iter" ) then
242         mprintf ( "%10s   %10s   %10s %17s\n" , "Iteration", "Func-count" , "min f(x)" , "Procedure" );
243     end
244
245     //check OutputFcn format
246     if type(OutputFcn) == 13 then
247         macroInfo = macrovar(OutputFcn);
248         if size(macroInfo(2), "*") <> 1 then
249             errMessage = fms_errheaderobsolete("outputfun(x,optimValues , state )", "stop=outputfun(x,optimValues , state )");
250             error(errMessage);
251         end
252     elseif type(OutputFcn) == 15 then
253         for i = 1 : size(OutputFcn)
254             if type(OutputFcn(i)) == 13 then
255                 macroInfo = macrovar(OutputFcn(i));
256                 if size(macroInfo(2), "*") <> 1 then
257                     errMessage = fms_errheaderobsolete("outputfun(x,optimValues , state )", "stop=outputfun(x,optimValues , state )");
258                     error(errMessage);
259                 end
260             end
261         end
262     elseif (OutputFcn <> [])
263         // The user did something wrong...
264         errmsg = msprintf(gettext("%s: The value of the ''OutputFcn'' option is neither a function nor a list."), "fminsearch");
265         error(errmsg)
266     end
267     //
268     // Check input arguments
269     assert_typecallable ( fun , "costf" , 1)
270     assert_typereal ( x0 , "x0" , 2 );
271     //
272     // Prepare the data structure to pass to the output function
273     fmsdata = tlist(["T_FMINSEARCH"
274     "Display"
275     "OutputFcn"
276     "PlotFcns"
277     ]);
278     fmsdata.Display = Display
279     fmsdata.OutputFcn = OutputFcn
280     fmsdata.PlotFcns = PlotFcns
281     // Prepare the data structure to pass to the cost function
282     fmsfundata = tlist(["T_FMINSEARCH"
283     "Fun"
284     ]);
285     fmsfundata.Fun = fun
286     // Perform Optimization
287     nm = neldermead_new ();
288     nm = neldermead_configure(nm,"-x0",x0);
289     nm = neldermead_configure(nm,"-numberofvariables",numberofvariables);
290     nm = neldermead_configure(nm,"-simplex0method","pfeffer");
291     nm = neldermead_configure(nm,"-simplex0deltausual",0.05);
292     nm = neldermead_configure(nm,"-simplex0deltazero",0.0075);
293     nm = neldermead_configure(nm,"-method","variable");
294     nm = neldermead_configure(nm,"-function",list(fminsearch_function,fmsfundata));
295     nm = neldermead_configure(nm,"-maxiter",MaxIter);
296     nm = neldermead_configure(nm,"-maxfunevals",MaxFunEvals);
297     nm = neldermead_configure(nm,"-tolxmethod",%f);
298     nm = neldermead_configure(nm,"-tolfunmethod",%f);
299     nm = neldermead_configure(nm,"-tolssizedeltafvmethod",%t);
300     nm = neldermead_configure(nm,"-tolsimplexizemethod",%f);
301     nm = neldermead_configure(nm,"-toldeltafv",TolFun);
302     nm = neldermead_configure(nm,"-tolsimplexizeabsolute",TolX);
303     nm = neldermead_configure(nm,"-checkcostfunction",%f);
304     nm = neldermead_configure(nm,"-outputcommand",list(fminsearch_outputfun,fmsdata));
305     //nm = neldermead_configure(nm,"-verbose",1);
306     //nm = neldermead_configure(nm,"-verbosetermination",1);
307     nm = neldermead_search(nm, "off");
308     x = neldermead_get(nm,"-xopt").';
309     fval = neldermead_get(nm,"-fopt");
310     status = neldermead_get(nm,"-status");
311     select status
312     case "maxiter" then
313         if ( ( Display == "notify" ) | ( Display == "iter" ) | ( Display == "final" ) ) then
314             msg = "%s: Exiting: Maximum number of iterations has been exceeded\n" + ...
315             "         - increase MaxIter option.\n" + ...
316             "         Current function value: %s\n"
317             mprintf(gettext(msg) , "fminsearch" , string(fval) )
318         end
319         exitflag = 0;
320     case "maxfuneval" then
321         if ( ( Display == "notify" ) | ( Display == "iter" ) | ( Display == "final" ) ) then
322             msg = "%s: Exiting: Maximum number of function evaluations has been exceeded\n" + ...
323             "          - increase MaxFunEvals option.\n" + ...
324             "         Current function value: %s\n"
325             mprintf(gettext(msg) , "fminsearch" , string(fval) )
326         end
327         exitflag = 0;
328     case "tolsizedeltafv" then
329         exitflag = 1;
330         msg = sprintf("%s\n%s %s\n%s %s", "Optimization terminated:",...
331         " the current x satisfies the termination criteria using OPTIONS.TolX of",...
332         string(TolX),...
333         " and F(X) satisfies the convergence criteria using OPTIONS.TolFun of",...
334         string(TolFun));
335     case "userstop" then
336         msg = sprintf("%s\n%s\n%s", "Optimization terminated:",...
337         " ",...
338         " User Stop");
339         exitflag = -1;
340     else
341         errmsg = msprintf(gettext("%s: Unknown status %s"), "fminsearch", status)
342         error(errmsg)
343     end
344     output = struct(...
345     "algorithm" ,[],...
346     "funcCount" ,[],...
347     "iterations" ,[],...
348     "message" , []);
349     output.algorithm = "Nelder-Mead simplex direct search";
350     output.funcCount = neldermead_get(nm,"-funevals");
351     output.iterations = neldermead_get(nm,"-iterations");
352     output.message = msg;
353     if ( ( Display == "final" ) | ( Display == "iter" ) ) then
354         if ( ( exitflag == 1 ) ) then
355             mprintf( "%s\n" , output.message(1) );
356             mprintf( "%s\n" , output.message(2) );
357             mprintf( "%s\n" , output.message(3) );
358         end
359     end
360     nm = neldermead_destroy(nm);
361 endfunction