Add a useful package (from Source forge) for octave

[CreaPhase.git] / octave_packages / optim-1.2.0 / private / __dfdp__.m

diff --git a/octave_packages/optim-1.2.0/private/__dfdp__.m b/octave_packages/optim-1.2.0/private/__dfdp__.m
new file mode 100644 (file)
index 0000000..02a622e
--- /dev/null
+++ b/octave_packages/optim-1.2.0/private/__dfdp__.m
@@ -0,0 +1,156 @@
+%% Copyright (C) 1992-1994 Richard Shrager
+%% Copyright (C) 1992-1994 Arthur Jutan
+%% Copyright (C) 1992-1994 Ray Muzic
+%% Copyright (C) 2010, 2011 Olaf Till <olaf.till@uni-jena.de>
+%%
+%% This program is free software; you can redistribute it and/or modify it under
+%% the terms of the GNU General Public License as published by the Free Software
+%% Foundation; either version 3 of the License, or (at your option) any later
+%% version.
+%%
+%% This program is distributed in the hope that it will be useful, but WITHOUT
+%% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+%% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+%% details.
+%%
+%% You should have received a copy of the GNU General Public License along with
+%% this program; if not, see <http://www.gnu.org/licenses/>.
+
+function prt = __dfdp__ (p, func, hook)
+
+  %% Meant to be called by interfaces 'dfxpdp.m' and 'dfpdp.m', see there.
+
+  
+  if (nargin > 2 && isfield (hook, 'f'))
+    f = hook.f;
+  else
+    f = func (p);
+    f = f(:);
+  end
+
+  m = length (f);
+  n = length (p);
+
+  if (nargin > 2)
+
+    if (isfield (hook, 'fixed'))
+      fixed = hook.fixed;
+    else
+      fixed = false (n, 1);
+    end
+
+    if (isfield (hook, 'diffp'))
+      diffp = hook.diffp;
+    else
+      diffp = .001 * ones (n, 1);
+    end
+
+    if (isfield (hook, 'diff_onesided'))
+      diff_onesided = hook.diff_onesided;
+    else
+      diff_onesided = false (n, 1);
+    end
+
+    if (isfield (hook, 'lbound'))
+      lbound = hook.lbound;
+    else
+      lbound = - Inf (n, 1);
+    end
+
+    if (isfield (hook, 'ubound'))
+      ubound = hook.ubound;
+    else
+      ubound = Inf (n, 1);
+    end
+
+    if (isfield (hook, 'plabels'))
+      plabels = hook.plabels;
+    else
+      plabels = num2cell (num2cell ((1:n).'));
+    end
+
+  else
+    fixed = false (n, 1);
+    diff_onesided = fixed;
+    diffp = .001 * ones (n, 1);
+    lbound = - Inf (n, 1);
+    ubound = Inf (n, 1);
+    plabels = num2cell (num2cell ((1:n).'));
+  end    
+
+  prt = zeros (m, n); % initialise Jacobian to Zero
+  del = diffp .* p;
+  idxa = p == 0;
+  del(idxa) = diffp(idxa);
+  del(diff_onesided) = - del(diff_onesided); % keep course of
+                                % optimization of previous versions
+  absdel = abs (del);
+  idxd = ~(diff_onesided | fixed); % double sided interval
+  p1 = zeros (n, 1);
+  p2 = p1;
+  idxvs = false (n, 1);
+  idx1g2w = idxvs;
+  idx1le2w = idxvs;
+
+  %% p may be slightly out of bounds due to inaccuracy, or exactly at
+  %% the bound -> single sided interval
+  idxvl = p <= lbound;
+  idxvg = p >= ubound;
+  p1(idxvl) = min (p(idxvl, 1) + absdel(idxvl, 1), ubound(idxvl, 1));
+  idxd(idxvl) = false;
+  p1(idxvg) = max (p(idxvg, 1) - absdel(idxvg, 1), lbound(idxvg, 1));
+  idxd(idxvg) = false;
+  idxs = ~(fixed | idxd); % single sided interval
+
+  idxnv = ~(idxvl | idxvg); % current paramters within bounds
+  idxnvs = idxs & idxnv; % within bounds, single sided interval
+  idxnvd = idxd & idxnv; % within bounds, double sided interval
+  %% remaining single sided intervals
+  p1(idxnvs) = p(idxnvs) + del(idxnvs); % don't take absdel, this could
+                                % change course of optimization without
+                                % bounds with respect to previous
+                                % versions
+  %% remaining single sided intervals, violating a bound -> take largest
+  %% possible direction of single sided interval
+  idxvs(idxnvs) = p1(idxnvs, 1) < lbound(idxnvs, 1) | ...
+      p1(idxnvs, 1) > ubound(idxnvs, 1);
+  del1 = p(idxvs, 1) - lbound(idxvs, 1);
+  del2 = ubound(idxvs, 1) - p(idxvs, 1);
+  idx1g2 = del1 > del2;
+  idx1g2w(idxvs) = idx1g2;
+  idx1le2w(idxvs) = ~idx1g2;
+  p1(idx1g2w) = max (p(idx1g2w, 1) - absdel(idx1g2w, 1), ...
+                     lbound(idx1g2w, 1));
+  p1(idx1le2w) = min (p(idx1le2w, 1) + absdel(idx1le2w, 1), ...
+                      ubound(idx1le2w, 1));
+  %% double sided interval
+  p1(idxnvd) = min (p(idxnvd, 1) + absdel(idxnvd, 1), ...
+                    ubound(idxnvd, 1));
+  p2(idxnvd) = max (p(idxnvd, 1) - absdel(idxnvd, 1), ...
+                    lbound(idxnvd, 1));
+
+  del(idxs) = p1(idxs) - p(idxs);
+  del(idxd) = p1(idxd) - p2(idxd);
+
+  info.f = f;
+  info.parallel = false;
+
+  for j = 1:n
+    if (~fixed(j))
+      info.plabels = plabels(j, :);
+      ps = p;
+      ps(j) = p1(j);
+      if (idxs(j))
+        info.side = 0; % onesided interval
+        tp1 = func (ps, info);
+        prt(:, j) = (tp1(:) - f) / del(j);
+      else
+        info.side = 1; % centered interval, side 1
+        tp1 = func (ps, info);
+        ps(j) = p2(j);
+        info.side = 2; % centered interval, side 2
+        tp2 = func (ps, info);
+        prt(:, j) = (tp1(:) - tp2(:)) / del(j);
+      end
+    end
+  end