Logo ROOT   6.08/07
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 //_______________________________________________________________________
29 //
30 // Implementation of Simulated Annealing fitter
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 {
67  fKernelTemperature = kIncreasingAdaptive;
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 /// destructor
123 
125 {
126 }
127 
128 ////////////////////////////////////////////////////////////////////////////////
129 /// random starting parameters
130 
131 void TMVA::SimulatedAnnealing::FillWithRandomValues( std::vector<Double_t>& parameters )
132 {
133  for (UInt_t rIter = 0; rIter < parameters.size(); rIter++) {
134  parameters[rIter] = fRandom->Uniform(0.0,1.0)*(fRanges[rIter]->GetMax() - fRanges[rIter]->GetMin()) + fRanges[rIter]->GetMin();
135  }
136 }
137 
138 ////////////////////////////////////////////////////////////////////////////////
139 /// copy parameters
140 
141 void TMVA::SimulatedAnnealing::ReWriteParameters( std::vector<Double_t>& from, std::vector<Double_t>& to)
142 {
143  for (UInt_t rIter = 0; rIter < from.size(); rIter++) to[rIter] = from[rIter];
144 }
145 
146 ////////////////////////////////////////////////////////////////////////////////
147 /// generate adjacent parameters
148 
149 void TMVA::SimulatedAnnealing::GenerateNeighbour( std::vector<Double_t>& parameters, std::vector<Double_t>& oldParameters,
150  Double_t currentTemperature )
151 {
152  ReWriteParameters( parameters, oldParameters );
153 
154  for (UInt_t rIter=0;rIter<parameters.size();rIter++) {
155  Double_t uni,distribution,sign;
156  do {
157  uni = fRandom->Uniform(0.0,1.0);
158  sign = (uni - 0.5 >= 0.0) ? (1.0) : (-1.0);
159  distribution = currentTemperature * (TMath::Power(1.0 + 1.0/currentTemperature, TMath::Abs(2.0*uni - 1.0)) -1.0)*sign;
160  parameters[rIter] = oldParameters[rIter] + (fRanges[rIter]->GetMax()-fRanges[rIter]->GetMin())*0.1*distribution;
161  }
162  while (parameters[rIter] < fRanges[rIter]->GetMin() || parameters[rIter] > fRanges[rIter]->GetMax() );
163  }
164 }
165 ////////////////////////////////////////////////////////////////////////////////
166 /// generate adjacent parameters
167 
168 std::vector<Double_t> TMVA::SimulatedAnnealing::GenerateNeighbour( std::vector<Double_t>& parameters, Double_t currentTemperature )
169 {
170  std::vector<Double_t> newParameters( fRanges.size() );
171 
172  for (UInt_t rIter=0; rIter<parameters.size(); rIter++) {
173  Double_t uni,distribution,sign;
174  do {
175  uni = fRandom->Uniform(0.0,1.0);
176  sign = (uni - 0.5 >= 0.0) ? (1.0) : (-1.0);
177  distribution = currentTemperature * (TMath::Power(1.0 + 1.0/currentTemperature, TMath::Abs(2.0*uni - 1.0)) -1.0)*sign;
178  newParameters[rIter] = parameters[rIter] + (fRanges[rIter]->GetMax()-fRanges[rIter]->GetMin())*0.1*distribution;
179  }
180  while (newParameters[rIter] < fRanges[rIter]->GetMin() || newParameters[rIter] > fRanges[rIter]->GetMax() );
181  }
182 
183  return newParameters;
184 }
185 
186 ////////////////////////////////////////////////////////////////////////////////
187 /// generate new temperature
188 
190 {
191  if (fKernelTemperature == kSqrt) {
192  currentTemperature = fInitialTemperature/(Double_t)TMath::Sqrt(Iter+2) * fTemperatureScale;
193  }
194  else if (fKernelTemperature == kLog) {
195  currentTemperature = fInitialTemperature/(Double_t)TMath::Log(Iter+2) * fTemperatureScale;
196  }
197  else if (fKernelTemperature == kHomo) {
198  currentTemperature = fInitialTemperature/(Double_t)(Iter+2) * fTemperatureScale;
199  }
200  else if (fKernelTemperature == kSin) {
201  currentTemperature = (TMath::Sin( (Double_t)Iter / fTemperatureScale ) + 1.0 )/ (Double_t)(Iter+1.0) * fInitialTemperature + fEps;
202  }
203  else if (fKernelTemperature == kGeo) {
204  currentTemperature = currentTemperature*fTemperatureScale;
205  }
208  }
210  currentTemperature = currentTemperature*fTemperatureScale;
211  }
212  else Log() << kFATAL << "No such kernel!" << Endl;
213 }
214 
215 ////////////////////////////////////////////////////////////////////////////////
216 /// result checker
217 
218 Bool_t TMVA::SimulatedAnnealing::ShouldGoIn( Double_t currentFit, Double_t localFit, Double_t currentTemperature )
219 {
220  if (currentTemperature < fEps) return kFALSE;
221  Double_t lim = TMath::Exp( -TMath::Abs( currentFit - localFit ) / currentTemperature );
222  Double_t prob = fRandom->Uniform(0.0, 1.0);
223  return (prob < lim) ? kTRUE : kFALSE;
224 }
225 
226 ////////////////////////////////////////////////////////////////////////////////
227 /// setting of default scale
228 
230 {
232  else if (fKernelTemperature == kLog) fTemperatureScale = 1.0;
233  else if (fKernelTemperature == kHomo) fTemperatureScale = 1.0;
234  else if (fKernelTemperature == kSin) fTemperatureScale = 20.0;
235  else if (fKernelTemperature == kGeo) fTemperatureScale = 0.99997;
237  fTemperatureScale = 1.0;
240  fTemperatureScale -= 0.000001;
241  }
242  }
243  else if (fKernelTemperature == kIncreasingAdaptive) fTemperatureScale = 0.15*( 1.0 / (Double_t)(fRanges.size() ) );
244  else Log() << kFATAL << "No such kernel!" << Endl;
245 }
246 
247 ////////////////////////////////////////////////////////////////////////////////
248 /// maximum temperature
249 
250 Double_t TMVA::SimulatedAnnealing::GenerateMaxTemperature( std::vector<Double_t>& parameters )
251 {
252  Int_t equilibrium;
253  Bool_t stopper = 0;
254  Double_t t, dT, cold, delta, deltaY, y, yNew, yBest, yOld;
255  std::vector<Double_t> x( fRanges.size() ), xNew( fRanges.size() ), xBest( fRanges.size() ), xOld( fRanges.size() );
256  t = fMinTemperature;
257  deltaY = cold = 0.0;
259  for (UInt_t rIter = 0; rIter < x.size(); rIter++)
260  x[rIter] = ( fRanges[rIter]->GetMax() + fRanges[rIter]->GetMin() ) / 2.0;
261  y = yBest = 1E10;
262  for (Int_t i=0; i<fMaxCalls/50; i++) {
263  if ((i>0) && (deltaY>0.0)) {
264  cold = deltaY;
265  stopper = 1;
266  }
267  t += dT*i;
268  x = xOld = GenerateNeighbour(x,t);
269  y = yOld = fFitterTarget.EstimatorFunction( xOld );
270  equilibrium = 0;
271  for ( Int_t k=0; (k<30) && (equilibrium<=12); k++ ) {
272  xNew = GenerateNeighbour(x,t);
273  //"energy"
274  yNew = fFitterTarget.EstimatorFunction( xNew );
275  deltaY = yNew - y;
276  if (deltaY < 0.0) { // keep xnew if energy is reduced
277  std::swap(x,xNew);
278  std::swap(y,yNew);
279  if (y < yBest) {
280  xBest = x;
281  yBest = y;
282  }
283  delta = TMath::Abs( deltaY );
284  if (y != 0.0) delta /= y;
285  else if (yNew != 0.0) delta /= y;
286 
287  // equilibrium is defined as a 10% or smaller change in 10 iterations
288  if (delta < 0.1) equilibrium++;
289  else equilibrium = 0;
290  }
291  else equilibrium++;
292  }
293 
294  // "energy"
295  yNew = fFitterTarget.EstimatorFunction( xNew );
296  deltaY = yNew - yOld;
297  if ( (deltaY < 0.0 )&&( yNew < yBest)) {
298  xBest=x;
299  yBest = yNew;
300  }
301  y = yNew;
302  if ((stopper) && (deltaY >= (100.0 * cold))) break; // phase transition with another parameter to change
303  }
304  parameters = xBest;
305  return t;
306 }
307 
308 ////////////////////////////////////////////////////////////////////////////////
309 /// minimisation algorithm
310 
311 Double_t TMVA::SimulatedAnnealing::Minimize( std::vector<Double_t>& parameters )
312 {
313  std::vector<Double_t> bestParameters(fRanges.size());
314  std::vector<Double_t> oldParameters (fRanges.size());
315 
316  Double_t currentTemperature, bestFit, currentFit;
317  Int_t optimizeCalls, generalCalls, equals;
318 
319  equals = 0;
320 
323  fMinTemperature = currentTemperature = 1e-06;
324  FillWithRandomValues( parameters );
325  }
326  else fInitialTemperature = currentTemperature = GenerateMaxTemperature( parameters );
327  }
328  else {
330  currentTemperature = fMinTemperature;
331  else
332  currentTemperature = fInitialTemperature;
333  FillWithRandomValues( parameters );
334  }
335 
337 
338  Log() << kINFO
339  << "Temperatur scale = " << fTemperatureScale
340  << ", current temperature = " << currentTemperature << Endl;
341 
342  bestParameters = parameters;
343  bestFit = currentFit = fFitterTarget.EstimatorFunction( bestParameters );
344 
345  optimizeCalls = fMaxCalls/100; //use 1% calls to optimize best founded minimum
346  generalCalls = fMaxCalls - optimizeCalls; //and 99% calls to found that one
347  fProgress = 0.0;
348 
349  Timer timer( fMaxCalls, fLogger->GetSource().c_str() );
350 
351  for (Int_t sample = 0; sample < generalCalls; sample++) {
352  if (fIPyCurrentIter) *fIPyCurrentIter = sample;
353  if (fExitFromTraining && *fExitFromTraining) break;
354  GenerateNeighbour( parameters, oldParameters, currentTemperature );
355  Double_t localFit = fFitterTarget.EstimatorFunction( parameters );
356 
357  if (localFit < currentFit || TMath::Abs(currentFit-localFit) < fEps) { // if not worse than last one
358  if (TMath::Abs(currentFit-localFit) < fEps) { // if the same as last one
359  equals++;
360  if (equals >= 3) //if we still at the same level, we should increase temperature
361  fProgress+=1.0;
362  }
363  else {
364  fProgress = 0.0;
365  equals = 0;
366  }
367 
368  currentFit = localFit;
369 
370  if (currentFit < bestFit) {
371  ReWriteParameters( parameters, bestParameters );
372  bestFit = currentFit;
373  }
374  }
375  else {
376  if (!ShouldGoIn(localFit, currentFit, currentTemperature))
377  ReWriteParameters( oldParameters, parameters );
378  else
379  currentFit = localFit;
380 
381  fProgress+=1.0;
382  equals = 0;
383  }
384 
385  GenerateNewTemperature( currentTemperature, sample );
386 
387  if ((fMaxCalls<100) || sample%Int_t(fMaxCalls/100.0) == 0) timer.DrawProgressBar( sample );
388  }
389 
390  // get elapsed time
391  Log() << kINFO << "Elapsed time: " << timer.GetElapsedTime()
392  << " " << Endl;
393 
394  // supose this minimum is the best one, now just try to improve it
395 
396  Double_t startingTemperature = fMinTemperature*(fRanges.size())*2.0;
397  currentTemperature = startingTemperature;
398 
399  Int_t changes = 0;
400  for (Int_t sample=0;sample<optimizeCalls;sample++) {
401  GenerateNeighbour( parameters, oldParameters, currentTemperature );
402  Double_t localFit = fFitterTarget.EstimatorFunction( parameters );
403 
404  if (localFit < currentFit) { //if better than last one
405  currentFit = localFit;
406  changes++;
407 
408  if (currentFit < bestFit) {
409  ReWriteParameters( parameters, bestParameters );
410  bestFit = currentFit;
411  }
412  }
413  else ReWriteParameters( oldParameters, parameters ); //we never try worse parameters
414 
415  currentTemperature-=(startingTemperature - fEps)/optimizeCalls;
416  }
417 
418  ReWriteParameters( bestParameters, parameters );
419 
420  return bestFit;
421 }
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:29
MsgLogger & Endl(MsgLogger &ml)
Definition: MsgLogger.h:162
const std::vector< TMVA::Interval * > & fRanges
Double_t Log(Double_t x)
Definition: TMath.h:526
void swap(TDirectoryEntry &e1, TDirectoryEntry &e2) noexcept
Basic string class.
Definition: TString.h:137
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
const Bool_t kFALSE
Definition: Rtypes.h:92
int equals(Double_t n1, Double_t n2, double ERRORLIMIT=1.E-10)
Definition: UnitTesting.cxx:24
STL namespace.
Short_t Abs(Short_t d)
Definition: TMathBase.h:110
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:501
TStopwatch timer
Definition: pirndm.C:37
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
Double_t Exp(Double_t x)
Definition: TMath.h:495
#define ClassImp(name)
Definition: Rtypes.h:279
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:606
Bool_t ShouldGoIn(Double_t currentFit, Double_t localFit, Double_t currentTemperature)
result checker
enum TMVA::SimulatedAnnealing::EKernelTemperature fKernelTemperature
Abstract ClassifierFactory template that handles arbitrary types.
std::string GetSource() const
Definition: MsgLogger.h:77
Double_t Minimize(std::vector< Double_t > &parameters)
minimisation algorithm
Double_t Sin(Double_t)
Definition: TMath.h:421
Double_t Sqrt(Double_t x)
Definition: TMath.h:464
const Bool_t kTRUE
Definition: Rtypes.h:91
void ReWriteParameters(std::vector< Double_t > &from, std::vector< Double_t > &to)
copy parameters