anneal_mt.h
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 #ifndef O2SCL_ANNEAL_MT_H
24 #define O2SCL_ANNEAL_MT_H
25 
26 #include <o2scl/rng_gsl.h>
27 #include <o2scl/multi_funct.h>
28 #include <o2scl/anneal.h>
29 #include <o2scl/vector.h>
30 
31 // for boost::bind()
32 #include <boost/bind.hpp>
33 // for boost::thread
34 #include <boost/thread/thread.hpp>
35 
36 #ifndef DOXYGEN_NO_O2NS
37 namespace o2scl {
38 #endif
39 
40  /** \brief Multidimensional minimization by simulated annealing
41  (Boost multi-threaded version)
42 
43  This header-only class additionally requires the Boost
44  libraries. It performs simulated annealing using an arbitrary
45  number of processors using <tt>boost::thread</tt>, which is
46  closely related to the standard Unix pthread library. It works
47  very similarly to \ref anneal_gsl, it performs \ref ntrial
48  evaluations over each processor, then applying the metropolis
49  algorithm to the results from all of the processors at the end.
50 
51  Because <tt>np</tt> function calls are happening simultaneously,
52  where <tt>np</tt> is the number of processors, <tt>np</tt>
53  copies of the function parameters of type <tt>param_t</tt> must
54  also be specified. The user-defined function to minimize must
55  also be thread-safe, allowing multiple threads to call the
56  function at once (albeit given different parameters). The
57  default type for these <tt>np</tt> copies of the parameters of
58  type <tt>param_t</tt> is <tt>std::vector<param_t></tt>.
59 
60  This works particularly well for functions which are not trivial
61  to evaluate, i.e. functions where the execution time is more
62  longer than the bookkeeping that \ref anneal_mt performs between
63  trials. For functions which satisfy this requirement, this
64  algorithm scales nearly linearly with the number of processors.
65 
66  Verbose I/O for this class happens only outside the theads
67  unless the user places I/O in the streams in the function that
68  is specified.
69 
70  \future There may be a good way to remove the function indirection
71  here to make this class a bit faster.
72  */
73  template<class func_t=multi_funct11,
75  class rng_t=int, class rng_dist_t=rng_gsl>
76  class anneal_mt : public anneal_base<func_t,vec_t,rng_t,rng_dist_t> {
77 
78  public:
79 
80  anneal_mt() {
81  ntrial=100;
82  nproc=1;
83  outs=&std::cout;
84  ins=&std::cin;
85  tolx=1.0e-6;
86  T_start=1.0;
87  T_dec=1.5;
88  step_dec=1.5;
89  min_step_ratio=100.0;
90  step_vec.resize(1);
91  step_vec[0]=1.0;
92  out_best=false;
93  out_step_changes=false;
94  }
95 
96  virtual ~anneal_mt() {
97  }
98 
99  /// \name Basic usage
100  //@{
101  /** \brief Calculate the minimum \c fmin of \c func w.r.t the
102  array \c x0 of size \c nv using \c np threads.
103  */
104  virtual int mmin(size_t nv, vec_t &x0, double &fmin,
105  func_t &func, size_t np) {
106 
107  if (nv==0) {
108  O2SCL_ERR2("Tried to minimize over zero variables ",
109  " in anneal_mt::mmin().",o2scl::exc_einval);
110  }
111  if (np==0) {
112  O2SCL_ERR2("Tried to use zero threads in ",
113  "anneal_mt::mmin().",o2scl::exc_einval);
114  }
115 
116  allocate(nv,np);
117  f=&func;
118  fmin=0.0;
119 
120  double E, best_E, T, old_E;
121  int i, iter=0;
122  size_t j;
123 
124  for(j=0;j<nv;j++) {
125  x[j]=x0[j];
126  best_x[j]=x0[j];
127  }
128 
129  E=func(nv,x);
130  best_E=E;
131 
132  // Setup initial temperature and step sizes
133  start(nv,T);
134 
135  bool done=false;
136 
137  boost::thread **thrd;
138 
139  while (!done) {
140 
141  // Copy old value of x for next() function
142 
143  for(j=0;j<nv;j++) old_x[j]=x[j];
144  old_E=E;
145 
146  size_t nmoves=0;
147 
148  for (i=0;i<ntrial;++i) {
149 
150  // Determine the stepsize, create and execute the new threads
151  thrd=new boost::thread *[np];
152  for(size_t ip=0;ip<np;ip++) {
153  new_x[ip]=x;
154  thrd[ip]=new boost::thread
155  (boost::bind(&anneal_mt::func_wrapper,this,ip));
156  }
157  // Wait until all the threads are done
158  for(size_t ip=0;ip<np;ip++) {
159  thrd[ip]->join();
160  }
161  // Delete the threads and continue
162  for(size_t ip=0;ip<np;ip++) {
163  delete thrd[ip];
164  }
165  delete[] thrd;
166  // Process the results from each thread
167  for(size_t ip=0;ip<np;ip++) {
168 
169  // Store best value obtained so far
170  if(new_E[ip]<=best_E){
171  for(j=0;j<nv;j++) best_x[j]=new_x[ip][j];
172  best_E=new_E[ip];
173  if (this->verbose>0 && out_best) {
174  std::cout << "Best: ";
175  o2scl::vector_out(std::cout,best_x);
176  std::cout << " " << best_E << std::endl;
177  }
178  }
179 
180  // Take the crucial step: see if the new point is accepted
181  // or not, as determined by the boltzmann probability
182  if (new_E[ip]<E) {
183  for(j=0;j<nv;j++) x[j]=new_x[ip][j];
184  E=new_E[ip];
185  nmoves++;
186  } else if (this->rng_dist(this->rng)
187  < exp(-(new_E[ip]-E)/(boltz*T)) ) {
188  for(j=0;j<nv;j++) x[j]=new_x[ip][j];
189  E=new_E[ip];
190  nmoves++;
191  }
192  }
193 
194  }
195 
196  if (this->verbose>0) {
197  print_iter(nv,best_x,best_E,iter,T,"anneal_mt");
198  iter++;
199  }
200 
201  // See if we're finished and proceed to the next step
202  next(nv,old_x,old_E,x,E,T,nmoves,done);
203 
204  }
205 
206  for(j=0;j<nv;j++) x0[j]=best_x[j];
207  fmin=best_E;
208 
209  return 0;
210  }
211  //@}
212 
213  /** \brief Desc
214  */
215  virtual int mmin(size_t nv, vec_t &x0, double &fmin,
216  func_t &func) {
217  return mmin(nv,x0,fmin,func,1);
218  }
219 
220  /// \name Iteration control
221  //@{
222  /// Determine how to change the minimization for the next iteration
223  virtual int next(size_t nv, vec_t &x_old, double min_old, vec_t &x_new,
224  double min_new, double &T, size_t n_moves,
225  bool &finished) {
226 
227  if (T/T_dec<this->tolx) {
228  finished=true;
229  return o2scl::success;
230  }
231  if (n_moves==0) {
232  bool changed=false;
233  for(size_t i=0;i<nv;i++) {
234  if (i<step_vec.size() && step_vec[i]>this->tolx*min_step_ratio) {
235  step_vec[i]/=step_dec;
236  changed=true;
237  }
238  }
239  if (changed && verbose>0 && out_step_changes) {
240  std::cout << "Step sizes changed: ";
241  o2scl::vector_out(std::cout,step_vec,true);
242  }
243  }
244  T/=T_dec;
245  return o2scl::success;
246  }
247 
248  /// Setup initial temperature and stepsize
249  virtual int start(size_t nv, double &T) {
250  T=T_start;
251  return o2scl::success;
252  }
253  //@}
254 
255  /// \name Parameters
256  //@{
257  /// Boltzmann factor (default 1.0).
258  double boltz;
259 
260  /// Number of iterations
261  int ntrial;
262 
263  /// Output control
264  int verbose;
265 
266  /// Output step size changes (default false)
268 
269  /// Output best point (default false)
270  bool out_best;
271 
272  /// The independent variable tolerance (default \f$ 10^{-6} \f$ )
273  double tolx;
274 
275  /// Initial temperature (default 1.0)
276  double T_start;
277 
278  /// Factor to decrease temperature by (default 1.5)
279  double T_dec;
280 
281  /// Factor to decrease step size by (default 1.5)
282  double step_dec;
283 
284  /// Ratio between minimum step size and \ref tolx (default 100.0)
286  //@}
287 
288  /// The default random number generator
289  rng_t def_rng;
290 
291  /// Return string denoting type ("anneal_mt")
292  virtual const char *type() { return "anneal_mt"; }
293 
294  /// Set streams for verbose I/O
295  int set_verbose_stream(std::ostream &out, std::istream &in) {
296  outs=&out;
297  ins=&in;
298  return 0;
299  }
300 
301  /** \brief Print out iteration information.
302 
303  Depending on the value of the variable verbose, this prints out
304  the iteration information. If verbose=0, then no information is
305  printed. If verbose>0, then after each iteration, the present
306  values of x and y are output to std::cout along with the
307  iteration number. Also, if verbose>0, every time a new smallest
308  function value is found, the location and the function value is
309  output. If verbose>=2 then each iteration waits for a character
310  between each trial.
311  */
312  virtual int print_iter(size_t nv, vec_t &xx, double y, int iter,
313  double tptr, std::string comment)
314  {
315  if (this->verbose<=0) return 0;
316 
317  size_t i;
318  char ch;
319 
320  (*this->outs) << comment << " Iteration: " << iter << std::endl;
321  std::cout << "x: ";
322  for(i=0;i<nv;i++) std::cout << x[i] << " ";
323  std::cout << std::endl;
324  (*this->outs) << "y: " << y << " Tptr: " << tptr << std::endl;
325  if (this->verbose>1) {
326  (*this->outs) << "Press a key and type enter to continue. ";
327  (*this->ins) >> ch;
328  }
329 
330  return 0;
331  }
332 
333  /// Set the step sizes
334  template<class vec2_t> int set_step(size_t nv, vec2_t &stepv) {
335  if (nv>0) {
336  step_vec.resize(nv);
337  for(size_t i=0;i<nv;i++) step_vec[i]=stepv[i];
338  }
339  return 0;
340  }
341 
342 #ifndef DOXYGEN_INTERNAL
343 
344  protected:
345 
346  /// The function wrapper executed by thread with index \c ip
347  void func_wrapper(size_t ip) {
348  step(nvar,new_x[ip]);
349  new_E[ip]=(*f)(nvar,new_x[ip]);
350  }
351 
352  /// Stream for verbose output
353  std::ostream *outs;
354 
355  /// Stream for verbose input
356  std::istream *ins;
357 
358  /// The number of threads to run
359  size_t nproc;
360 
361  /// The number of variables over which we minimize
362  size_t nvar;
363 
364  /// The function to minimize
365  func_t *f;
366 
367  /// \name Storage for present, next, and best vectors
368  //@{
369  vec_t x, best_x, new_E, old_x;
370  std::vector<vec_t> new_x;
371  //@}
372 
373  /// Vector of step sizes
374  vec_t step_vec;
375 
376  /** \brief Allocate memory for a minimizer over \c n dimensions
377  with stepsize \c step
378  */
379  virtual int allocate(size_t nv, size_t np) {
380  nvar=nv;
381  nproc=np;
382  x.resize(nvar);
383  new_x.resize(np);
384  old_x.resize(nvar);
385  for(size_t i=0;i<np;i++) {
386  new_x[i].resize(nvar);
387  }
388  new_E.resize(np);
389  best_x.resize(nvar);
390  return 0;
391  }
392 
393  /// Make a step to a new attempted minimum
394  virtual int step(size_t nv, vec_t &sx) {
395  size_t nstep=step_vec.size();
396  for(size_t i=0;i<nv;i++) {
397  double u=this->rng_dist(this->rng);
398  sx[i]=(2.0*u-1.0)*step_vec[i%nstep]+sx[i];
399  }
400  return 0;
401  }
402 
403 #endif
404 
405  };
406 
407 #ifndef DOXYGEN_NO_O2NS
408 }
409 #endif
410 
411 #endif
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
size_t nproc
The number of threads to run.
Definition: anneal_mt.h:359
rng_t def_rng
The default random number generator.
Definition: anneal_mt.h:289
double T_dec
Factor to decrease temperature by (default 1.5)
Definition: anneal_mt.h:279
double T_start
Initial temperature (default 1.0)
Definition: anneal_mt.h:276
invalid argument supplied by user
Definition: err_hnd.h:59
virtual int start(size_t nv, double &T)
Setup initial temperature and stepsize.
Definition: anneal_mt.h:249
std::istream * ins
Stream for verbose input.
Definition: anneal_mt.h:356
int verbose
Output control.
Definition: anneal_mt.h:264
virtual int print_iter(size_t nv, vec_t &xx, double y, int iter, double tptr, std::string comment)
Print out iteration information.
Definition: anneal_mt.h:312
double tolx
The independent variable tolerance (default )
Definition: anneal_mt.h:273
virtual int mmin(size_t nv, vec_t &x0, double &fmin, func_t &func, size_t np)
Calculate the minimum fmin of func w.r.t the array x0 of size nv using np threads.
Definition: anneal_mt.h:104
double min_step_ratio
Ratio between minimum step size and tolx (default 100.0)
Definition: anneal_mt.h:285
virtual int step(size_t nv, vec_t &sx)
Make a step to a new attempted minimum.
Definition: anneal_mt.h:394
double boltz
Boltzmann factor (default 1.0).
Definition: anneal_mt.h:258
int set_step(size_t nv, vec2_t &stepv)
Set the step sizes.
Definition: anneal_mt.h:334
virtual const char * type()
Return string denoting type ("anneal_mt")
Definition: anneal_mt.h:292
Simulated annealing base.
Definition: anneal.h:67
virtual int allocate(size_t nv, size_t np)
Allocate memory for a minimizer over n dimensions with stepsize step.
Definition: anneal_mt.h:379
void func_wrapper(size_t ip)
The function wrapper executed by thread with index ip.
Definition: anneal_mt.h:347
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
bool out_best
Output best point (default false)
Definition: anneal_mt.h:270
int ntrial
Number of iterations.
Definition: anneal_mt.h:261
double step_dec
Factor to decrease step size by (default 1.5)
Definition: anneal_mt.h:282
int set_verbose_stream(std::ostream &out, std::istream &in)
Set streams for verbose I/O.
Definition: anneal_mt.h:295
size_t nvar
The number of variables over which we minimize.
Definition: anneal_mt.h:362
bool out_step_changes
Output step size changes (default false)
Definition: anneal_mt.h:267
void vector_out(std::ostream &os, size_t n, const vec_t &v, bool endline=false)
Output the first n elements of a vector to a stream, os.
Definition: vector.h:2319
rng_t rng
The default random number generator.
Definition: anneal.h:120
Multidimensional minimization by simulated annealing (Boost multi-threaded version) ...
Definition: anneal_mt.h:76
virtual int next(size_t nv, vec_t &x_old, double min_old, vec_t &x_new, double min_new, double &T, size_t n_moves, bool &finished)
Determine how to change the minimization for the next iteration.
Definition: anneal_mt.h:223
vec_t step_vec
Vector of step sizes.
Definition: anneal_mt.h:374
std::ostream * outs
Stream for verbose output.
Definition: anneal_mt.h:353
std::function< double(size_t, const boost::numeric::ublas::vector< double > &)> multi_funct11
Multi-dimensional function typedef.
Definition: multi_funct.h:45
func_t * f
The function to minimize.
Definition: anneal_mt.h:365
virtual int mmin(size_t nv, vec_t &x0, double &fmin, func_t &func)
Desc.
Definition: anneal_mt.h:215
Success.
Definition: err_hnd.h:47

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).