1 ## Copyright (C) 2012 Olaf Till <i7tiol@t-online.de>
3 ## This program is free software; you can redistribute it and/or modify
4 ## it under the terms of the GNU General Public License as published by
5 ## the Free Software Foundation; either version 3 of the License, or
6 ## (at your option) any later version.
8 ## This program is distributed in the hope that it will be useful,
9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 ## GNU General Public License for more details.
13 ## You should have received a copy of the GNU General Public License
14 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
16 ## The simulated annealing code is translated and adapted from siman.c,
17 ## written by Mark Galassi, of the GNU Scientific Library.
19 function [p_res, objf, cvg, outp] = __sqp__ (f, pin, hook)
21 ## needed for some anonymous functions
22 if (exist ("ifelse") != 5)
23 ifelse = @ scalar_ifelse;
29 mc = hook.mc; # matrix of linear constraints
30 vc = hook.vc; # vector of linear constraints
31 f_cstr = hook.f_cstr; # function of all constraints
32 df_cstr = hook.df_cstr; # function of derivatives of all constraints
33 n_gencstr = hook.n_gencstr; # number of non-linear constraints
34 eq_idx = hook.eq_idx; # logical index of equality constraints in all
36 lbound = hook.lbound; # bounds, subset of linear inequality
37 ubound = hook.ubound; # constraints in mc and vc
39 ## passed values of constraints for initial parameters
40 pin_cstr = hook.pin_cstr;
42 ## passed return value of f for initial parameters
45 ## passed function for gradient of objective function
48 ## passed function for hessian of objective function
49 if (isempty (hessian = hook.hessian))
56 ## passed function for complementary pivoting
62 if (isempty (niter)) niter = 20; endif
65 ## some useful variables derived from passed variables
69 nz = 20 * eps; # This is arbitrary. Accuracy of equality constraints.
71 ## backend-specific checking of options and constraints
73 if (any (pin < lbound | pin > ubound) ||
74 any (pin_cstr.inequ.lin_except_bounds < 0) ||
75 any (pin_cstr.inequ.gen < 0) ||
76 any (abs (pin_cstr.equ.lin)) >= nz ||
77 any (abs (pin_cstr.equ.gen)) >= nz)
78 error ("Initial parameters violate constraints.");
82 error ("no free parameters");
85 ## fill constant fields of hook for derivative-functions; some fields
86 ## may be backend-specific
87 dfdp_hook.fixed = fixed; # this may be handled by the frontend, but
88 # the backend still may add to it
90 ## set up for iterations
105 if (any (isnan (H(:))))
106 error ("some second derivatives undefined by user function");
109 error ("second derivatives given by user function not real");
111 if (! issymmetric (H))
112 error ("Hessian returned by user function not symmetric");
115 R = directional_discrimination (H);
126 function R = directional_discrimination (A)
128 ## A is expected to be real and symmetric without checking
130 ## "Directional discrimination" (Bard, Nonlinear Parameter Estimation,
131 ## Academic Press, 1974). Compute R, which is "similar" to computing
132 ## inv(H), but succeeds even if H is singular; R is positive definite
133 ## and is tuned to avoid large steps of one parameter with respect to
134 ## the steps of the others.
136 ## make matrix Binv for scaling
138 nidx = ! (idx = Binv == 0);
139 Binv(nidx) = 1 ./ sqrt (abs (Binv(nidx)));
143 ## eigendecomposition of scaled hessian
144 [V, L] = eig (Binv * A * Binv);
146 ## A is symmetric, so V and L are real, delete any imaginary parts,
147 ## which might occur due to inaccuracy
151 ## actual directional discrimination, does not exactly follow Bard
152 L = abs (diag (L)); # R should get positive definite
153 L = max (L, .001 * max (L)); # avoids relatively large steps of
158 R = G * diag (1 ./ L) * G.';