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] = __siman__ (f, pin, hook)
21 ## needed for some anonymous functions
22 if (exist ("ifelse") != 5)
23 ifelse = @ scalar_ifelse;
27 mc = hook.mc; # matrix of linear constraints
28 vc = hook.vc; # vector of linear constraints
29 f_cstr = hook.f_cstr; # function of all constraints
30 df_cstr = hook.df_cstr; # function of derivatives of all constraints
31 n_gencstr = hook.n_gencstr; # number of non-linear constraints
32 eq_idx = hook.eq_idx; # logical index of equality constraints in all
34 lbound = hook.lbound; # bounds, subset of linear inequality
35 ubound = hook.ubound; # constraints in mc and vc
37 ## passed values of constraints for initial parameters
38 pin_cstr = hook.pin_cstr;
40 ## passed return value of f for initial parameters
43 ## passed function for complementary pivoting, currently sqp is used
48 ## passed simulated annealing parameters
49 T_init = hook.siman.T_init;
50 T_min = hook.siman.T_min;
51 mu_T = hook.siman.mu_T;
52 iters_fixed_T = hook.siman.iters_fixed_T;
53 max_rand_step = hook.max_rand_step;
57 verbose = strcmp (hook.Display, "iter");
58 regain_constraints = hook.stoch_regain_constr;
59 if ((siman_log = hook.siman_log))
62 if ((trace_steps = hook.trace_steps))
63 trace = [0, 0, f_pin, pin.'];
66 ## some useful variables derived from passed variables
68 sqp_hessian = 2 * eye (n);
69 n_lconstr = length (vc);
70 n_bounds = sum (lbound != -Inf) + sum (ubound != Inf);
71 bidx = false (n_lconstr + n_gencstr, 1);
72 bidx(1 : n_bounds) = true;
73 ac_idx = true (n_lconstr + n_gencstr, 1);
75 leq_idx = eq_idx(1:n_lconstr);
76 lineq_idx = ineq_idx(1:n_lconstr);
77 lfalse_idx = false(n_lconstr, 1);
79 nz = 20 * eps; # This is arbitrary. Accuracy of equality constraints.
81 ## backend-specific checking of options and constraints
83 ## equality constraints can not be met by chance
84 if ((any (eq_idx) || any (lbound == ubound)) && ! regain_constraints)
85 error ("If 'stoch_regain_constr' is not set, equality constraints or identical lower and upper bounds are not allowed by simulated annealing backend.");
88 if (any (pin < lbound | pin > ubound) ||
89 any (pin_cstr.inequ.lin_except_bounds < 0) ||
90 any (pin_cstr.inequ.gen < 0) ||
91 any (abs (pin_cstr.equ.lin)) >= nz ||
92 any (abs (pin_cstr.equ.gen)) >= nz)
93 error ("Initial parameters violate constraints.");
97 error ("no free parameters");
100 idx = isna (max_rand_step);
101 max_rand_step(idx) = 0.005 * pin(idx);
103 ## fill constant fields of hook for derivative-functions; some fields
104 ## may be backend-specific
105 dfdp_hook.fixed = fixed; # this may be handled by the frontend, but
106 # the backend still may add to it
108 ## set up for iterations
113 n_evals = 1; # one has been done by frontend
119 ## simulated annealing
124 n_accepts = n_rejects = n_eless = 0;
126 for id = 1 : iters_fixed_T
128 new_p = p + max_rand_step .* (2 * rand (sizep) - 1);
131 if (regain_constraints)
132 evidx = (abs ((ac = f_cstr (new_p, ac_idx))(eq_idx)) >= nz);
133 ividx = (ac(ineq_idx) < 0);
134 if (any (evidx) || any (ividx))
135 nv = sum (evidx) + sum (ividx);
136 if (sum (lbvidx = (new_p < lbound)) + \
137 sum (ubvidx = (new_p > ubound)) == \
139 ## special case only bounds violated, set back to bound
140 new_p(lbvidx) = lbound(lbvidx);
141 new_p(ubvidx) = ubound(ubvidx);
143 sum (t_eq = (abs (ac(leq_idx)) >= nz)) + \
144 sum (t_inequ = (ac(lineq_idx) < 0)) == 1)
145 ## special case only one linear constraint violated, set
146 ## back perpendicularly to constraint
148 tidx(leq_idx) = t_eq;
149 tidx(lineq_idx) = t_inequ;
152 new_p -= c * (d / (c.' * c));
154 ## other cases, set back keeping the distance to original
155 ## 'new_p' minimal, using quadratic programming, or
156 ## sequential quadratic programming for nonlinear
158 [new_p, discarded, sqp_info] = \
160 {@(x)sumsq(x-new_p), \
163 {@(x)f_cstr(x,eq_idx), \
164 @(x)df_cstr(x,eq_idx, \
166 f_cstr(x,ac_idx)))}, \
167 {@(x)f_cstr(x,ineq_idx), \
168 @(x)df_cstr(x,ineq_idx, \
170 f_cstr(x,ac_idx)))});
180 while (any (abs ((ac = f_cstr (new_p, ac_idx))(eq_idx)) >= nz) \
181 || any (ac(ineq_idx) < 0))
182 new_p = p + max_rand_step .* (2 * rand (sizep) - 1);
185 if (verbose && n_retry_constr)
186 printf ("%i additional tries of random step to meet constraints\n",
204 trace(end + 1, :) = [n_iter, id, E, p.'];
206 elseif (rand (1) < exp (- (new_E - E) / T))
212 trace(end + 1, :) = [n_iter, id, E, p.'];
218 endfor # iters_fixed_T
221 printf ("temperature no. %i: %e, energy %e,\n", n_iter, T, E);
222 printf ("tries with energy less / not less but accepted / rejected:\n");
223 printf ("%i / %i / %i\n", n_eless, n_accepts, n_rejects);
227 log(end + 1, :) = [T, E, n_eless, n_accepts, n_rejects];