ROOT  6.06/09
Reference Guide
GSLSimAnnealing.cxx
Go to the documentation of this file.
1 // @(#)root/mathmore:$Id$
2 // Author: L. Moneta Thu Jan 25 11:13:48 2007
3 
4 /**********************************************************************
5  * *
6  * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT *
7  * *
8  * *
9  **********************************************************************/
10 
11 // Implementation file for class GSLSimAnnealing
12 
13 #include "Math/GSLSimAnnealing.h"
14 
15 #include "gsl/gsl_siman.h"
16 
17 #include "Math/IFunction.h"
18 #include "Math/GSLRndmEngines.h"
19 #include "GSLRngWrapper.h"
20 
21 
22 #include <cassert>
23 #include <iostream>
24 #include <cmath>
25 #include <vector>
26 
27 namespace ROOT {
28 
29  namespace Math {
30 
31 
32 // implementation of GSLSimAnFunc
33 
35  fX( std::vector<double>(x, x + func.NDim() ) ),
36  fScale( std::vector<double>(func.NDim() )),
37  fFunc(&func)
38 {
39  // set scale factors to 1
40  fScale.assign(fScale.size(), 1.);
41 }
42 
43 GSLSimAnFunc::GSLSimAnFunc(const ROOT::Math::IMultiGenFunction & func, const double * x, const double * scale) :
44  fX( std::vector<double>(x, x + func.NDim() ) ),
45  fScale( std::vector<double>(scale, scale + func.NDim() ) ),
46  fFunc(&func)
47 {}
48 
49 
50 double GSLSimAnFunc::Energy() const {
51  // evaluate the energy
52  return (*fFunc)(&fX.front() );
53 }
54 
55 void GSLSimAnFunc::Step(const GSLRandomEngine & random, double maxstep) {
56  // x -> x + Random[-step,step] for each coordinate
57  unsigned int ndim = NDim();
58  for (unsigned int i = 0; i < ndim; ++i) {
59  double urndm = random();
60  double sstep = maxstep * fScale[i];
61  fX[i] += 2 * sstep * urndm - sstep;
62  }
63 }
64 
65 
66 double GSLSimAnFunc::Distance(const GSLSimAnFunc & f) const {
67  // calculate the distance with respect onother configuration
68  const std::vector<double> & x = fX;
69  const std::vector<double> & y = f.X();
70  unsigned int n = x.size();
71  assert (n == y.size());
72  if (n > 1) {
73  double d2 = 0;
74  for (unsigned int i = 0; i < n; ++i)
75  d2 += ( x[i] - y[i] ) * ( x[i] - y[i] );
76  return std::sqrt(d2);
77  }
78  else
79  // avoid doing a sqrt for 1 dim
80  return std::abs( x[0] - y[0] );
81 }
82 
84  // print the position x in standard std::ostream
85  // GSL prints also niter- ntrials - temperature and then the energy and energy min value (from 1.10)
86  std::cout << "\tx = ( ";
87  unsigned n = NDim();
88  for (unsigned int i = 0; i < n-1; ++i) {
89  std::cout << fX[i] << " , ";
90  }
91  std::cout << fX.back() << " )\t";
92  // energy us printed by GSL (and also end-line)
93  std::cout << "E / E_best = "; // GSL print then E and E best
94 }
95 
97  // copy only the information which is changed during the process
98  // in this case only the x values
99  std::copy(rhs.fX.begin(), rhs.fX.end(), fX.begin() );
100  return *this;
101 }
102 
103 
104 
105 // definition and implementations of the static functions required by GSL
106 
107 namespace GSLSimAn {
108 
109 
110  double E( void * xp) {
111  // evaluate the energy given a state xp
112  GSLSimAnFunc * fx = reinterpret_cast<GSLSimAnFunc *> (xp);
113  assert (fx != 0);
114  return fx->Energy();
115  }
116 
117  void Step( const gsl_rng * r, void * xp, double step_size) {
118  // change xp according to the random step
119  GSLSimAnFunc * fx = reinterpret_cast<GSLSimAnFunc *> (xp);
120  assert (fx != 0);
121  // create GSLRandomEngine class
122  // cast away constness (we make sure we don't delete (no call to Terminate() )
123  GSLRngWrapper rng(const_cast<gsl_rng *>(r));
124  GSLRandomEngine random(&rng);
125  // wrapper classes random and rng must exist during call to Step()
126  fx->Step(random, step_size);
127  }
128 
129  double Dist( void * xp, void * yp) {
130  // calculate the distance between two configuration
131  GSLSimAnFunc * fx = reinterpret_cast<GSLSimAnFunc *> (xp);
132  GSLSimAnFunc * fy = reinterpret_cast<GSLSimAnFunc *> (yp);
133 
134  assert (fx != 0);
135  assert (fy != 0);
136  return fx->Distance(*fy);
137  }
138 
139  void Print(void * xp) {
140  // print the position xp
141  // GSL prints also first niter- ntrials - temperature and then the energy
142  GSLSimAnFunc * fx = reinterpret_cast<GSLSimAnFunc *> (xp);
143  assert (fx != 0);
144  fx->Print();
145  }
146 
147 // static function to pass to GSL copy - create and destroy the object
148 
149  void Copy( void * source, void * dest) {
150  GSLSimAnFunc * fx = reinterpret_cast<GSLSimAnFunc *> (source);
151  assert (fx != 0);
152  GSLSimAnFunc * gx = reinterpret_cast<GSLSimAnFunc *> (dest);
153  assert (gx != 0);
154  gx->FastCopy(*fx);
155  }
156 
157  void * CopyCtor( void * xp) {
158  GSLSimAnFunc * fx = reinterpret_cast<GSLSimAnFunc *> (xp);
159  assert (fx != 0);
160  return static_cast<void *> ( fx->Clone() );
161  }
162 
163  void Destroy( void * xp) {
164  GSLSimAnFunc * fx = reinterpret_cast<GSLSimAnFunc *> (xp);
165  assert (fx != 0);
166  delete fx;
167  }
168 
169 }
170 
171 // implementation of GSLSimAnnealing class
172 
173 
175 {
176  // Default constructor implementation.
177 }
178 
179 
180 
181 // function for solving (from a Genfunction interface)
182 
183 int GSLSimAnnealing::Solve(const ROOT::Math::IMultiGenFunction & func, const double * x0, const double * scale, double * xmin, bool debug) {
184  // solve the simulated annealing problem given starting point and objective function interface
185 
186 
187  // initial conditions
188  GSLSimAnFunc fx(func, x0, scale);
189 
190  int iret = Solve(fx, debug);
191 
192  if (iret == 0) {
193  // copy value of the minimum in xmin
194  std::copy(fx.X().begin(), fx.X().end(), xmin);
195  }
196  return iret;
197 
198 }
199 
201  // solve the simulated annealing problem given starting point and GSLSimAnfunc object
202 
203  gsl_rng * r = gsl_rng_alloc(gsl_rng_mt19937);
204 
205 
206 
207  gsl_siman_params_t simanParams;
208 
209  // parameters for the simulated annealing
210  // copy them in GSL structure
211 
212  simanParams.n_tries = fParams.n_tries; /* how many points to try for each step */
213  simanParams.iters_fixed_T = fParams.iters_fixed_T; /* how many iterations at each temperature? */
214  simanParams.step_size = fParams.step_size; /* max step size in the random walk */
215  // the following parameters are for the Boltzmann distribution */
216  simanParams.k = fParams.k;
217  simanParams.t_initial = fParams.t_initial;
218  simanParams.mu_t = fParams.mu;
219  simanParams.t_min = fParams.t_min;
220 
221 
222  if (debug)
223  gsl_siman_solve(r, &fx, &GSLSimAn::E, &GSLSimAn::Step, &GSLSimAn::Dist,
225 
226  else
227  gsl_siman_solve(r, &fx, &GSLSimAn::E, &GSLSimAn::Step, &GSLSimAn::Dist,
228  0, &GSLSimAn::Copy, &GSLSimAn::CopyCtor , &GSLSimAn::Destroy, 0, simanParams );
229 
230  return 0;
231 
232 }
233 
234 
235 
236 
237  } // end namespace Math
238 
239 } // end namespace ROOT
240 
GSLSimAnnealing()
Default constructor.
float xmin
Definition: THbookFile.cxx:93
Namespace for new ROOT classes and functions.
Definition: ROOT.py:1
virtual void Step(const GSLRandomEngine &r, double maxstep)
change the x[i] value using a random value urndm generated between [0,1] up to a maximum value maxste...
const ROOT::Math::IMultiGenFunction * fFunc
#define assert(cond)
Definition: unittest.h:542
GSLSimAnFunc class description.
virtual double Distance(const GSLSimAnFunc &func) const
calculate the distance (metric) between this one and another configuration Presently a cartesian metr...
void Step(const gsl_rng *r, void *xp, double step_size)
unsigned int NDim() const
STL namespace.
double sqrt(double)
int Solve(const ROOT::Math::IMultiGenFunction &func, const double *x0, const double *scale, double *xmin, bool debug=false)
solve the simulated annealing given a multi-dim function, the initial vector parameters and a vector ...
Double_t x[n]
Definition: legend1.C:17
double k
parameters for the Boltzman distribution
GSLRandomEngine Base class for all GSL random engines, normally user instantiate the derived classes ...
GSLRngWrapper class to wrap gsl_rng structure.
Definition: GSLRngWrapper.h:25
static Vc_ALWAYS_INLINE Vector< T > abs(const Vector< T > &x)
Definition: vector.h:450
virtual void Print()
print the position in the standard output std::ostream GSL prints in addition n iteration, n function calls, temperature and energy re-implement by derived classes if necessary
double Dist(void *xp, void *yp)
virtual double Energy() const
evaluate the energy ( objective function value) re-implement by derived classes if needed to be modif...
std::vector< double > fScale
ROOT::R::TRInterface & r
Definition: Object.C:4
virtual GSLSimAnFunc & FastCopy(const GSLSimAnFunc &f)
fast copy method called by GSL simuated annealing internally copy only the things which have been cha...
void Copy(void *source, void *dest)
double f(double x)
GSLSimAnFunc()
derived classes might need to re-define completely the class
Double_t y[n]
Definition: legend1.C:17
double func(double *x, double *p)
Definition: stressTF1.cxx:213
Namespace for new Math classes and functions.
virtual GSLSimAnFunc * Clone() const
clone method.
#define dest(otri, vertexptr)
Definition: triangle.c:1040
void * CopyCtor(void *xp)
bool debug
std::vector< double > fX
Documentation for the abstract class IBaseFunctionMultiDim.
Definition: IFunction.h:63
const Int_t n
Definition: legend1.C:16
double X(unsigned int i) const