octave_packages/optim-1.2.0/test_minimize_1.m

   1 ## Copyright (C) 2002 Etienne Grossmann <etienne@egdn.net>
   2 ##
   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
   6 ## version.
   7 ##
   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
  11 ## details.
  12 ##
  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/>.
  15
  16 ## ok = test_minimize           - Test that minimize works
  17
  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);
  22
  23 if ! exist ("verbose"), verbose = 0; end
  24
  25 N = 2;
  26
  27 x0 = randn(N,1) ;
  28 y0 = randn(N,1) ;
  29
  30 ## Return value
  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;
  34 endfunction
  35
  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
  44 endfunction
  45
  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
  56 endfunction
  57
  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
  68 endfunction
  69
  70 ## PRint Now
  71 function prn (varargin), printf (varargin{:}); fflush (stdout); end
  72
  73
  74 if verbose
  75   prn ("\n   Testing that minimize() works as it should\n\n");
  76   prn ("  Nparams = N = %i\n",N);
  77   fflush (stdout);
  78 end
  79
  80 ## Plain run, just to make sure ######################################
  81 ## Minimum wrt 'x' is y0
  82 [xlev,vlev,nlev] = minimize ("ff",{x0,y0,1});
  83
  84 cnt++;
  85 if max (abs (xlev-y0)) > 100*sqrt (eps)
  86   if verbose
  87     prn ("Error is too big : %8.3g\n", max (abs (xlev-y0)));
  88   end
  89   ok = 0;
  90 elseif verbose,  prn ("ok %i\n",cnt);
  91 end
  92
  93 ## See what 'backend' gives in that last case ########################
  94 [method,ctl] = minimize ("ff",{x0,y0,1},"order",0,"backend");
  95
  96 cnt++;
  97 if ! ischar (method) || ! strcmp (method,"nelder_mead_min")
  98   if verbose
  99     if ischar (method)
 100       prn ("Wrong method '%s' != 'nelder_mead_min' was chosen\n", method);
 101     else
 102       prn ("minimize pretends to use a method that isn't a string\n");
 103     end
 104     return
 105   end
 106   ok = 0;
 107 elseif verbose,  prn ("ok %i\n",cnt);
 108 end
 109
 110 [xle2,vle2,nle2] = feval (method, "ff", {x0,y0,1}, ctl);
 111 cnt++;
 112                                 # nelder_mead_min is not very repeatable
 113                                 # because of restarts from random positions
 114 if max (abs (xlev-xle2)) > 100*sqrt (eps)
 115   if verbose
 116     prn ("Error is too big : %8.3g\n", max (abs (xlev-xle2)));
 117   end
 118   ok = 0;
 119 elseif verbose,  prn ("ok %i\n",cnt);
 120 end
 121
 122
 123 ## Run, w/ differential, just to make sure ###########################
 124 ## Minimum wrt 'x' is y0
 125
 126 # [xlev,vlev,nlev] = minimize ("ff",{x0,y0,1},"df","dff");
 127
 128 # cnt++;
 129 # if max (abs (xlev-y0)) > 100*sqrt (eps)
 130 #   if verbose
 131 #     prn ("Error is too big : %8.3g\n", max (abs (xlev-y0)));
 132 #   end
 133 #   ok = 0;
 134 # elseif verbose,  prn ("ok %i\n",cnt);
 135 # en
 136
 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");
 140
 141 # cnt++;
 142 # if max (abs (xlev-y0)) > 100*sqrt (eps)
 143 #   if verbose
 144 #     prn ("Error is too big : %8.3g\n", max (abs (xlev-y0)));
 145 #   end
 146 #   ok = 0;
 147 # elseif verbose,  prn ("ok %i\n",cnt);
 148 # end
 149
 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");
 153
 154 cnt++;
 155 if max (abs (xlev-y0)) > 100*sqrt (eps)
 156   if verbose
 157     prn ("Error is too big : %8.3g\n", max (abs (xlev-y0)));
 158   end
 159   ok = 0;
 160 elseif verbose,  prn ("ok %i\n",cnt);
 161 end
 162
 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");
 166
 167 cnt++;
 168 if max (abs (xlev-y0)) > 100*sqrt (eps)
 169   if verbose
 170     prn ("Error is too big : %8.3g\n", max (abs (xlev-y0)));
 171   end
 172   ok = 0;
 173 elseif verbose,  prn ("ok %i\n",cnt);
 174 end
 175
 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");
 179
 180 cnt++;
 181 if max (abs (xlev-y0)) > 100*sqrt (eps)
 182   if verbose
 183     prn ("Error is too big : %8.3g\n", max (abs (xlev-y0)));
 184   end
 185   ok = 0;
 186 elseif verbose,  prn ("ok %i\n",cnt);
 187 end
 188
 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");
 192
 193 cnt++;
 194 if max (abs (xlev-y0)) > 100*sqrt (eps)
 195   if verbose
 196     prn ("Error is too big : %8.3g\n", max (abs (xlev-y0)));
 197   end
 198   ok = 0;
 199 elseif verbose,  prn ("ok %i\n",cnt);
 200 end
 201
 202 ## Run, w/ numerical differential ####################################
 203 ## Minimum wrt 'x' is y0
 204 [xlev,vlev,nlev] = minimize ("ff",{x0,y0,1},"ndiff");
 205
 206 cnt++;
 207 if max (abs (xlev-y0)) > 100*sqrt (eps)
 208   if verbose
 209     prn ("Error is too big : %8.3g\n", max (abs (xlev-y0)));
 210   end
 211   ok = 0;
 212 elseif verbose,  prn ("ok %i\n",cnt);
 213 end
 214
 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);
 218
 219 cnt++;
 220 if max (abs (xlev-y0)) > 100*sqrt (eps)
 221   if verbose
 222     prn ("Error is too big : %8.3g\n", max (abs (xlev-y0)));
 223   end
 224   ok = 0;
 225 elseif verbose,  prn ("ok %i\n",cnt);
 226 end
 227
 228 # ## See what 'backend' gives in that last case ########################
 229 # [method,ctl] = minimize ("ff",{x0,y0,1},"order",1,"backend");
 230
 231 # cnt++;
 232 # if ! strcmp (method,"bfgsmin")
 233 #   if verbose
 234 #     prn ("Wrong method '%s' != 'bfgsmin' was chosen\n", method);
 235 #   end
 236 #   ok = 0;
 237 # elseif verbose,  prn ("ok %i\n",cnt);
 238 # end
 239
 240 ## [xle2,vle2,nle2] = feval (method, "ff",{x0,y0,1}, ctl);
 241 [xle2,vle2,nle2] = minimize ("ff",{x0,y0,1},"order",1);
 242 cnt++;
 243 if max (abs (xlev-xle2)) > 100*eps
 244   if verbose
 245     prn ("Error is too big : %8.3g\n", max (abs (xlev-y0)));
 246   end
 247   ok = 0;
 248 elseif verbose,  prn ("ok %i\n",cnt);
 249 end
 250
 251
 252 if verbose && ok
 253   prn ( "All tests ok\n");
 254 end
 255