1 ## Copyright (C) 2002 Etienne Grossmann <etienne@egdn.net>
3 ## This program is free software; you can redistribute it and/or modify it under
4 ## the terms of the GNU General Public License as published by the Free Software
5 ## Foundation; either version 3 of the License, or (at your option) any later
8 ## This program is distributed in the hope that it will be useful, but WITHOUT
9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 ## You should have received a copy of the GNU General Public License along with
14 ## this program; if not, see <http://www.gnu.org/licenses/>.
16 ## ok = test_minimize - Test that minimize works
18 ok = 1; # Remains set if all ok. Set to 0 otherwise
19 cnt = 0; # Test counter
20 page_screen_output (0);
21 page_output_immediately (1);
23 if ! exist ("verbose"), verbose = 0; end
31 function v = ff(x,y,t)
32 A = [1 -1;1 1]; M = A'*diag([100,1])*A;
33 v = (x(1:2) - y(1:2))'*M*(x(1:2)-y(1:2)) + 1;
36 ## Return differential
37 function dv = dff(x,y,t)
38 if nargin < 3, t = 1; end
39 if t == 1, N = length (x); else N = length (y); end
40 A = [1 -1;1 1]; M = A'*diag([100,1])*A;
41 dv = 2*(x(1:2)-y(1:2))'*M;
42 if N>2, dv = [dv, zeros(1,N-2)]; end
43 if t == 2, dv = -dv; end
46 ## Return value, diff and 2nd diff
47 function [v,dv,d2v] = d2ff(x,y,t)
48 if nargin < 3, t = 1; end
49 if t == 1, N = length (x); else N = length (y); end
50 A = [1 -1;1 1]; M = A'*diag([100,1])*A;
51 v = (x(1:2) - y(1:2))'*M*(x(1:2)-y(1:2)) + 1;
52 dv = 2*(x(1:2)-y(1:2))'*M;
53 d2v = zeros (N); d2v(1:2,1:2) = 2*M;
54 if N>2, dv = [dv, zeros(1,N-2)]; end
55 if t == 2, dv = -dv; end
58 ## Return value, diff and inv of 2nd diff
59 function [v,dv,d2v] = d2iff(x,y,t)
60 if nargin < 3, t = 1; end
61 if t == 1, N = length (x); else N = length (y); end
62 A = [1 -1;1 1]; M = A'*diag([100,1])*A;
63 v = (x(1:2) - y(1:2))'*M*(x(1:2)-y(1:2)) + 1;
64 dv = 2*(x(1:2)-y(1:2))'*M;
65 d2v = zeros (N); d2v(1:2,1:2) = inv (2*M);
66 if N>2, dv = [dv, zeros(1,N-2)]; end
67 if t == 2, dv = -dv; end
71 function prn (varargin), printf (varargin{:}); fflush (stdout); end
75 prn ("\n Testing that minimize() works as it should\n\n");
76 prn (" Nparams = N = %i\n",N);
80 ## Plain run, just to make sure ######################################
81 ## Minimum wrt 'x' is y0
82 [xlev,vlev,nlev] = minimize ("ff",{x0,y0,1});
85 if max (abs (xlev-y0)) > 100*sqrt (eps)
87 prn ("Error is too big : %8.3g\n", max (abs (xlev-y0)));
90 elseif verbose, prn ("ok %i\n",cnt);
93 ## See what 'backend' gives in that last case ########################
94 [method,ctl] = minimize ("ff",{x0,y0,1},"order",0,"backend");
97 if ! ischar (method) || ! strcmp (method,"nelder_mead_min")
100 prn ("Wrong method '%s' != 'nelder_mead_min' was chosen\n", method);
102 prn ("minimize pretends to use a method that isn't a string\n");
107 elseif verbose, prn ("ok %i\n",cnt);
110 [xle2,vle2,nle2] = feval (method, "ff", {x0,y0,1}, ctl);
112 # nelder_mead_min is not very repeatable
113 # because of restarts from random positions
114 if max (abs (xlev-xle2)) > 100*sqrt (eps)
116 prn ("Error is too big : %8.3g\n", max (abs (xlev-xle2)));
119 elseif verbose, prn ("ok %i\n",cnt);
123 ## Run, w/ differential, just to make sure ###########################
124 ## Minimum wrt 'x' is y0
126 # [xlev,vlev,nlev] = minimize ("ff",{x0,y0,1},"df","dff");
129 # if max (abs (xlev-y0)) > 100*sqrt (eps)
131 # prn ("Error is too big : %8.3g\n", max (abs (xlev-y0)));
134 # elseif verbose, prn ("ok %i\n",cnt);
137 ## Run, w/ differential returned by function ('jac' option) ##########
138 ## Minimum wrt 'x' is y0
139 # [xlev,vlev,nlev] = minimize ("d2ff",{x0,y0,1},"jac");
142 # if max (abs (xlev-y0)) > 100*sqrt (eps)
144 # prn ("Error is too big : %8.3g\n", max (abs (xlev-y0)));
147 # elseif verbose, prn ("ok %i\n",cnt);
150 ## Run, w/ 2nd differential, just to make sure #######################
151 ## Minimum wrt 'x' is y0
152 [xlev,vlev,nlev] = minimize ("ff",{x0,y0,1},"d2f","d2ff");
155 if max (abs (xlev-y0)) > 100*sqrt (eps)
157 prn ("Error is too big : %8.3g\n", max (abs (xlev-y0)));
160 elseif verbose, prn ("ok %i\n",cnt);
163 ## Use the 'hess' option, when f can return 2nd differential #########
164 ## Minimum wrt 'x' is y0
165 [xlev,vlev,nlev] = minimize ("d2ff", {x0,y0,1},"hess");
168 if max (abs (xlev-y0)) > 100*sqrt (eps)
170 prn ("Error is too big : %8.3g\n", max (abs (xlev-y0)));
173 elseif verbose, prn ("ok %i\n",cnt);
176 ## Run, w/ inverse of 2nd differential, just to make sure ############
177 ## Minimum wrt 'x' is y0
178 [xlev,vlev,nlev] = minimize ("ff", {x0,y0,1},"d2i","d2iff");
181 if max (abs (xlev-y0)) > 100*sqrt (eps)
183 prn ("Error is too big : %8.3g\n", max (abs (xlev-y0)));
186 elseif verbose, prn ("ok %i\n",cnt);
189 ## Use the 'ihess' option, when f can return pinv of 2nd differential
190 ## Minimum wrt 'x' is y0
191 [xlev,vlev,nlev] = minimize ("d2iff", {x0,y0,1},"ihess");
194 if max (abs (xlev-y0)) > 100*sqrt (eps)
196 prn ("Error is too big : %8.3g\n", max (abs (xlev-y0)));
199 elseif verbose, prn ("ok %i\n",cnt);
202 ## Run, w/ numerical differential ####################################
203 ## Minimum wrt 'x' is y0
204 [xlev,vlev,nlev] = minimize ("ff",{x0,y0,1},"ndiff");
207 if max (abs (xlev-y0)) > 100*sqrt (eps)
209 prn ("Error is too big : %8.3g\n", max (abs (xlev-y0)));
212 elseif verbose, prn ("ok %i\n",cnt);
215 ## Run, w/ numerical differential, specified by "order" ##############
216 ## Minimum wrt 'x' is y0
217 [xlev,vlev,nlev] = minimize ("ff",{x0,y0,1},"order",1);
220 if max (abs (xlev-y0)) > 100*sqrt (eps)
222 prn ("Error is too big : %8.3g\n", max (abs (xlev-y0)));
225 elseif verbose, prn ("ok %i\n",cnt);
228 # ## See what 'backend' gives in that last case ########################
229 # [method,ctl] = minimize ("ff",{x0,y0,1},"order",1,"backend");
232 # if ! strcmp (method,"bfgsmin")
234 # prn ("Wrong method '%s' != 'bfgsmin' was chosen\n", method);
237 # elseif verbose, prn ("ok %i\n",cnt);
240 ## [xle2,vle2,nle2] = feval (method, "ff",{x0,y0,1}, ctl);
241 [xle2,vle2,nle2] = minimize ("ff",{x0,y0,1},"order",1);
243 if max (abs (xlev-xle2)) > 100*eps
245 prn ("Error is too big : %8.3g\n", max (abs (xlev-y0)));
248 elseif verbose, prn ("ok %i\n",cnt);
253 prn ( "All tests ok\n");