octave_packages/optim-1.2.0/private/__siman__.m

   1 ## Copyright (C) 2012 Olaf Till <i7tiol@t-online.de>
   2 ##
   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.
   7 ##
   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.
  12 ##
  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/>.
  15
  16 ## The simulated annealing code is translated and adapted from siman.c,
  17 ## written by Mark Galassi, of the GNU Scientific Library.
  18
  19 function [p_res, objf, cvg, outp] = __siman__ (f, pin, hook)
  20
  21   ## needed for some anonymous functions
  22   if (exist ("ifelse") != 5)
  23     ifelse = @ scalar_ifelse;
  24   endif
  25
  26   ## passed constraints
  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
  33                                 # constraints
  34   lbound = hook.lbound; # bounds, subset of linear inequality
  35   ubound = hook.ubound; # constraints in mc and vc
  36
  37   ## passed values of constraints for initial parameters
  38   pin_cstr = hook.pin_cstr;
  39
  40   ## passed return value of f for initial parameters
  41   f_pin = hook.f_pin;
  42
  43   ## passed function for complementary pivoting, currently sqp is used
  44   ## instead
  45   ##
  46   ## cpiv = hook.cpiv;
  47
  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;
  54
  55   ## passed options
  56   fixed = hook.fixed;
  57   verbose = strcmp (hook.Display, "iter");
  58   regain_constraints = hook.stoch_regain_constr;
  59   if ((siman_log = hook.siman_log))
  60     log = zeros (0, 5);
  61   endif
  62   if ((trace_steps = hook.trace_steps))
  63     trace = [0, 0, f_pin, pin.'];
  64   endif
  65
  66   ## some useful variables derived from passed variables
  67   n = length (pin);
  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);
  74   ineq_idx = ! eq_idx;
  75   leq_idx = eq_idx(1:n_lconstr);
  76   lineq_idx = ineq_idx(1:n_lconstr);
  77   lfalse_idx = false(n_lconstr, 1);
  78
  79   nz = 20 * eps; # This is arbitrary. Accuracy of equality constraints.
  80
  81   ## backend-specific checking of options and constraints
  82   ##
  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.");
  86   endif
  87   ##
  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.");
  94   endif
  95   ##
  96   if (all (fixed))
  97     error ("no free parameters");
  98   endif
  99   ##
 100   idx = isna (max_rand_step);
 101   max_rand_step(idx) = 0.005 * pin(idx);
 102
 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
 107
 108   ## set up for iterations
 109   sizep = size (pin);
 110   p = best_p = pin;
 111   E = best_E = f_pin;
 112   T = T_init;
 113   n_evals = 1; # one has been done by frontend
 114   n_iter = 0;
 115   done = false;
 116
 117   cvg = 1;
 118
 119   ## simulated annealing
 120   while (! done)
 121
 122     n_iter++;
 123
 124     n_accepts = n_rejects = n_eless = 0;
 125
 126     for id = 1 : iters_fixed_T
 127
 128       new_p = p + max_rand_step .* (2 * rand (sizep) - 1);
 129
 130       ## apply constraints
 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)) == \
 138               nv)
 139             ## special case only bounds violated, set back to bound
 140             new_p(lbvidx) = lbound(lbvidx);
 141             new_p(ubvidx) = ubound(ubvidx);
 142           elseif (nv == 1 && \
 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
 147             tidx = lfalse_idx;
 148             tidx(leq_idx) = t_eq;
 149             tidx(lineq_idx) = t_inequ;
 150             c = mc(:, tidx);
 151             d = ac(tidx);
 152             new_p -= c * (d / (c.' * c));
 153           else
 154             ## other cases, set back keeping the distance to original
 155             ## 'new_p' minimal, using quadratic programming, or
 156             ## sequential quadratic programming for nonlinear
 157             ## constraints
 158             [new_p, discarded, sqp_info] = \
 159                 sqp (new_p, \
 160                      {@(x)sumsq(x-new_p), \
 161                       @(x)2*(x-new_p), \
 162                       @(x)sqp_hessian}, \
 163                      {@(x)f_cstr(x,eq_idx), \
 164                       @(x)df_cstr(x,eq_idx, \
 165                                   setfield(hook,"f", \
 166                                            f_cstr(x,ac_idx)))}, \
 167                      {@(x)f_cstr(x,ineq_idx), \
 168                       @(x)df_cstr(x,ineq_idx, \
 169                                   setfield(hook,"f", \
 170                                            f_cstr(x,ac_idx)))});
 171             if (sqp_info != 101)
 172               cvg = 0;
 173               done = true;
 174               break;
 175             endif
 176           endif
 177         endif
 178       else
 179         n_retry_constr = 0;
 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);
 183           n_retry_constr++;
 184         endwhile
 185         if (verbose && n_retry_constr)
 186           printf ("%i additional tries of random step to meet constraints\n",
 187                   n_retry_constr);
 188         endif
 189       endif
 190
 191       new_E = f (new_p);
 192       n_evals++;
 193
 194       if (new_E < best_E)
 195         best_p = new_p;
 196         best_E = new_E;
 197       endif
 198       if (new_E < E)
 199         ## take a step
 200         p = new_p;
 201         E = new_E;
 202         n_eless++;
 203         if (trace_steps)
 204           trace(end + 1, :) = [n_iter, id, E, p.'];
 205         endif
 206       elseif (rand (1) < exp (- (new_E - E) / T))
 207         ## take a step
 208         p = new_p;
 209         E = new_E;
 210         n_accepts++;
 211         if (trace_steps)
 212           trace(end + 1, :) = [n_iter, id, E, p.'];
 213         endif
 214       else
 215         n_rejects++;
 216       endif
 217
 218     endfor # iters_fixed_T
 219
 220     if (verbose)
 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);
 224     endif
 225
 226     if (siman_log)
 227       log(end + 1, :) = [T, E, n_eless, n_accepts, n_rejects];
 228     endif
 229
 230     ## cooling
 231     T /= mu_T;
 232     if (T < T_min)
 233       done = true;
 234     endif
 235
 236   endwhile
 237
 238   ## return result
 239   p_res = best_p;
 240   objf = best_E;
 241   outp.niter = n_iter;
 242   if (trace_steps)
 243     outp.trace = trace;
 244   endif
 245   if (siman_log)
 246     outp.log = log;
 247   endif
 248
 249 endfunction