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>
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
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
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/>.
19 function prt = __dfdp__ (p, func, hook)
21 %% Meant to be called by interfaces 'dfxpdp.m' and 'dfpdp.m', see there.
24 if (nargin > 2 && isfield (hook, 'f'))
36 if (isfield (hook, 'fixed'))
42 if (isfield (hook, 'diffp'))
45 diffp = .001 * ones (n, 1);
48 if (isfield (hook, 'diff_onesided'))
49 diff_onesided = hook.diff_onesided;
51 diff_onesided = false (n, 1);
54 if (isfield (hook, 'lbound'))
57 lbound = - Inf (n, 1);
60 if (isfield (hook, 'ubound'))
66 if (isfield (hook, 'plabels'))
67 plabels = hook.plabels;
69 plabels = num2cell (num2cell ((1:n).'));
74 diff_onesided = fixed;
75 diffp = .001 * ones (n, 1);
76 lbound = - Inf (n, 1);
78 plabels = num2cell (num2cell ((1:n).'));
81 prt = zeros (m, n); % initialise Jacobian to Zero
84 del(idxa) = diffp(idxa);
85 del(diff_onesided) = - del(diff_onesided); % keep course of
86 % optimization of previous versions
88 idxd = ~(diff_onesided | fixed); % double sided interval
95 %% p may be slightly out of bounds due to inaccuracy, or exactly at
96 %% the bound -> single sided interval
99 p1(idxvl) = min (p(idxvl, 1) + absdel(idxvl, 1), ubound(idxvl, 1));
101 p1(idxvg) = max (p(idxvg, 1) - absdel(idxvg, 1), lbound(idxvg, 1));
103 idxs = ~(fixed | idxd); % single sided interval
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
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), ...
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), ...
129 p2(idxnvd) = max (p(idxnvd, 1) - absdel(idxnvd, 1), ...
132 del(idxs) = p1(idxs) - p(idxs);
133 del(idxd) = p1(idxd) - p2(idxd);
136 info.parallel = false;
140 info.plabels = plabels(j, :);
144 info.side = 0; % onesided interval
145 tp1 = func (ps, info);
146 prt(:, j) = (tp1(:) - f) / del(j);
148 info.side = 1; % centered interval, side 1
149 tp1 = func (ps, info);
151 info.side = 2; % centered interval, side 2
152 tp2 = func (ps, info);
153 prt(:, j) = (tp1(:) - tp2(:)) / del(j);