Logo ROOT   6.14/05
Reference Guide
SimulatedAnnealing.cxx
Go to the documentation of this file.
1 // @(#)root/tmva $Id$
2 // Author: Andreas Hoecker, Joerg Stelzer, Helge Voss, Krzysztof Danielowski, Kamil Kraszewski, Maciej Kruk
3 
4 /**********************************************************************************
5  * Project: TMVA - a Root-integrated toolkit for multivariate data analysis *
6  * Package: TMVA *
7  * Class : SimulatedAnnealing *
8  * Web : http://tmva.sourceforge.net *
9  * *
10  * Description: *
11  * Implementation (see header for description) *
12  * *
13  * Authors (alphabetical): *
14  * Krzysztof Danielowski <danielow@cern.ch> - IFJ & AGH, Poland *
15  * Kamil Kraszewski <kalq@cern.ch> - IFJ & UJ, Poland *
16  * Maciej Kruk <mkruk@cern.ch> - IFJ & AGH, Poland *
17  * *
18  * Copyright (c) 2008: *
19  * IFJ-Krakow, Poland *
20  * CERN, Switzerland *
21  * MPI-K Heidelberg, Germany *
22  * *
23  * Redistribution and use in source and binary forms, with or without *
24  * modification, are permitted according to the terms listed in LICENSE *
25  * (http://tmva.sourceforge.net/LICENSE) *
26  **********************************************************************************/
27 
28 /*! \class TMVA::SimulatedAnnealing
29 \ingroup TMVA
30 Base implementation of simulated annealing fitting procedure.
31 */
32 
34 
35 #include "TMVA/IFitterTarget.h"
36 #include "TMVA/Interval.h"
37 #include "TMVA/GeneticRange.h"
38 #include "TMVA/MsgLogger.h"
39 #include "TMVA/Timer.h"
40 #include "TMVA/Types.h"
41 
42 #include "TRandom3.h"
43 #include "TMath.h"
44 
46 
47 ////////////////////////////////////////////////////////////////////////////////
48 /// constructor
49 
50 TMVA::SimulatedAnnealing::SimulatedAnnealing( IFitterTarget& target, const std::vector<Interval*>& ranges )
51 : fKernelTemperature (kIncreasingAdaptive),
52  fFitterTarget ( target ),
53  fRandom ( new TRandom3(100) ),
54  fRanges ( ranges ),
55  fMaxCalls ( 500000 ),
56  fInitialTemperature ( 1000 ),
57  fMinTemperature ( 0 ),
58  fEps ( 1e-10 ),
59  fTemperatureScale ( 0.06 ),
60  fAdaptiveSpeed ( 1.0 ),
61  fTemperatureAdaptiveStep( 0.0 ),
62  fUseDefaultScale ( kFALSE ),
63  fUseDefaultTemperature ( kFALSE ),
64  fLogger( new MsgLogger("SimulatedAnnealing") ),
65  fProgress(0.0)
66 {
68 }
69 
70 ////////////////////////////////////////////////////////////////////////////////
71 /// option setter
72 
74  Double_t initialTemperature,
75  Double_t minTemperature,
76  Double_t eps,
77  TString kernelTemperatureS,
78  Double_t temperatureScale,
79  Double_t adaptiveSpeed,
80  Double_t temperatureAdaptiveStep,
81  Bool_t useDefaultScale,
82  Bool_t useDefaultTemperature)
83 {
84  fMaxCalls = maxCalls;
85  fInitialTemperature = initialTemperature;
86  fMinTemperature = minTemperature;
87  fEps = eps;
88 
89  if (kernelTemperatureS == "IncreasingAdaptive") {
91  Log() << kINFO << "Using increasing adaptive algorithm" << Endl;
92  }
93  else if (kernelTemperatureS == "DecreasingAdaptive") {
95  Log() << kINFO << "Using decreasing adaptive algorithm" << Endl;
96  }
97  else if (kernelTemperatureS == "Sqrt") {
99  Log() << kINFO << "Using \"Sqrt\" algorithm" << Endl;
100  }
101  else if (kernelTemperatureS == "Homo") {
103  Log() << kINFO << "Using \"Homo\" algorithm" << Endl;
104  }
105  else if (kernelTemperatureS == "Log") {
107  Log() << kINFO << "Using \"Log\" algorithm" << Endl;
108  }
109  else if (kernelTemperatureS == "Sin") {
111  Log() << kINFO << "Using \"Sin\" algorithm" << Endl;
112  }
113 
114  fTemperatureScale = temperatureScale;
115  fAdaptiveSpeed = adaptiveSpeed;
116  fTemperatureAdaptiveStep = temperatureAdaptiveStep;
117 
118  fUseDefaultScale = useDefaultScale;
119  fUseDefaultTemperature = useDefaultTemperature;
120 }
121 
122 ////////////////////////////////////////////////////////////////////////////////
123 /// destructor
124 
126 {
127 }
128 
129 ////////////////////////////////////////////////////////////////////////////////
130 /// random starting parameters
131 
132 void TMVA::SimulatedAnnealing::FillWithRandomValues( std::vector<Double_t>& parameters )
133 {
134  for (UInt_t rIter = 0; rIter < parameters.size(); rIter++) {
135  parameters[rIter] = fRandom->Uniform(0.0,1.0)*(fRanges[rIter]->GetMax() - fRanges[rIter]->GetMin()) + fRanges[rIter]->GetMin();
136  }
137 }
138 
139 ////////////////////////////////////////////////////////////////////////////////
140 /// copy parameters
141 
142 void TMVA::SimulatedAnnealing::ReWriteParameters( std::vector<Double_t>& from, std::vector<Double_t>& to)
143 {
144  for (UInt_t rIter = 0; rIter < from.size(); rIter++) to[rIter] = from[rIter];
145 }
146 
147 ////////////////////////////////////////////////////////////////////////////////
148 /// generate adjacent parameters
149 
150 void TMVA::SimulatedAnnealing::GenerateNeighbour( std::vector<Double_t>& parameters, std::vector<Double_t>& oldParameters,
151  Double_t currentTemperature )
152 {
153  ReWriteParameters( parameters, oldParameters );
154 
155  for (UInt_t rIter=0;rIter<parameters.size();rIter++) {
156  Double_t uni,distribution,sign;
157  do {
158  uni = fRandom->Uniform(0.0,1.0);
159  sign = (uni - 0.5 >= 0.0) ? (1.0) : (-1.0);
160  distribution = currentTemperature * (TMath::Power(1.0 + 1.0/currentTemperature, TMath::Abs(2.0*uni - 1.0)) -1.0)*sign;
161  parameters[rIter] = oldParameters[rIter] + (fRanges[rIter]->GetMax()-fRanges[rIter]->GetMin())*0.1*distribution;
162  }
163  while (parameters[rIter] < fRanges[rIter]->GetMin() || parameters[rIter] > fRanges[rIter]->GetMax() );
164  }
165 }
166 ////////////////////////////////////////////////////////////////////////////////
167 /// generate adjacent parameters
168 
169 std::vector<Double_t> TMVA::SimulatedAnnealing::GenerateNeighbour( std::vector<Double_t>& parameters, Double_t currentTemperature )
170 {
171  std::vector<Double_t> newParameters( fRanges.size() );
172 
173  for (UInt_t rIter=0; rIter<parameters.size(); rIter++) {
174  Double_t uni,distribution,sign;
175  do {
176  uni = fRandom->Uniform(0.0,1.0);
177  sign = (uni - 0.5 >= 0.0) ? (1.0) : (-1.0);
178  distribution = currentTemperature * (TMath::Power(1.0 + 1.0/currentTemperature, TMath::Abs(2.0*uni - 1.0)) -1.0)*sign;
179  newParameters[rIter] = parameters[rIter] + (fRanges[rIter]->GetMax()-fRanges[rIter]->GetMin())*0.1*distribution;
180  }
181  while (newParameters[rIter] < fRanges[rIter]->GetMin() || newParameters[rIter] > fRanges[rIter]->GetMax() );
182  }
183 
184  return newParameters;
185 }
186 
187 ////////////////////////////////////////////////////////////////////////////////
188 /// generate new temperature
189 
191 {
192  if (fKernelTemperature == kSqrt) {
193  currentTemperature = fInitialTemperature/(Double_t)TMath::Sqrt(Iter+2) * fTemperatureScale;
194  }
195  else if (fKernelTemperature == kLog) {
196  currentTemperature = fInitialTemperature/(Double_t)TMath::Log(Iter+2) * fTemperatureScale;
197  }
198  else if (fKernelTemperature == kHomo) {
199  currentTemperature = fInitialTemperature/(Double_t)(Iter+2) * fTemperatureScale;
200  }
201  else if (fKernelTemperature == kSin) {
202  currentTemperature = (TMath::Sin( (Double_t)Iter / fTemperatureScale ) + 1.0 )/ (Double_t)(Iter+1.0) * fInitialTemperature + fEps;
203  }
204  else if (fKernelTemperature == kGeo) {
205  currentTemperature = currentTemperature*fTemperatureScale;
206  }
209  }
211  currentTemperature = currentTemperature*fTemperatureScale;
212  }
213  else Log() << kFATAL << "No such kernel!" << Endl;
214 }
215 
216 ////////////////////////////////////////////////////////////////////////////////
217 /// result checker
218 
219 Bool_t TMVA::SimulatedAnnealing::ShouldGoIn( Double_t currentFit, Double_t localFit, Double_t currentTemperature )
220 {
221  if (currentTemperature < fEps) return kFALSE;
222  Double_t lim = TMath::Exp( -TMath::Abs( currentFit - localFit ) / currentTemperature );
223  Double_t prob = fRandom->Uniform(0.0, 1.0);
224  return (prob < lim) ? kTRUE : kFALSE;
225 }
226 
227 ////////////////////////////////////////////////////////////////////////////////
228 /// setting of default scale
229 
231 {
233  else if (fKernelTemperature == kLog) fTemperatureScale = 1.0;
234  else if (fKernelTemperature == kHomo) fTemperatureScale = 1.0;
235  else if (fKernelTemperature == kSin) fTemperatureScale = 20.0;
236  else if (fKernelTemperature == kGeo) fTemperatureScale = 0.99997;
238  fTemperatureScale = 1.0;
241  fTemperatureScale -= 0.000001;
242  }
243  }
244  else if (fKernelTemperature == kIncreasingAdaptive) fTemperatureScale = 0.15*( 1.0 / (Double_t)(fRanges.size() ) );
245  else Log() << kFATAL << "No such kernel!" << Endl;
246 }
247 
248 ////////////////////////////////////////////////////////////////////////////////
249 /// maximum temperature
250 
251 Double_t TMVA::SimulatedAnnealing::GenerateMaxTemperature( std::vector<Double_t>& parameters )
252 {
253  Int_t equilibrium;
254  Bool_t stopper = 0;
255  Double_t t, dT, cold, delta, deltaY, y, yNew, yBest, yOld;
256  std::vector<Double_t> x( fRanges.size() ), xNew( fRanges.size() ), xBest( fRanges.size() ), xOld( fRanges.size() );
257  t = fMinTemperature;
258  deltaY = cold = 0.0;
260  for (UInt_t rIter = 0; rIter < x.size(); rIter++)
261  x[rIter] = ( fRanges[rIter]->GetMax() + fRanges[rIter]->GetMin() ) / 2.0;
262  y = yBest = 1E10;
263  for (Int_t i=0; i<fMaxCalls/50; i++) {
264  if ((i>0) && (deltaY>0.0)) {
265  cold = deltaY;
266  stopper = 1;
267  }
268  t += dT*i;
269  x = xOld = GenerateNeighbour(x,t);
270  y = yOld = fFitterTarget.EstimatorFunction( xOld );
271  equilibrium = 0;
272  for ( Int_t k=0; (k<30) && (equilibrium<=12); k++ ) {
273  xNew = GenerateNeighbour(x,t);
274  //"energy"
275  yNew = fFitterTarget.EstimatorFunction( xNew );
276  deltaY = yNew - y;
277  if (deltaY < 0.0) { // keep xnew if energy is reduced
278  std::swap(x,xNew);
279  std::swap(y,yNew);
280  if (y < yBest) {
281  xBest = x;
282  yBest = y;
283  }
284  delta = TMath::Abs( deltaY );
285  if (y != 0.0) delta /= y;
286  else if (yNew != 0.0) delta /= y;
287 
288  // equilibrium is defined as a 10% or smaller change in 10 iterations
289  if (delta < 0.1) equilibrium++;
290  else equilibrium = 0;
291  }
292  else equilibrium++;
293  }
294 
295  // "energy"
296  yNew = fFitterTarget.EstimatorFunction( xNew );
297  deltaY = yNew - yOld;
298  if ( (deltaY < 0.0 )&&( yNew < yBest)) {
299  xBest=x;
300  yBest = yNew;
301  }
302  y = yNew;
303  if ((stopper) && (deltaY >= (100.0 * cold))) break; // phase transition with another parameter to change
304  }
305  parameters = xBest;
306  return t;
307 }
308 
309 ////////////////////////////////////////////////////////////////////////////////
310 /// minimisation algorithm
311 
312 Double_t TMVA::SimulatedAnnealing::Minimize( std::vector<Double_t>& parameters )
313 {
314  std::vector<Double_t> bestParameters(fRanges.size());
315  std::vector<Double_t> oldParameters (fRanges.size());
316 
317  Double_t currentTemperature, bestFit, currentFit;
318  Int_t optimizeCalls, generalCalls, equals;
319 
320  equals = 0;
321 
324  fMinTemperature = currentTemperature = 1e-06;
325  FillWithRandomValues( parameters );
326  }
327  else fInitialTemperature = currentTemperature = GenerateMaxTemperature( parameters );
328  }
329  else {
331  currentTemperature = fMinTemperature;
332  else
333  currentTemperature = fInitialTemperature;
334  FillWithRandomValues( parameters );
335  }
336 
338 
339  Log() << kINFO
340  << "Temperatur scale = " << fTemperatureScale
341  << ", current temperature = " << currentTemperature << Endl;
342 
343  bestParameters = parameters;
344  bestFit = currentFit = fFitterTarget.EstimatorFunction( bestParameters );
345 
346  optimizeCalls = fMaxCalls/100; //use 1% calls to optimize best founded minimum
347  generalCalls = fMaxCalls - optimizeCalls; //and 99% calls to found that one
348  fProgress = 0.0;
349 
350  Timer timer( fMaxCalls, fLogger->GetSource().c_str() );
351 
352  for (Int_t sample = 0; sample < generalCalls; sample++) {
353  if (fIPyCurrentIter) *fIPyCurrentIter = sample;
354  if (fExitFromTraining && *fExitFromTraining) break;
355  GenerateNeighbour( parameters, oldParameters, currentTemperature );
356  Double_t localFit = fFitterTarget.EstimatorFunction( parameters );
357 
358  if (localFit < currentFit || TMath::Abs(currentFit-localFit) < fEps) { // if not worse than last one
359  if (TMath::Abs(currentFit-localFit) < fEps) { // if the same as last one
360  equals++;
361  if (equals >= 3) //if we still at the same level, we should increase temperature
362  fProgress+=1.0;
363  }
364  else {
365  fProgress = 0.0;
366  equals = 0;
367  }
368 
369  currentFit = localFit;
370 
371  if (currentFit < bestFit) {
372  ReWriteParameters( parameters, bestParameters );
373  bestFit = currentFit;
374  }
375  }
376  else {
377  if (!ShouldGoIn(localFit, currentFit, currentTemperature))
378  ReWriteParameters( oldParameters, parameters );
379  else
380  currentFit = localFit;
381 
382  fProgress+=1.0;
383  equals = 0;
384  }
385 
386  GenerateNewTemperature( currentTemperature, sample );
387 
388  if ((fMaxCalls<100) || sample%Int_t(fMaxCalls/100.0) == 0) timer.DrawProgressBar( sample );
389  }
390 
391  // get elapsed time
392  Log() << kINFO << "Elapsed time: " << timer.GetElapsedTime()
393  << " " << Endl;
394 
395  // suppose this minimum is the best one, now just try to improve it
396 
397  Double_t startingTemperature = fMinTemperature*(fRanges.size())*2.0;
398  currentTemperature = startingTemperature;
399 
400  Int_t changes = 0;
401  for (Int_t sample=0;sample<optimizeCalls;sample++) {
402  GenerateNeighbour( parameters, oldParameters, currentTemperature );
403  Double_t localFit = fFitterTarget.EstimatorFunction( parameters );
404 
405  if (localFit < currentFit) { //if better than last one
406  currentFit = localFit;
407  changes++;
408 
409  if (currentFit < bestFit) {
410  ReWriteParameters( parameters, bestParameters );
411  bestFit = currentFit;
412  }
413  }
414  else ReWriteParameters( oldParameters, parameters ); //we never try worse parameters
415 
416  currentTemperature-=(startingTemperature - fEps)/optimizeCalls;
417  }
418 
419  ReWriteParameters( bestParameters, parameters );
420 
421  return bestFit;
422 }
Double_t GenerateMaxTemperature(std::vector< Double_t > &parameters)
maximum temperature
void FillWithRandomValues(std::vector< Double_t > &parameters)
random starting parameters
virtual Double_t EstimatorFunction(std::vector< Double_t > &parameters)=0
void GenerateNewTemperature(Double_t &currentTemperature, Int_t Iter)
generate new temperature
Random number generator class based on M.
Definition: TRandom3.h:27
MsgLogger & Endl(MsgLogger &ml)
Definition: MsgLogger.h:158
const std::vector< TMVA::Interval * > & fRanges
Double_t Log(Double_t x)
Definition: TMath.h:759
void swap(TDirectoryEntry &e1, TDirectoryEntry &e2) noexcept
SimulatedAnnealing(IFitterTarget &target, const std::vector< TMVA::Interval *> &ranges)
constructor
Basic string class.
Definition: TString.h:131
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
Short_t Abs(Short_t d)
Definition: TMathBase.h:108
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:734
Double_t x[n]
Definition: legend1.C:17
void GenerateNeighbour(std::vector< Double_t > &parameters, std::vector< Double_t > &oldParameters, Double_t currentTemperature)
generate adjacent parameters
virtual ~SimulatedAnnealing()
destructor
void SetOptions(Int_t maxCalls, Double_t initialTemperature, Double_t minTemperature, Double_t eps, TString kernelTemperatureS, Double_t temperatureScale, Double_t adaptiveSpeed, Double_t temperatureAdaptiveStep, Bool_t useDefaultScale, Bool_t useDefaultTemperature)
option setter
unsigned int UInt_t
Definition: RtypesCore.h:42
Base implementation of simulated annealing fitting procedure.
const Bool_t kFALSE
Definition: RtypesCore.h:88
Double_t Exp(Double_t x)
Definition: TMath.h:726
#define ClassImp(name)
Definition: Rtypes.h:359
void SetDefaultScale()
setting of default scale
double Double_t
Definition: RtypesCore.h:55
Double_t y[n]
Definition: legend1.C:17
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition: TRandom.cxx:627
ostringstream derivative to redirect and format output
Definition: MsgLogger.h:59
Bool_t ShouldGoIn(Double_t currentFit, Double_t localFit, Double_t currentTemperature)
result checker
enum TMVA::SimulatedAnnealing::EKernelTemperature fKernelTemperature
std::string GetSource() const
Definition: MsgLogger.h:73
Double_t Minimize(std::vector< Double_t > &parameters)
minimisation algorithm
Double_t Sin(Double_t)
Definition: TMath.h:636
Double_t Sqrt(Double_t x)
Definition: TMath.h:690
Interface for a fitter &#39;target&#39;.
Definition: IFitterTarget.h:44
const Bool_t kTRUE
Definition: RtypesCore.h:87
Timing information for training and evaluation of MVA methods.
Definition: Timer.h:58
void ReWriteParameters(std::vector< Double_t > &from, std::vector< Double_t > &to)
copy parameters