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
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
27namespace 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
43GSLSimAnFunc::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
50double GSLSimAnFunc::Energy() const {
51 // evaluate the energy
52 return (*fFunc)(&fX.front() );
53}
54
55void 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
66double 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
107namespace 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
183int 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_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
#define f(i)
Definition: RSha256.hxx:104
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t dest
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
float xmin
Definition: THbookFile.cxx:95
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
GSLSimAnFunc class description.
double X(unsigned int i) const
virtual double Distance(const GSLSimAnFunc &func) const
calculate the distance (metric) between this one and another configuration Presently a cartesian metr...
std::vector< double > fX
virtual void Print()
print the position in the standard output std::ostream GSL prints in addition n iteration,...
GSLSimAnFunc()
derived classes might need to re-define completely the class
unsigned int NDim() const
virtual GSLSimAnFunc & FastCopy(const GSLSimAnFunc &f)
fast copy method called by GSL simuated annealing internally copy only the things which have been cha...
virtual double Energy() const
evaluate the energy ( objective function value) re-implement by derived classes if needed to be modif...
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...
virtual GSLSimAnFunc * Clone() const
clone method.
std::vector< double > fScale
const ROOT::Math::IMultiGenFunction * fFunc
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 ...
GSLSimAnnealing()
Default constructor.
Documentation for the abstract class IBaseFunctionMultiDim.
Definition: IFunction.h:62
RVec< PromoteType< T > > abs(const RVec< T > &v)
Definition: RVec.hxx:1739
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
Namespace for new Math classes and functions.
double Dist(void *xp, void *yp)
void * CopyCtor(void *xp)
void Copy(void *source, void *dest)
void Step(const gsl_rng *r, void *xp, double step_size)
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
double k
parameters for the Boltzman distribution