octave_packages/optim-1.2.0/private/__dfdp__.m

   1 %% Copyright (C) 1992-1994 Richard Shrager
   2 %% Copyright (C) 1992-1994 Arthur Jutan
   3 %% Copyright (C) 1992-1994 Ray Muzic
   4 %% Copyright (C) 2010, 2011 Olaf Till <olaf.till@uni-jena.de>
   5 %%
   6 %% This program is free software; you can redistribute it and/or modify it under
   7 %% the terms of the GNU General Public License as published by the Free Software
   8 %% Foundation; either version 3 of the License, or (at your option) any later
   9 %% version.
  10 %%
  11 %% This program is distributed in the hope that it will be useful, but WITHOUT
  12 %% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  13 %% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  14 %% details.
  15 %%
  16 %% You should have received a copy of the GNU General Public License along with
  17 %% this program; if not, see <http://www.gnu.org/licenses/>.
  18
  19 function prt = __dfdp__ (p, func, hook)
  20
  21   %% Meant to be called by interfaces 'dfxpdp.m' and 'dfpdp.m', see there.
  22
  23
  24   if (nargin > 2 && isfield (hook, 'f'))
  25     f = hook.f;
  26   else
  27     f = func (p);
  28     f = f(:);
  29   end
  30
  31   m = length (f);
  32   n = length (p);
  33
  34   if (nargin > 2)
  35
  36     if (isfield (hook, 'fixed'))
  37       fixed = hook.fixed;
  38     else
  39       fixed = false (n, 1);
  40     end
  41
  42     if (isfield (hook, 'diffp'))
  43       diffp = hook.diffp;
  44     else
  45       diffp = .001 * ones (n, 1);
  46     end
  47
  48     if (isfield (hook, 'diff_onesided'))
  49       diff_onesided = hook.diff_onesided;
  50     else
  51       diff_onesided = false (n, 1);
  52     end
  53
  54     if (isfield (hook, 'lbound'))
  55       lbound = hook.lbound;
  56     else
  57       lbound = - Inf (n, 1);
  58     end
  59
  60     if (isfield (hook, 'ubound'))
  61       ubound = hook.ubound;
  62     else
  63       ubound = Inf (n, 1);
  64     end
  65
  66     if (isfield (hook, 'plabels'))
  67       plabels = hook.plabels;
  68     else
  69       plabels = num2cell (num2cell ((1:n).'));
  70     end
  71
  72   else
  73     fixed = false (n, 1);
  74     diff_onesided = fixed;
  75     diffp = .001 * ones (n, 1);
  76     lbound = - Inf (n, 1);
  77     ubound = Inf (n, 1);
  78     plabels = num2cell (num2cell ((1:n).'));
  79   end
  80
  81   prt = zeros (m, n); % initialise Jacobian to Zero
  82   del = diffp .* p;
  83   idxa = p == 0;
  84   del(idxa) = diffp(idxa);
  85   del(diff_onesided) = - del(diff_onesided); % keep course of
  86                                 % optimization of previous versions
  87   absdel = abs (del);
  88   idxd = ~(diff_onesided | fixed); % double sided interval
  89   p1 = zeros (n, 1);
  90   p2 = p1;
  91   idxvs = false (n, 1);
  92   idx1g2w = idxvs;
  93   idx1le2w = idxvs;
  94
  95   %% p may be slightly out of bounds due to inaccuracy, or exactly at
  96   %% the bound -> single sided interval
  97   idxvl = p <= lbound;
  98   idxvg = p >= ubound;
  99   p1(idxvl) = min (p(idxvl, 1) + absdel(idxvl, 1), ubound(idxvl, 1));
 100   idxd(idxvl) = false;
 101   p1(idxvg) = max (p(idxvg, 1) - absdel(idxvg, 1), lbound(idxvg, 1));
 102   idxd(idxvg) = false;
 103   idxs = ~(fixed | idxd); % single sided interval
 104
 105   idxnv = ~(idxvl | idxvg); % current paramters within bounds
 106   idxnvs = idxs & idxnv; % within bounds, single sided interval
 107   idxnvd = idxd & idxnv; % within bounds, double sided interval
 108   %% remaining single sided intervals
 109   p1(idxnvs) = p(idxnvs) + del(idxnvs); % don't take absdel, this could
 110                                 % change course of optimization without
 111                                 % bounds with respect to previous
 112                                 % versions
 113   %% remaining single sided intervals, violating a bound -> take largest
 114   %% possible direction of single sided interval
 115   idxvs(idxnvs) = p1(idxnvs, 1) < lbound(idxnvs, 1) | ...
 116       p1(idxnvs, 1) > ubound(idxnvs, 1);
 117   del1 = p(idxvs, 1) - lbound(idxvs, 1);
 118   del2 = ubound(idxvs, 1) - p(idxvs, 1);
 119   idx1g2 = del1 > del2;
 120   idx1g2w(idxvs) = idx1g2;
 121   idx1le2w(idxvs) = ~idx1g2;
 122   p1(idx1g2w) = max (p(idx1g2w, 1) - absdel(idx1g2w, 1), ...
 123                      lbound(idx1g2w, 1));
 124   p1(idx1le2w) = min (p(idx1le2w, 1) + absdel(idx1le2w, 1), ...
 125                       ubound(idx1le2w, 1));
 126   %% double sided interval
 127   p1(idxnvd) = min (p(idxnvd, 1) + absdel(idxnvd, 1), ...
 128                     ubound(idxnvd, 1));
 129   p2(idxnvd) = max (p(idxnvd, 1) - absdel(idxnvd, 1), ...
 130                     lbound(idxnvd, 1));
 131
 132   del(idxs) = p1(idxs) - p(idxs);
 133   del(idxd) = p1(idxd) - p2(idxd);
 134
 135   info.f = f;
 136   info.parallel = false;
 137
 138   for j = 1:n
 139     if (~fixed(j))
 140       info.plabels = plabels(j, :);
 141       ps = p;
 142       ps(j) = p1(j);
 143       if (idxs(j))
 144         info.side = 0; % onesided interval
 145         tp1 = func (ps, info);
 146         prt(:, j) = (tp1(:) - f) / del(j);
 147       else
 148         info.side = 1; % centered interval, side 1
 149         tp1 = func (ps, info);
 150         ps(j) = p2(j);
 151         info.side = 2; % centered interval, side 2
 152         tp2 = func (ps, info);
 153         prt(:, j) = (tp1(:) - tp2(:)) / del(j);
 154       end
 155     end
 156   end