Logo ROOT  
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
30Base 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
50TMVA::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") {
90 fKernelTemperature = kIncreasingAdaptive;
91 Log() << kINFO << "Using increasing adaptive algorithm" << Endl;
92 }
93 else if (kernelTemperatureS == "DecreasingAdaptive") {
94 fKernelTemperature = kDecreasingAdaptive;
95 Log() << kINFO << "Using decreasing adaptive algorithm" << Endl;
96 }
97 else if (kernelTemperatureS == "Sqrt") {
98 fKernelTemperature = kSqrt;
99 Log() << kINFO << "Using \"Sqrt\" algorithm" << Endl;
100 }
101 else if (kernelTemperatureS == "Homo") {
102 fKernelTemperature = kHomo;
103 Log() << kINFO << "Using \"Homo\" algorithm" << Endl;
104 }
105 else if (kernelTemperatureS == "Log") {
106 fKernelTemperature = kLog;
107 Log() << kINFO << "Using \"Log\" algorithm" << Endl;
108 }
109 else if (kernelTemperatureS == "Sin") {
110 fKernelTemperature = kSin;
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
132void 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
142void 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
150void 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
169std::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 }
207 else if (fKernelTemperature == kIncreasingAdaptive) {
208 currentTemperature = fMinTemperature + fTemperatureScale*TMath::Log(1.0+fProgress*fAdaptiveSpeed);
209 }
210 else if (fKernelTemperature == kDecreasingAdaptive) {
211 currentTemperature = currentTemperature*fTemperatureScale;
212 }
213 else Log() << kFATAL << "No such kernel!" << Endl;
214}
215
216////////////////////////////////////////////////////////////////////////////////
217/// result checker
218
219Bool_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{
232 if (fKernelTemperature == kSqrt) fTemperatureScale = 1.0;
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;
237 else if (fKernelTemperature == kDecreasingAdaptive) {
238 fTemperatureScale = 1.0;
239 while (TMath::Abs(TMath::Power(fTemperatureScale,fMaxCalls) * fInitialTemperature - fMinTemperature) >
240 TMath::Abs(TMath::Power(fTemperatureScale-0.000001,fMaxCalls) * fInitialTemperature - fMinTemperature)) {
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
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;
259 dT = fTemperatureAdaptiveStep;
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
312Double_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
322 if (fUseDefaultTemperature) {
323 if (fKernelTemperature == kIncreasingAdaptive) {
324 fMinTemperature = currentTemperature = 1e-06;
325 FillWithRandomValues( parameters );
326 }
327 else fInitialTemperature = currentTemperature = GenerateMaxTemperature( parameters );
328 }
329 else {
330 if (fKernelTemperature == kIncreasingAdaptive)
331 currentTemperature = fMinTemperature;
332 else
333 currentTemperature = fInitialTemperature;
334 FillWithRandomValues( parameters );
335 }
336
337 if (fUseDefaultScale) SetDefaultScale();
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}
#define e(i)
Definition: RSha256.hxx:103
int Int_t
Definition: RtypesCore.h:45
unsigned int UInt_t
Definition: RtypesCore.h:46
const Bool_t kFALSE
Definition: RtypesCore.h:101
bool Bool_t
Definition: RtypesCore.h:63
double Double_t
Definition: RtypesCore.h:59
const Bool_t kTRUE
Definition: RtypesCore.h:100
#define ClassImp(name)
Definition: Rtypes.h:364
Interface for a fitter 'target'.
Definition: IFitterTarget.h:44
ostringstream derivative to redirect and format output
Definition: MsgLogger.h:57
Base implementation of simulated annealing fitting procedure.
void GenerateNeighbour(std::vector< Double_t > &parameters, std::vector< Double_t > &oldParameters, Double_t currentTemperature)
generate adjacent parameters
Double_t Minimize(std::vector< Double_t > &parameters)
minimisation algorithm
virtual ~SimulatedAnnealing()
destructor
Double_t GenerateMaxTemperature(std::vector< Double_t > &parameters)
maximum temperature
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
void SetDefaultScale()
setting of default scale
void GenerateNewTemperature(Double_t &currentTemperature, Int_t Iter)
generate new temperature
void FillWithRandomValues(std::vector< Double_t > &parameters)
random starting parameters
void ReWriteParameters(std::vector< Double_t > &from, std::vector< Double_t > &to)
copy parameters
enum TMVA::SimulatedAnnealing::EKernelTemperature fKernelTemperature
SimulatedAnnealing(IFitterTarget &target, const std::vector< TMVA::Interval * > &ranges)
constructor
Bool_t ShouldGoIn(Double_t currentFit, Double_t localFit, Double_t currentTemperature)
result checker
Timing information for training and evaluation of MVA methods.
Definition: Timer.h:58
TString GetElapsedTime(Bool_t Scientific=kTRUE)
returns pretty string with elapsed time
Definition: Timer.cxx:146
void DrawProgressBar(Int_t, const TString &comment="")
draws progress bar in color or B&W caution:
Definition: Timer.cxx:202
@ kINFO
Definition: Types.h:58
@ kFATAL
Definition: Types.h:61
Random number generator class based on M.
Definition: TRandom3.h:27
Basic string class.
Definition: TString.h:136
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
void swap(RDirectoryEntry &e1, RDirectoryEntry &e2) noexcept
MsgLogger & Endl(MsgLogger &ml)
Definition: MsgLogger.h:148
Double_t Exp(Double_t x)
Definition: TMath.h:677
Double_t Log(Double_t x)
Definition: TMath.h:710
Double_t Sqrt(Double_t x)
Definition: TMath.h:641
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:685
Double_t Sin(Double_t)
Definition: TMath.h:589
Short_t Abs(Short_t d)
Definition: TMathBase.h:120