Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[10510] trunk/octave-forge/main/optim
|
From: <i7...@us...> - 2012-05-24 05:50:50
|
Revision: 10510
http://octave.svn.sourceforge.net/octave/?rev=10510&view=rev
Author: i7tiol
Date: 2012-05-24 05:50:43 +0000 (Thu, 24 May 2012)
Log Message:
-----------
Untabify and some cleanup.
Modified Paths:
--------------
trunk/octave-forge/main/optim/DESCRIPTION
trunk/octave-forge/main/optim/PKG_ADD
trunk/octave-forge/main/optim/doc/development/interfaces.txt
trunk/octave-forge/main/optim/inst/cpiv_bard.m
trunk/octave-forge/main/optim/inst/dcdp.m
trunk/octave-forge/main/optim/inst/dfdp.m
trunk/octave-forge/main/optim/inst/dfpdp.m
trunk/octave-forge/main/optim/inst/dfxpdp.m
trunk/octave-forge/main/optim/inst/gjp.m
trunk/octave-forge/main/optim/inst/leasqr.m
trunk/octave-forge/main/optim/inst/nonlin_residmin.m
trunk/octave-forge/main/optim/inst/optim_problems.m
trunk/octave-forge/main/optim/inst/private/__dfdp__.m
trunk/octave-forge/main/optim/inst/private/__lm_svd__.m
trunk/octave-forge/main/optim/inst/private/__nonlin_residmin__.m
Modified: trunk/octave-forge/main/optim/DESCRIPTION
===================================================================
--- trunk/octave-forge/main/optim/DESCRIPTION 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/DESCRIPTION 2012-05-24 05:50:43 UTC (rev 10510)
@@ -1,11 +1,11 @@
Name: Optim
Version: 1.1.0
-Date: 2012-XX-XX
+Date: 2012-05-24
Author: various authors
Maintainer: Octave-Forge community <oct...@li...>
Title: Optimization.
Description: Non-linear optimization toolkit.
-Depends: octave (>= 2.9.7), miscellaneous (>= 1.0.10), struct (>= 1.0.9)
+Depends: octave (>= 2.9.7), miscellaneous (>= 1.0.10), struct (>= 1.0.10)
Autoload: no
License: GFDL, GPLv3+, modified BSD, public domain
Url: http://octave.sf.net
Modified: trunk/octave-forge/main/optim/PKG_ADD
===================================================================
--- trunk/octave-forge/main/optim/PKG_ADD 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/PKG_ADD 2012-05-24 05:50:43 UTC (rev 10510)
@@ -3,4 +3,5 @@
## sometime in 3.3.54+, if I remember right
__all_opts__ ("nonlin_residmin");
__all_opts__ ("residmin_stat");
+ ## __all_opts__ ("nonlin_min");
endif
Modified: trunk/octave-forge/main/optim/doc/development/interfaces.txt
===================================================================
--- trunk/octave-forge/main/optim/doc/development/interfaces.txt 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/doc/development/interfaces.txt 2012-05-24 05:50:43 UTC (rev 10510)
@@ -80,6 +80,9 @@
'vc' (required): vector (possibly empty) of the function "mc.' *
parameters + vc" of linear constraints,
+If bounds have been specified, they are contained in 'mc' and 'vc'
+_before_ all other linear constraints.
+
'n_gencstr' (required): number of general constraints (except the
linear constraints given by 'mc' and 'vc',
Modified: trunk/octave-forge/main/optim/inst/cpiv_bard.m
===================================================================
--- trunk/octave-forge/main/optim/inst/cpiv_bard.m 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/inst/cpiv_bard.m 2012-05-24 05:50:43 UTC (rev 10510)
@@ -1,4 +1,4 @@
-%% Copyright (C) 2010, 2011 Olaf Till <ola...@un...>
+%% Copyright (C) 2010, 2011 Olaf Till <i7...@t-...>
%%
%% 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
Modified: trunk/octave-forge/main/optim/inst/dcdp.m
===================================================================
--- trunk/octave-forge/main/optim/inst/dcdp.m 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/inst/dcdp.m 2012-05-24 05:50:43 UTC (rev 10510)
@@ -1,4 +1,4 @@
-%% Copyright (C) 2010, 2011 Olaf Till <ola...@un...>
+%% Copyright (C) 2010, 2011 Olaf Till <i7...@t-...>
%%
%% 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
@@ -40,4 +40,3 @@
hook.diff_onesided = dp < 0;
prt = __dfdp__ (p, func, hook);
-end
Modified: trunk/octave-forge/main/optim/inst/dfdp.m
===================================================================
--- trunk/octave-forge/main/optim/inst/dfdp.m 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/inst/dfdp.m 2012-05-24 05:50:43 UTC (rev 10510)
@@ -1,7 +1,7 @@
%% Copyright (C) 1992-1994 Richard Shrager
%% Copyright (C) 1992-1994 Arthur Jutan
%% Copyright (C) 1992-1994 Ray Muzic
-%% Copyright (C) 2010, 2011 Olaf Till <ola...@un...>
+%% Copyright (C) 2010, 2011 Olaf Till <i7...@t-...>
%%
%% 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
Modified: trunk/octave-forge/main/optim/inst/dfpdp.m
===================================================================
--- trunk/octave-forge/main/optim/inst/dfpdp.m 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/inst/dfpdp.m 2012-05-24 05:50:43 UTC (rev 10510)
@@ -48,5 +48,3 @@
end
ret = __dfdp__ (varargin{:});
-
-end
Modified: trunk/octave-forge/main/optim/inst/dfxpdp.m
===================================================================
--- trunk/octave-forge/main/optim/inst/dfxpdp.m 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/inst/dfxpdp.m 2012-05-24 05:50:43 UTC (rev 10510)
@@ -1,4 +1,4 @@
-%% Copyright (C) 2010, 2011 Olaf Till <ola...@un...>
+%% Copyright (C) 2010, 2011 Olaf Till <i7...@t-...>
%%
%% 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
@@ -51,5 +51,3 @@
end
ret = __dfdp__ (varargin{2:end});
-
-end
Modified: trunk/octave-forge/main/optim/inst/gjp.m
===================================================================
--- trunk/octave-forge/main/optim/inst/gjp.m 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/inst/gjp.m 2012-05-24 05:50:43 UTC (rev 10510)
@@ -1,4 +1,4 @@
-%% Copyright (C) 2010, 2011 Olaf Till <ola...@un...>
+%% Copyright (C) 2010, 2011 Olaf Till <i7...@t-...>
%%
%% 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
@@ -50,4 +50,3 @@
m([1:k-1, k+1:end], l) * m(k, [1:l-1, l+1:end]);
m([1:k-1, k+1:end], l) = - m([1:k-1, k+1:end], l) / p; % pivot column
m(k, l) = 1 / p;
-end
Modified: trunk/octave-forge/main/optim/inst/leasqr.m
===================================================================
--- trunk/octave-forge/main/optim/inst/leasqr.m 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/inst/leasqr.m 2012-05-24 05:50:43 UTC (rev 10510)
@@ -1,7 +1,7 @@
%% Copyright (C) 1992-1994 Richard Shrager
%% Copyright (C) 1992-1994 Arthur Jutan <ju...@ch...>
%% Copyright (C) 1992-1994 Ray Muzic <rf...@ds...>
-%% Copyright (C) 2010, 2011 Olaf Till <ola...@un...>
+%% Copyright (C) 2010, 2011 Olaf Till <i7...@t-...>
%%
%% 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
Modified: trunk/octave-forge/main/optim/inst/nonlin_residmin.m
===================================================================
--- trunk/octave-forge/main/optim/inst/nonlin_residmin.m 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/inst/nonlin_residmin.m 2012-05-24 05:50:43 UTC (rev 10510)
@@ -90,7 +90,8 @@
## intervals should be used by jacobian functions performing finite
## differencing. Default: @code{false (size (parameters))}.
##
-## @code{complex_step_derivative}, @code{complex_step_derivative_inequc},
+## @code{complex_step_derivative_f},
+## @code{complex_step_derivative_inequc},
## @code{complex_step_derivative_equc}: logical scalars, default: false.
## Estimate Jacobian of model function, general inequality constraints,
## and general equality constraints, respectively, with complex step
@@ -176,7 +177,7 @@
##
## @code{plot_cmd}: Function enabling backend to plot results or
## intermediate results. Will be called with current computed
-## residualse. Default: plot nothing.
+## residuals. Default: plot nothing.
##
## @code{debug}: Logical scalar, default: @code{false}. Will be passed
## to the backend, which might print debugging information if true.
Modified: trunk/octave-forge/main/optim/inst/optim_problems.m
===================================================================
--- trunk/octave-forge/main/optim/inst/optim_problems.m 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/inst/optim_problems.m 2012-05-24 05:50:43 UTC (rev 10510)
@@ -1,6 +1,6 @@
%% Copyright (C) 2007 Paul Kienzle (sort-based lookup in ODE solver)
-%% Copyright (C) 2009 Thomas Treichl <tho...@gm...>
-%% Copyright (C) 2010 Olaf Till <ola...@un...>
+%% Copyright (C) 2009 Thomas Treichl <tho...@gm...> (ode23 code)
+%% Copyright (C) 2010 Olaf Till <i7...@t-...>
%%
%% 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
Modified: trunk/octave-forge/main/optim/inst/private/__dfdp__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/private/__dfdp__.m 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/inst/private/__dfdp__.m 2012-05-24 05:50:43 UTC (rev 10510)
@@ -83,7 +83,7 @@
idxa = p == 0;
del(idxa) = diffp(idxa);
del(diff_onesided) = - del(diff_onesided); % keep course of
- % optimization of previous versions
+ % optimization of previous versions
absdel = abs (del);
idxd = ~(diff_onesided | fixed); % double sided interval
p1 = zeros (n, 1);
@@ -107,9 +107,9 @@
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
+ % 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) | ...
@@ -120,14 +120,14 @@
idx1g2w(idxvs) = idx1g2;
idx1le2w(idxvs) = ~idx1g2;
p1(idx1g2w) = max (p(idx1g2w, 1) - absdel(idx1g2w, 1), ...
- lbound(idx1g2w, 1));
+ lbound(idx1g2w, 1));
p1(idx1le2w) = min (p(idx1le2w, 1) + absdel(idx1le2w, 1), ...
- ubound(idx1le2w, 1));
+ ubound(idx1le2w, 1));
%% double sided interval
p1(idxnvd) = min (p(idxnvd, 1) + absdel(idxnvd, 1), ...
- ubound(idxnvd, 1));
+ ubound(idxnvd, 1));
p2(idxnvd) = max (p(idxnvd, 1) - absdel(idxnvd, 1), ...
- lbound(idxnvd, 1));
+ lbound(idxnvd, 1));
del(idxs) = p1(idxs) - p(idxs);
del(idxd) = p1(idxd) - p2(idxd);
@@ -141,16 +141,16 @@
ps = p;
ps(j) = p1(j);
if (idxs(j))
- info.side = 0; % onesided interval
- tp1 = func (ps, info);
- prt(:, j) = (tp1(:) - f) / del(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);
+ 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
Modified: trunk/octave-forge/main/optim/inst/private/__lm_svd__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/private/__lm_svd__.m 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/inst/private/__lm_svd__.m 2012-05-24 05:50:43 UTC (rev 10510)
@@ -37,7 +37,7 @@
df_cstr = hook.df_cstr; % function of derivatives of all constraints
n_gencstr = hook.n_gencstr; % number of non-linear constraints
eq_idx = hook.eq_idx; % logical index of equality constraints in all
- % constraints
+ % constraints
lbound = hook.lbound; % bounds, subset of linear inequality
ubound = hook.ubound; % constraints in mc and vc
@@ -86,16 +86,16 @@
wtl = wt(:);
nz = 20 * eps; % This is arbitrary. Constraint function will be
- % regarded as <= zero if less than nz.
+ % regarded as <= zero if less than nz.
%% backend-specific checking of options and constraints
if (have_constraints_except_bounds)
if (any (pin_cstr.inequ.lin_except_bounds < 0) || ...
- (n_gencstr > 0 && any (pin_cstr.inequ.gen < 0)))
+ (n_gencstr > 0 && any (pin_cstr.inequ.gen < 0)))
warning ('initial parameters violate inequality constraints');
end
if (any (abs (pin_cstr.equ.lin) >= nz) || ...
- (n_gencstr > 0 && any (abs (pin_cstr.equ.gen) >= nz)))
+ (n_gencstr > 0 && any (abs (pin_cstr.equ.gen) >= nz)))
warning ('initial parameters violate equality constraints');
end
end
@@ -128,7 +128,7 @@
%% fill constant fields of hook for derivative-functions; some fields
%% may be backend-specific
dfdp_hook.fixed = fixed; % this may be handled by the frontend, but
- % the backend still may add to it
+ % the backend still may add to it
%% set up for iterations
%%
@@ -156,18 +156,18 @@
v_cstr = f_cstr (p, ac_idx);
%% index of active constraints
c_act = v_cstr < nz | eq_idx; # equality constraints might be
- # violated at start
+ # violated at start
if (any (c_act))
if (n_gencstr > 0)
- dct = df_cstr (p, ac_idx, ...
- setfield (dfdp_hook, 'f', v_cstr));
- dct(:, fixed) = 0; % for user supplied dfdp; necessary?
- dc = dct.';
- dcat = dct(c_act, :);
+ dct = df_cstr (p, ac_idx, ...
+ setfield (dfdp_hook, 'f', v_cstr));
+ dct(:, fixed) = 0; % for user supplied dfdp; necessary?
+ dc = dct.';
+ dcat = dct(c_act, :);
else
- dcat = df_cstr (p, c_act, ...
- setfield (dfdp_hook, 'f', v_cstr));
- dcat(:, fixed) = 0; % for user supplied dfdp; necessary?
+ dcat = df_cstr (p, c_act, ...
+ setfield (dfdp_hook, 'f', v_cstr));
+ dcat(:, fixed) = 0; % for user supplied dfdp; necessary?
end
dca = dcat.';
end
@@ -207,46 +207,46 @@
%%% ser = 1 ./ sqrt((s.*s)+epsL);
se = sqrt ((s.*s) + epsL);
if (new_s)
- %% for testing
- ser = (1 / (1 + epsL^2)) * (1 ./ se + epsL * s);
+ %% for testing
+ ser = (1 / (1 + epsL^2)) * (1 ./ se + epsL * s);
else
- ser = 1 ./ se;
+ ser = 1 ./ se;
end
tp1 = (v * (g .* ser)) .* nrm;
if (any (c_act))
- deb_printf (testing, 'constraints are active:\n');
- deb_printf (testing, '%i\n', c_act);
- %% calculate chg by 'quadratic programming'
- nrme= diag (nrm);
- ser2 = diag (ser .* ser);
- mfc1 = nrme * v * ser2 * v.' * nrme;
- tp2 = mfc1 * dca;
- a_eq_idx = eq_idx(c_act);
- [lb, bidx, ridx, tbl] = cpiv (dcat * tp1, dcat * tp2, a_eq_idx);
- chg = tp1 + tp2(:, bidx) * lb; % if a parameter is 'fixed',
- % the respective component of chg should
- % be zero too, even here (with active
- % constraints)
- deb_printf (testing, 'change:\n');
- deb_printf (testing, '%e\n', chg);
- deb_printf (testing, '\n');
- %% indices for different types of constraints
- c_inact = ~c_act; % inactive constraints
- c_binding = nc_idx;
- c_binding(c_act) = bidx; % constraints selected binding
- c_unbinding = nc_idx;
- c_unbinding(c_act) = ridx; % constraints unselected binding
- c_nonbinding = c_act & ~(c_binding | c_unbinding); % constraints
- % selected non-binding
+ deb_printf (testing, 'constraints are active:\n');
+ deb_printf (testing, '%i\n', c_act);
+ %% calculate chg by 'quadratic programming'
+ nrme= diag (nrm);
+ ser2 = diag (ser .* ser);
+ mfc1 = nrme * v * ser2 * v.' * nrme;
+ tp2 = mfc1 * dca;
+ a_eq_idx = eq_idx(c_act);
+ [lb, bidx, ridx, tbl] = cpiv (dcat * tp1, dcat * tp2, a_eq_idx);
+ chg = tp1 + tp2(:, bidx) * lb; % if a parameter is 'fixed',
+ % the respective component of chg should
+ % be zero too, even here (with active
+ % constraints)
+ deb_printf (testing, 'change:\n');
+ deb_printf (testing, '%e\n', chg);
+ deb_printf (testing, '\n');
+ %% indices for different types of constraints
+ c_inact = ~c_act; % inactive constraints
+ c_binding = nc_idx;
+ c_binding(c_act) = bidx; % constraints selected binding
+ c_unbinding = nc_idx;
+ c_unbinding(c_act) = ridx; % constraints unselected binding
+ c_nonbinding = c_act & ~(c_binding | c_unbinding); % constraints
+ % selected non-binding
else
- %% chg is the Levenberg/Marquardt step
- chg = tp1;
- %% indices for different types of constraints
- c_inact = ac_idx; % inactive constraints consist of all
- % constraints
- c_binding = nc_idx;
- c_unbinding = nc_idx;
- c_nonbinding = nc_idx;
+ %% chg is the Levenberg/Marquardt step
+ chg = tp1;
+ %% indices for different types of constraints
+ c_inact = ac_idx; % inactive constraints consist of all
+ % constraints
+ c_binding = nc_idx;
+ c_unbinding = nc_idx;
+ c_nonbinding = nc_idx;
end
%% apply constraints to step width (since this is a
%% Levenberg/Marquardt algorithm, no line-search is performed
@@ -258,200 +258,200 @@
hstep = mcit * chg;
idx = hstep < 0;
if (any (idx))
- k = min (1, min (- (vci(idx) + mcit(idx, :) * pprev) ./ ...
- hstep(idx)));
+ k = min (1, min (- (vci(idx) + mcit(idx, :) * pprev) ./ ...
+ hstep(idx)));
end
if (k < 1)
- deb_printf (testing, 'stepwidth: linear constraints\n');
+ deb_printf (testing, 'stepwidth: linear constraints\n');
end
if (n_gencstr > 0)
- c_tp = gc_idx & (c_nonbinding | c_inact);
- if (any (c_tp) && any (f_cstr (pprev + k * chg, c_tp) < 0))
- [k, fval, info] = ...
- fzero (@ (x) min (cat (1, ...
- f_cstr (pprev + x * chg, c_tp), ...
- k - x, ...
- ifelse (x < 0, -Inf, Inf))), ...
- 0);
- if (info ~= 1 || abs (fval) >= nz)
- error ('could not find stepwidth to satisfy inactive and non-binding general inequality constraints');
- end
- deb_printf (testing, 'general constraints limit stepwidth\n');
- end
+ c_tp = gc_idx & (c_nonbinding | c_inact);
+ if (any (c_tp) && any (f_cstr (pprev + k * chg, c_tp) < 0))
+ [k, fval, info] = ...
+ fzero (@ (x) min (cat (1, ...
+ f_cstr (pprev + x * chg, c_tp), ...
+ k - x, ...
+ ifelse (x < 0, -Inf, Inf))), ...
+ 0);
+ if (info ~= 1 || abs (fval) >= nz)
+ error ('could not find stepwidth to satisfy inactive and non-binding general inequality constraints');
+ end
+ deb_printf (testing, 'general constraints limit stepwidth\n');
+ end
end
chg = k * chg;
if (any (gc_idx & c_binding)) % none selected binding =>
- % none unselected binding
- deb_printf (testing, 'general binding constraints must be regained:\n');
- %% regain binding constraints and one of the possibly active
- %% previously inactive or non-binding constraints
- ptp1 = pprev + chg;
+ % none unselected binding
+ deb_printf (testing, 'general binding constraints must be regained:\n');
+ %% regain binding constraints and one of the possibly active
+ %% previously inactive or non-binding constraints
+ ptp1 = pprev + chg;
- tp = true;
- nt_nosuc = true;
- lim = 20;
- while (nt_nosuc && lim >= 0)
- deb_printf (testing, 'starting from new value of p in regaining:\n');
- deb_printf (testing, '%e\n', ptp1);
- %% we keep d_p.' * inv (mfc1) * d_p minimal in each step of
- %% the inner loop; this is both sensible (this metric
- %% considers a guess of curvature of sum of squared residuals)
- %% and convenient (we have useful matrices available for it)
- c_tp0 = c_inact | c_nonbinding;
- c_tp1 = c_inact | (gc_idx & c_nonbinding);
- btbl = tbl(bidx, bidx);
- c_tp2 = c_binding;
- if (any (tp)) % if none before, does not get true again
- tp = f_cstr (ptp1, c_tp1) < nz;
- if (any (tp)) % could be less clumsy, but ml-compatibility..
- %% keep only the first true entry in tp
- tp(tp) = logical (cat (1, 1, zeros (sum (tp) - 1, 1)));
- %% supplement binding index with one (the first) getting
- %% binding in c_tp1
- c_tp2(c_tp1) = tp;
- %% gradient of this added constraint
- caddt = dct(c_tp2 & ~c_binding, :);
- cadd = caddt.';
- C = dct(c_binding, :) * mfc1 * cadd;
- Ct = C.';
- G = [btbl, btbl * C; ...
- -Ct * btbl, caddt * mfc1 * cadd - Ct * btbl * C];
- btbl = gjp (G, size (G, 1));
- end
- end
- dcbt = dct(c_tp2, :);
- mfc = - mfc1 * dcbt.' * btbl;
- deb_printf (testing, 'constraints to regain:\n');
- deb_printf (testing, '%i\n', c_tp2);
+ tp = true;
+ nt_nosuc = true;
+ lim = 20;
+ while (nt_nosuc && lim >= 0)
+ deb_printf (testing, 'starting from new value of p in regaining:\n');
+ deb_printf (testing, '%e\n', ptp1);
+ %% we keep d_p.' * inv (mfc1) * d_p minimal in each step of
+ %% the inner loop; this is both sensible (this metric
+ %% considers a guess of curvature of sum of squared residuals)
+ %% and convenient (we have useful matrices available for it)
+ c_tp0 = c_inact | c_nonbinding;
+ c_tp1 = c_inact | (gc_idx & c_nonbinding);
+ btbl = tbl(bidx, bidx);
+ c_tp2 = c_binding;
+ if (any (tp)) % if none before, does not get true again
+ tp = f_cstr (ptp1, c_tp1) < nz;
+ if (any (tp)) % could be less clumsy, but ml-compatibility..
+ %% keep only the first true entry in tp
+ tp(tp) = logical (cat (1, 1, zeros (sum (tp) - 1, 1)));
+ %% supplement binding index with one (the first) getting
+ %% binding in c_tp1
+ c_tp2(c_tp1) = tp;
+ %% gradient of this added constraint
+ caddt = dct(c_tp2 & ~c_binding, :);
+ cadd = caddt.';
+ C = dct(c_binding, :) * mfc1 * cadd;
+ Ct = C.';
+ G = [btbl, btbl * C; ...
+ -Ct * btbl, caddt * mfc1 * cadd - Ct * btbl * C];
+ btbl = gjp (G, size (G, 1));
+ end
+ end
+ dcbt = dct(c_tp2, :);
+ mfc = - mfc1 * dcbt.' * btbl;
+ deb_printf (testing, 'constraints to regain:\n');
+ deb_printf (testing, '%i\n', c_tp2);
- ptp2 = ptp1;
- nt_niter_start = 100;
- nt_niter = nt_niter_start;
- while (nt_nosuc && nt_niter >= 0)
- hv = f_cstr (ptp2, c_tp2);
- if (all (abs (hv) < nz))
- nt_nosuc = false;
- chg = ptp2 - pprev;
- else
- ptp2 = ptp2 + mfc * hv; % step should be zero for each
- % component for which the parameter is
- % 'fixed'
- end
- nt_niter = nt_niter - 1;
- end
- deb_printf (testing, 'constraints after regaining:\n');
- deb_printf (testing, '%e\n', hv);
- if (nt_nosuc || ...
- any (abs (chg) > abs (pprev .* maxstep)) || ...
- any (f_cstr (ptp2, c_tp0) < -nz))
- if (nt_nosuc)
- deb_printf (testing, 'regaining did not converge\n');
- else
- deb_printf (testing, 'regaining violated type 3 and 4\n');
- end
- nt_nosuc = true;
- ptp1 = (pprev + ptp1) / 2;
- end
- if (~nt_nosuc)
- tp = f_cstr (ptp2, c_unbinding);
- if (any (tp) < 0) % again ml-compatibility clumsyness..
- [discarded, id] = min(tp);
- tid = find (ridx);
- id = tid(id); % index within active constraints
- unsuccessful_exchange = false;
- if (abs (tbl(id, id)) < nz) % Bard: not absolute value
- %% exchange this unselected binding constraint against a
- %% binding constraint, but not against an equality
- %% constraint
- tbidx = bidx & ~a_eq_idx;
- if (~any (tbidx))
- unsuccessful_exchange = true;
- else
- [discarded, idm] = max (abs (tbl(tbidx, id)));
- tid = find (tbidx);
- idm = tid(idm); % -> index within active constraints
- tbl = gjp (tbl, idm);
- bidx(idm) = false;
- ridx(idm) = true;
- end
- end
- if (unsuccessful_exchange)
- %% It probably doesn't look good now; this desperate
- %% last attempt is not in the original algortithm, since
- %% that didn't account for equality constraints.
- ptp1 = (pprev + ptp1) / 2;
- else
- tbl = gjp (tbl, id);
- bidx(id) = true;
- ridx(id) = false;
- c_binding = nc_idx;
- c_binding(c_act) = bidx;
- c_unbinding = nc_idx;
- c_unbinding(c_act) = ridx;
- end
- nt_nosuc = true;
- deb_printf (testing, 'regaining violated type 2\n');
- end
- end
- if (~nt_nosuc)
- deb_printf (testing, 'regaining successful, converged with %i iterations:\n', ...
- nt_niter_start - nt_niter);
- deb_printf (testing, '%e\n', ptp2);
- end
- lim = lim - 1;
- end
- if (nt_nosuc)
- error ('could not regain binding constraints');
- end
+ ptp2 = ptp1;
+ nt_niter_start = 100;
+ nt_niter = nt_niter_start;
+ while (nt_nosuc && nt_niter >= 0)
+ hv = f_cstr (ptp2, c_tp2);
+ if (all (abs (hv) < nz))
+ nt_nosuc = false;
+ chg = ptp2 - pprev;
+ else
+ ptp2 = ptp2 + mfc * hv; % step should be zero for each
+ % component for which the parameter is
+ % 'fixed'
+ end
+ nt_niter = nt_niter - 1;
+ end
+ deb_printf (testing, 'constraints after regaining:\n');
+ deb_printf (testing, '%e\n', hv);
+ if (nt_nosuc || ...
+ any (abs (chg) > abs (pprev .* maxstep)) || ...
+ any (f_cstr (ptp2, c_tp0) < -nz))
+ if (nt_nosuc)
+ deb_printf (testing, 'regaining did not converge\n');
+ else
+ deb_printf (testing, 'regaining violated type 3 and 4\n');
+ end
+ nt_nosuc = true;
+ ptp1 = (pprev + ptp1) / 2;
+ end
+ if (~nt_nosuc)
+ tp = f_cstr (ptp2, c_unbinding);
+ if (any (tp) < 0) % again ml-compatibility clumsyness..
+ [discarded, id] = min(tp);
+ tid = find (ridx);
+ id = tid(id); % index within active constraints
+ unsuccessful_exchange = false;
+ if (abs (tbl(id, id)) < nz) % Bard: not absolute value
+ %% exchange this unselected binding constraint against a
+ %% binding constraint, but not against an equality
+ %% constraint
+ tbidx = bidx & ~a_eq_idx;
+ if (~any (tbidx))
+ unsuccessful_exchange = true;
+ else
+ [discarded, idm] = max (abs (tbl(tbidx, id)));
+ tid = find (tbidx);
+ idm = tid(idm); % -> index within active constraints
+ tbl = gjp (tbl, idm);
+ bidx(idm) = false;
+ ridx(idm) = true;
+ end
+ end
+ if (unsuccessful_exchange)
+ %% It probably doesn't look good now; this desperate
+ %% last attempt is not in the original algortithm, since
+ %% that didn't account for equality constraints.
+ ptp1 = (pprev + ptp1) / 2;
+ else
+ tbl = gjp (tbl, id);
+ bidx(id) = true;
+ ridx(id) = false;
+ c_binding = nc_idx;
+ c_binding(c_act) = bidx;
+ c_unbinding = nc_idx;
+ c_unbinding(c_act) = ridx;
+ end
+ nt_nosuc = true;
+ deb_printf (testing, 'regaining violated type 2\n');
+ end
+ end
+ if (~nt_nosuc)
+ deb_printf (testing, 'regaining successful, converged with %i iterations:\n', ...
+ nt_niter_start - nt_niter);
+ deb_printf (testing, '%e\n', ptp2);
+ end
+ lim = lim - 1;
+ end
+ if (nt_nosuc)
+ error ('could not regain binding constraints');
+ end
else
- %% check the maximal stepwidth and apply as necessary
- ochg=chg;
- idx = ~isinf(maxstep);
- limit = abs(maxstep(idx).*pprev(idx));
- chg(idx) = min(max(chg(idx),-limit),limit);
- if (verbose && any(ochg ~= chg))
- disp(['Change in parameter(s): ', ...
- sprintf('%d ',find(ochg ~= chg)), 'maximal fractional stepwidth enforced']);
- end
+ %% check the maximal stepwidth and apply as necessary
+ ochg=chg;
+ idx = ~isinf(maxstep);
+ limit = abs(maxstep(idx).*pprev(idx));
+ chg(idx) = min(max(chg(idx),-limit),limit);
+ if (verbose && any(ochg ~= chg))
+ disp(['Change in parameter(s): ', ...
+ sprintf('%d ',find(ochg ~= chg)), 'maximal fractional stepwidth enforced']);
+ end
end
aprec=abs(pprec.*pbest); %---
%% ss=scalar sum of squares=sum((wt.*f)^2).
if (any(abs(chg) > 0.1*aprec))%--- % only worth evaluating
- % function if there is some non-miniscule
- % change
- p=chg+pprev;
- %% since the projection method may have slightly violated
- %% constraints due to inaccuracy, correct parameters to bounds
- %% --- but only if no further constraints are given, otherwise
- %% the inaccuracy in honoring them might increase by this
- if (~have_constraints_except_bounds)
- lidx = p < lbound;
- uidx = p > ubound;
- p(lidx, 1) = lbound(lidx, 1);
- p(uidx, 1) = ubound(uidx, 1);
- chg(lidx, 1) = p(lidx, 1) - pprev(lidx, 1);
- chg(uidx, 1) = p(uidx, 1) - pprev(uidx, 1);
- end
- %%
- f = F (p);
- r = wt .* f;
- r = r(:);
- if (~isreal (r))
- error ('weighted residuals are not real');
- end
- ss = r.' * r;
- deb_printf (testing, 'sbest: %.16e\n', sbest);
- deb_printf (testing, 'sgoal: %.16e\n', sgoal);
- deb_printf (testing, ' ss: %.16e\n', ss);
- if (ss<sbest)
+ % function if there is some non-miniscule
+ % change
+ p=chg+pprev;
+ %% since the projection method may have slightly violated
+ %% constraints due to inaccuracy, correct parameters to bounds
+ %% --- but only if no further constraints are given, otherwise
+ %% the inaccuracy in honoring them might increase by this
+ if (~have_constraints_except_bounds)
+ lidx = p < lbound;
+ uidx = p > ubound;
+ p(lidx, 1) = lbound(lidx, 1);
+ p(uidx, 1) = ubound(uidx, 1);
+ chg(lidx, 1) = p(lidx, 1) - pprev(lidx, 1);
+ chg(uidx, 1) = p(uidx, 1) - pprev(uidx, 1);
+ end
+ %%
+ f = F (p);
+ r = wt .* f;
+ r = r(:);
+ if (~isreal (r))
+ error ('weighted residuals are not real');
+ end
+ ss = r.' * r;
+ deb_printf (testing, 'sbest: %.16e\n', sbest);
+ deb_printf (testing, 'sgoal: %.16e\n', sgoal);
+ deb_printf (testing, ' ss: %.16e\n', ss);
+ if (ss<sbest)
pbest=p;
fbest=f;
sbest=ss;
- end
- if (ss<=sgoal)
+ end
+ if (ss<=sgoal)
break;
- end
+ end
end %---
end
%% printf ('epsL no.: %i\n', jjj); % for testing
@@ -472,7 +472,7 @@
if (all(abs(chg) <= aprec) && all(abs(chgprev) <= aprec))
cvg = 2;
if (verbose)
- fprintf('Parameter changes converged to specified precision\n');
+ fprintf('Parameter changes converged to specified precision\n');
end
break;
else
Modified: trunk/octave-forge/main/optim/inst/private/__nonlin_residmin__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/private/__nonlin_residmin__.m 2012-05-24 04:38:19 UTC (rev 10509)
+++ trunk/octave-forge/main/optim/inst/private/__nonlin_residmin__.m 2012-05-24 05:50:43 UTC (rev 10510)
@@ -31,6 +31,12 @@
optimget = @ __optimget__;
endif
+ if (compare_versions (version (), "3.2.4", "<="))
+ ## For bug #31484; but Octave 3.6... shows bug #36288 due to this
+ ## workaround. Octave 3.7... seems to be all right.
+ __dfdp__ = @ __dfdp__;
+ endif
+
## some scalar defaults; some defaults are backend specific, so
## lacking elements in respective constructed vectors will be set to
## NA here in the frontend
@@ -56,7 +62,7 @@
"fract_prec", [], \
"diffp", [], \
"diff_onesided", [], \
- "complex_step_derivative", false, \
+ "complex_step_derivative_f", false, \
"complex_step_derivative_inequc", false, \
"complex_step_derivative_equc", false, \
"cstep", cstep_default, \
@@ -112,10 +118,10 @@
diffp = optimget (settings, "diffp");
diff_onesided = optimget (settings, "diff_onesided");
fixed = optimget (settings, "fixed");
- do_cstep = optimget (settings, "complex_step_derivative", false);
+ do_cstep = optimget (settings, "complex_step_derivative_f", false);
cstep = optimget (settings, "cstep", cstep_default);
if (do_cstep && ! isempty (dfdp))
- error ("both 'complex_step_derivative' and 'dfdp' are set");
+ error ("both 'complex_step_derivative_f' and 'dfdp' are set");
endif
do_cstep_inequc = \
optimget (settings, "complex_step_derivative_inequc", false);
@@ -458,7 +464,6 @@
if (do_cstep)
dfdp = @ (p, hook) jacobs (p, f, hook);
else
- __dfdp__ = @ __dfdp__; # for bug #31484 (Octave <= 3.2.4)
dfdp = @ (p, hook) __dfdp__ (p, f, hook);
endif
endif
@@ -814,8 +819,8 @@
tp = eye (sum (nonfixed));
lidx = lbound != - Inf;
uidx = ubound != Inf;
- mc = cat (2, mc, tp(:, lidx), - tp(:, uidx));
- vc = cat (1, vc, - lbound(lidx, 1), ubound(uidx, 1));
+ mc = cat (2, tp(:, lidx), - tp(:, uidx), mc);
+ vc = cat (1, - lbound(lidx, 1), ubound(uidx, 1), vc);
## concatenate linear inequality and equality constraints
mc = cat (2, mc, emc);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
