octave_packages/optim-1.2.0/test_d2_min_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 ## Test whether d2_min() functions correctly
  17 ##
  18 ## Gives a simple quadratic programming problem (function ff below).
  19 ##
  20 ## Sets a ok variable to 1 in case of success, 0 in case of failure
  21 ##
  22 ## If a variables "verbose" is set, then some comments are output.
  23
  24 1 ;
  25
  26 if ! exist ("verbose"), verbose = 0; end
  27
  28 if verbose
  29   printf ("\n   Testing d2_min () on a quadratic programming problem\n\n");
  30 end
  31
  32 P = 10+floor(30*rand(1)) ;      # Nparams
  33 R = P+floor(30*rand(1)) ;       # Nobses
  34 noise = 0 ;
  35 global obsmat ;
  36 obsmat = randn(R,P) ;
  37 global truep ;
  38 truep = randn(P,1) ;
  39 xinit = randn(P,1) ;
  40
  41 global obses ;
  42 obses = obsmat*truep ;
  43 if noise, obses = adnois(obses,noise); end
  44
  45
  46
  47 function v = ff(x)
  48   global obsmat;
  49   global obses;
  50   v = msq (obses - obsmat*x ) ;
  51 endfunction
  52
  53 function [v,dv,d2v] = d2ff(x) # Return pseudo-inverse
  54   global obsmat;
  55   global obses;
  56   er = -obses + obsmat*x ;
  57   dv = er'*obsmat ;
  58   v = msq(er ) ;
  59   d2v = pinv (obsmat'*obsmat ) ;
  60 endfunction
  61
  62 function [v,dv,d2v] = d2ff_2(x) # Return 2nd derivs, not pseudo-inv
  63   global obsmat;
  64   global obses;
  65   er = -obses + obsmat*x ;
  66   dv = er'*obsmat ;
  67   v = msq(er ) ;
  68   d2v = obsmat'*obsmat ;
  69 endfunction
  70
  71 ##       dt = mytic()
  72 ##
  73 ## Returns the cputime since last call to 'mytic'.
  74
  75 function dt = mytic()
  76    persistent last_mytic = 0 ;
  77    [t,u,s] = cputime() ;
  78    dt = t - last_mytic ;
  79    last_mytic = t ;
  80 endfunction
  81
  82 ## s = msq(x)                   - Mean squared value, ignoring nans
  83 ##
  84 ## s == mean(x(:).^2) , but ignores NaN's
  85
  86
  87 function s = msq(x)
  88 try
  89   s = mean(x(find(!isnan(x))).^2);
  90 catch
  91   s = nan;
  92 end
  93 endfunction
  94
  95 cnt = 1;
  96 ok = 1;
  97
  98 ctl = nan*zeros(1,5); ctl(5) = 1;
  99
 100 if verbose
 101   printf ("Going to call d2_min\n");
 102 end
 103 mytic() ;
 104 [xlev,vlev,nev] = d2_min ("ff","d2ff",xinit,ctl);
 105 tlev = mytic() ;
 106
 107 if verbose,
 108   printf("d2_min should find in one iteration + one more to check\n");
 109   printf(["d2_min : niter=%-4d  nev=%-4d  nobs=%-4d,nparams=%-4d\n",...
 110           "  time=%-8.3g errx=%-8.3g   minv=%-8.3g\n"],...
 111          nev([2,1]),R,P,tlev,max(abs(xlev-truep )),vlev);
 112 end
 113
 114
 115
 116 if nev(2) != 2,
 117   if verbose
 118       printf ("Too many iterations for this function\n");
 119   end
 120   ok = 0;
 121 else
 122   if verbose
 123       printf ("Ok: single iteration (%i)\n",cnt);
 124   end
 125 end
 126
 127 if max (abs(xlev-truep )) > sqrt (eps),
 128   if verbose
 129       printf ("Error is too big : %-8.3g\n", max (abs (xlev-truep)));
 130   end
 131   ok = 0;
 132 else
 133   if verbose
 134       printf ("Ok: single error amplitude (%i)\n",cnt);
 135   end
 136 end
 137
 138 cnt++;
 139
 140 if verbose
 141   printf ("Going to call d2_min() \n");
 142 end
 143 mytic() ;
 144 [xlev,vlev,nev] = d2_min("ff","d2ff_2",xinit) ;
 145 tlev = mytic() ;
 146
 147 if verbose,
 148   printf("d2_min should find in one iteration + one more to check\n");
 149   printf(["d2_min : niter=%-4d  nev=%-4d  nobs=%-4d,nparams=%-4d\n",...
 150           "  time=%-8.3g errx=%-8.3g   minv=%-8.3g\n"],...
 151          nev([2,1]),R,P,tlev,max(abs(xlev-truep )),vlev);
 152 end
 153
 154
 155 if nev(2) != 2,
 156   if verbose
 157       printf ("Too many iterations for this function\n");
 158   end
 159   ok = 0;
 160 else
 161   if verbose
 162       printf ("Ok: single iteration (%i)\n",cnt);
 163   end
 164 end
 165
 166 if max (abs(xlev-truep )) > sqrt (eps),
 167   if verbose
 168       printf ("Error is too big : %-8.3g\n", max (abs (xlev-truep)));
 169   end
 170   ok = 0;
 171 else
 172   if verbose
 173       printf ("Ok: single error amplitude (%i)\n",cnt);
 174   end
 175 end
 176
 177 if verbose
 178   if ok
 179     printf ("All tests ok\n");
 180   else
 181     printf ("Some tests failed\n");
 182   end
 183 end
 184