Logo ROOT  
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 
200 int GSLSimAnnealing::Solve(GSLSimAnFunc & fx, bool debug) {
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_t;
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 
ROOT::Math::GSLSimAnFunc::Energy
virtual double Energy() const
evaluate the energy ( objective function value) re-implement by derived classes if needed to be modif...
Definition: GSLSimAnnealing.cxx:50
n
const Int_t n
Definition: legend1.C:16
ROOT::Math::GSLSimAnnealing::GSLSimAnnealing
GSLSimAnnealing()
Default constructor.
Definition: GSLSimAnnealing.cxx:174
ROOT::Math::GSLSimAnFunc::GSLSimAnFunc
GSLSimAnFunc()
derived classes might need to re-define completely the class
Definition: GSLSimAnnealing.h:70
ROOT::Math::GSLSimAnParams::k
double k
parameters for the Boltzman distribution
Definition: GSLSimAnnealing.h:180
IFunction.h
f
#define f(i)
Definition: RSha256.hxx:104
dest
#define dest(otri, vertexptr)
Definition: triangle.c:1040
ROOT::Math::GSLSimAn::CopyCtor
void * CopyCtor(void *xp)
Definition: GSLSimAnnealing.cxx:157
ROOT::Math::GSLSimAnFunc::Distance
virtual double Distance(const GSLSimAnFunc &func) const
calculate the distance (metric) between this one and another configuration Presently a cartesian metr...
Definition: GSLSimAnnealing.cxx:66
ROOT::Math::GSLSimAnFunc::fScale
std::vector< double > fScale
Definition: GSLSimAnnealing.h:150
ROOT::Math::GSLSimAn::Copy
void Copy(void *source, void *dest)
Definition: GSLSimAnnealing.cxx:149
ROOT::Math::GSLSimAnFunc::fX
std::vector< double > fX
Definition: GSLSimAnnealing.h:149
r
ROOT::R::TRInterface & r
Definition: Object.C:4
ROOT::Math::GSLRngWrapper
GSLRngWrapper class to wrap gsl_rng structure.
Definition: GSLRngWrapper.h:25
ROOT::Math::GSLSimAnFunc::FastCopy
virtual GSLSimAnFunc & FastCopy(const GSLSimAnFunc &f)
fast copy method called by GSL simuated annealing internally copy only the things which have been cha...
Definition: GSLSimAnnealing.cxx:96
ROOT::Math::GSLSimAnnealing::Solve
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 ...
Definition: GSLSimAnnealing.cxx:183
ROOT::Math::GSLSimAn::Step
void Step(const gsl_rng *r, void *xp, double step_size)
Definition: GSLSimAnnealing.cxx:117
x
Double_t x[n]
Definition: legend1.C:17
ROOT::Math::GSLSimAnFunc::fFunc
const ROOT::Math::IMultiGenFunction * fFunc
Definition: GSLSimAnnealing.h:151
ROOT::Math::GSLSimAnParams::n_tries
int n_tries
Definition: GSLSimAnnealing.h:176
ROOT::Math::GSLSimAnFunc::NDim
unsigned int NDim() const
Definition: GSLSimAnnealing.h:135
ROOT::Math::GSLSimAn::E
double E(void *xp)
Definition: GSLSimAnnealing.cxx:110
ROOT::Math::GSLSimAnParams::step_size
double step_size
Definition: GSLSimAnnealing.h:178
GSLSimAnnealing.h
xmin
float xmin
Definition: THbookFile.cxx:95
ROOT::Math::GSLSimAnnealing::fParams
GSLSimAnParams fParams
Definition: GSLSimAnnealing.h:249
ROOT::Math::GSLSimAnParams::t_initial
double t_initial
Definition: GSLSimAnnealing.h:181
ROOT::Math::GSLSimAn::Print
void Print(void *xp)
Definition: GSLSimAnnealing.cxx:139
ROOT::Math::GSLSimAnFunc::Step
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...
Definition: GSLSimAnnealing.cxx:55
ROOT::Math::GSLSimAnFunc::Print
virtual void Print()
print the position in the standard output std::ostream GSL prints in addition n iteration,...
Definition: GSLSimAnnealing.cxx:83
ROOT::Math::GSLSimAnParams::t_min
double t_min
Definition: GSLSimAnnealing.h:183
double
double
Definition: Converters.cxx:921
y
Double_t y[n]
Definition: legend1.C:17
sqrt
double sqrt(double)
ROOT::Math::GSLSimAn::Dist
double Dist(void *xp, void *yp)
Definition: GSLSimAnnealing.cxx:129
ROOT::Math::GSLSimAn::Destroy
void Destroy(void *xp)
Definition: GSLSimAnnealing.cxx:163
ROOT::Math::GSLSimAnParams::iters_fixed_T
int iters_fixed_T
Definition: GSLSimAnnealing.h:177
ROOT::Math::GSLSimAnFunc
GSLSimAnFunc class description.
Definition: GSLSimAnnealing.h:50
GSLRndmEngines.h
ROOT::Math::GSLSimAnFunc::X
double X(unsigned int i) const
Definition: GSLSimAnnealing.h:137
ROOT::Math::IBaseFunctionMultiDimTempl
Documentation for the abstract class IBaseFunctionMultiDim.
Definition: IFunction.h:62
ROOT::Math::GSLRandomEngine
GSLRandomEngine Base class for all GSL random engines, normally user instantiate the derived classes ...
Definition: GSLRndmEngines.h:65
ROOT
VSD Structures.
Definition: StringConv.hxx:21
GSLRngWrapper.h
Math
Namespace for new Math classes and functions.
ROOT::Math::GSLSimAnFunc::Clone
virtual GSLSimAnFunc * Clone() const
clone method.
Definition: GSLSimAnnealing.h:92
ROOT::Math::GSLSimAnParams::mu_t
double mu_t
Definition: GSLSimAnnealing.h:182