00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #include "TMVA/SimulatedAnnealing.h"
00033
00034 #include "TRandom3.h"
00035 #include "TMath.h"
00036
00037 #include "TMVA/Interval.h"
00038 #include "TMVA/IFitterTarget.h"
00039 #include "TMVA/GeneticRange.h"
00040 #include "TMVA/Timer.h"
00041 #include "TMVA/MsgLogger.h"
00042
00043 ClassImp(TMVA::SimulatedAnnealing)
00044
00045
00046 TMVA::SimulatedAnnealing::SimulatedAnnealing( IFitterTarget& target, const std::vector<Interval*>& ranges )
00047 : fKernelTemperature (kIncreasingAdaptive),
00048 fFitterTarget ( target ),
00049 fRandom ( new TRandom3(100) ),
00050 fRanges ( ranges ),
00051 fMaxCalls ( 500000 ),
00052 fInitialTemperature ( 1000 ),
00053 fMinTemperature ( 0 ),
00054 fEps ( 1e-10 ),
00055 fTemperatureScale ( 0.06 ),
00056 fAdaptiveSpeed ( 1.0 ),
00057 fTemperatureAdaptiveStep( 0.0 ),
00058 fUseDefaultScale ( kFALSE ),
00059 fUseDefaultTemperature ( kFALSE ),
00060 fLogger( new MsgLogger("SimulatedAnnealing") ),
00061 fProgress(0.0)
00062 {
00063
00064 fKernelTemperature = kIncreasingAdaptive;
00065 }
00066
00067
00068 void TMVA::SimulatedAnnealing::SetOptions( Int_t maxCalls,
00069 Double_t initialTemperature,
00070 Double_t minTemperature,
00071 Double_t eps,
00072 TString kernelTemperatureS,
00073 Double_t temperatureScale,
00074 Double_t adaptiveSpeed,
00075 Double_t temperatureAdaptiveStep,
00076 Bool_t useDefaultScale,
00077 Bool_t useDefaultTemperature)
00078 {
00079
00080
00081 fMaxCalls = maxCalls;
00082 fInitialTemperature = initialTemperature;
00083 fMinTemperature = minTemperature;
00084 fEps = eps;
00085
00086 if (kernelTemperatureS == "IncreasingAdaptive") {
00087 fKernelTemperature = kIncreasingAdaptive;
00088 Log() << kINFO << "Using increasing adaptive algorithm" << Endl;
00089 }
00090 else if (kernelTemperatureS == "DecreasingAdaptive") {
00091 fKernelTemperature = kDecreasingAdaptive;
00092 Log() << kINFO << "Using decreasing adaptive algorithm" << Endl;
00093 }
00094 else if (kernelTemperatureS == "Sqrt") {
00095 fKernelTemperature = kSqrt;
00096 Log() << kINFO << "Using \"Sqrt\" algorithm" << Endl;
00097 }
00098 else if (kernelTemperatureS == "Homo") {
00099 fKernelTemperature = kHomo;
00100 Log() << kINFO << "Using \"Homo\" algorithm" << Endl;
00101 }
00102 else if (kernelTemperatureS == "Log") {
00103 fKernelTemperature = kLog;
00104 Log() << kINFO << "Using \"Log\" algorithm" << Endl;
00105 }
00106 else if (kernelTemperatureS == "Sin") {
00107 fKernelTemperature = kSin;
00108 Log() << kINFO << "Using \"Sin\" algorithm" << Endl;
00109 }
00110
00111 fTemperatureScale = temperatureScale;
00112 fAdaptiveSpeed = adaptiveSpeed;
00113 fTemperatureAdaptiveStep = temperatureAdaptiveStep;
00114
00115 fUseDefaultScale = useDefaultScale;
00116 fUseDefaultTemperature = useDefaultTemperature;
00117 }
00118
00119 TMVA::SimulatedAnnealing::~SimulatedAnnealing()
00120 {
00121
00122 }
00123
00124
00125 void TMVA::SimulatedAnnealing::FillWithRandomValues( std::vector<Double_t>& parameters )
00126 {
00127
00128 for (UInt_t rIter = 0; rIter < parameters.size(); rIter++) {
00129 parameters[rIter] = fRandom->Uniform(0.0,1.0)*(fRanges[rIter]->GetMax() - fRanges[rIter]->GetMin()) + fRanges[rIter]->GetMin();
00130 }
00131 }
00132
00133
00134 void TMVA::SimulatedAnnealing::ReWriteParameters( std::vector<Double_t>& from, std::vector<Double_t>& to)
00135 {
00136
00137 for (UInt_t rIter = 0; rIter < from.size(); rIter++) to[rIter] = from[rIter];
00138 }
00139
00140
00141 void TMVA::SimulatedAnnealing::GenerateNeighbour( std::vector<Double_t>& parameters, std::vector<Double_t>& oldParameters,
00142 Double_t currentTemperature )
00143 {
00144
00145 ReWriteParameters( parameters, oldParameters );
00146
00147 for (UInt_t rIter=0;rIter<parameters.size();rIter++) {
00148 Double_t uni,distribution,sign;
00149 do {
00150 uni = fRandom->Uniform(0.0,1.0);
00151 sign = (uni - 0.5 >= 0.0) ? (1.0) : (-1.0);
00152 distribution = currentTemperature * (TMath::Power(1.0 + 1.0/currentTemperature, TMath::Abs(2.0*uni - 1.0)) -1.0)*sign;
00153 parameters[rIter] = oldParameters[rIter] + (fRanges[rIter]->GetMax()-fRanges[rIter]->GetMin())*0.1*distribution;
00154 }
00155 while (parameters[rIter] < fRanges[rIter]->GetMin() || parameters[rIter] > fRanges[rIter]->GetMax() );
00156 }
00157 }
00158
00159 std::vector<Double_t> TMVA::SimulatedAnnealing::GenerateNeighbour( std::vector<Double_t>& parameters, Double_t currentTemperature )
00160 {
00161
00162 std::vector<Double_t> newParameters( fRanges.size() );
00163
00164 for (UInt_t rIter=0; rIter<parameters.size(); rIter++) {
00165 Double_t uni,distribution,sign;
00166 do {
00167 uni = fRandom->Uniform(0.0,1.0);
00168 sign = (uni - 0.5 >= 0.0) ? (1.0) : (-1.0);
00169 distribution = currentTemperature * (TMath::Power(1.0 + 1.0/currentTemperature, TMath::Abs(2.0*uni - 1.0)) -1.0)*sign;
00170 newParameters[rIter] = parameters[rIter] + (fRanges[rIter]->GetMax()-fRanges[rIter]->GetMin())*0.1*distribution;
00171 }
00172 while (newParameters[rIter] < fRanges[rIter]->GetMin() || newParameters[rIter] > fRanges[rIter]->GetMax() );
00173 }
00174
00175 return newParameters;
00176 }
00177
00178
00179 void TMVA::SimulatedAnnealing::GenerateNewTemperature( Double_t& currentTemperature, Int_t Iter )
00180 {
00181
00182 if (fKernelTemperature == kSqrt) {
00183 currentTemperature = fInitialTemperature/(Double_t)TMath::Sqrt(Iter+2) * fTemperatureScale;
00184 }
00185 else if (fKernelTemperature == kLog) {
00186 currentTemperature = fInitialTemperature/(Double_t)TMath::Log(Iter+2) * fTemperatureScale;
00187 }
00188 else if (fKernelTemperature == kHomo) {
00189 currentTemperature = fInitialTemperature/(Double_t)(Iter+2) * fTemperatureScale;
00190 }
00191 else if (fKernelTemperature == kSin) {
00192 currentTemperature = (TMath::Sin( (Double_t)Iter / fTemperatureScale ) + 1.0 )/ (Double_t)(Iter+1.0) * fInitialTemperature + fEps;
00193 }
00194 else if (fKernelTemperature == kGeo) {
00195 currentTemperature = currentTemperature*fTemperatureScale;
00196 }
00197 else if (fKernelTemperature == kIncreasingAdaptive) {
00198 currentTemperature = fMinTemperature + fTemperatureScale*TMath::Log(1.0+fProgress*fAdaptiveSpeed);
00199 }
00200 else if (fKernelTemperature == kDecreasingAdaptive) {
00201 currentTemperature = currentTemperature*fTemperatureScale;
00202 }
00203 else Log() << kFATAL << "No such kernel!" << Endl;
00204 }
00205
00206
00207 Bool_t TMVA::SimulatedAnnealing::ShouldGoIn( Double_t currentFit, Double_t localFit, Double_t currentTemperature )
00208 {
00209
00210 if (currentTemperature < fEps) return kFALSE;
00211 Double_t lim = TMath::Exp( -TMath::Abs( currentFit - localFit ) / currentTemperature );
00212 Double_t prob = fRandom->Uniform(0.0, 1.0);
00213 return (prob < lim) ? kTRUE : kFALSE;
00214 }
00215
00216
00217 void TMVA::SimulatedAnnealing::SetDefaultScale()
00218 {
00219
00220 if (fKernelTemperature == kSqrt) fTemperatureScale = 1.0;
00221 else if (fKernelTemperature == kLog) fTemperatureScale = 1.0;
00222 else if (fKernelTemperature == kHomo) fTemperatureScale = 1.0;
00223 else if (fKernelTemperature == kSin) fTemperatureScale = 20.0;
00224 else if (fKernelTemperature == kGeo) fTemperatureScale = 0.99997;
00225 else if (fKernelTemperature == kDecreasingAdaptive) {
00226 fTemperatureScale = 1.0;
00227 while (TMath::Abs(TMath::Power(fTemperatureScale,fMaxCalls) * fInitialTemperature - fMinTemperature) >
00228 TMath::Abs(TMath::Power(fTemperatureScale-0.000001,fMaxCalls) * fInitialTemperature - fMinTemperature)) {
00229 fTemperatureScale -= 0.000001;
00230 }
00231 }
00232 else if (fKernelTemperature == kIncreasingAdaptive) fTemperatureScale = 0.15*( 1.0 / (Double_t)(fRanges.size() ) );
00233 else Log() << kFATAL << "No such kernel!" << Endl;
00234 }
00235
00236
00237 Double_t TMVA::SimulatedAnnealing::GenerateMaxTemperature( std::vector<Double_t>& parameters )
00238 {
00239
00240 Int_t equilibrium;
00241 Bool_t stopper = 0;
00242 Double_t t, dT, cold, delta, deltaY, y, yNew, yBest, yOld;
00243 std::vector<Double_t> x( fRanges.size() ), xNew( fRanges.size() ), xBest( fRanges.size() ), xOld( fRanges.size() );
00244 t = fMinTemperature;
00245 deltaY = cold = 0.0;
00246 dT = fTemperatureAdaptiveStep;
00247 for (UInt_t rIter = 0; rIter < x.size(); rIter++)
00248 x[rIter] = ( fRanges[rIter]->GetMax() + fRanges[rIter]->GetMin() ) / 2.0;
00249 y = yBest = 1E10;
00250 for (Int_t i=0; i<fMaxCalls/50; i++) {
00251 if ((i>0) && (deltaY>0.0)) {
00252 cold = deltaY;
00253 stopper = 1;
00254 }
00255 t += dT*i;
00256 x = xOld = GenerateNeighbour(x,t);
00257 y = yOld = fFitterTarget.EstimatorFunction( xOld );
00258 equilibrium = 0;
00259 for ( Int_t k=0; (k<30) && (equilibrium<=12); k++ ) {
00260 xNew = GenerateNeighbour(x,t);
00261
00262 yNew = fFitterTarget.EstimatorFunction( xNew );
00263 deltaY = yNew - y;
00264 if (deltaY < 0.0) {
00265 std::swap(x,xNew);
00266 std::swap(y,yNew);
00267 if (y < yBest) {
00268 xBest = x;
00269 yBest = y;
00270 }
00271 delta = TMath::Abs( deltaY );
00272 if (y != 0.0) delta /= y;
00273 else if (yNew != 0.0) delta /= y;
00274
00275
00276 if (delta < 0.1) equilibrium++;
00277 else equilibrium = 0;
00278 }
00279 else equilibrium++;
00280 }
00281
00282
00283 yNew = fFitterTarget.EstimatorFunction( xNew );
00284 deltaY = yNew - yOld;
00285 if ( (deltaY < 0.0 )&&( yNew < yBest)) {
00286 xBest=x;
00287 yBest = yNew;
00288 }
00289 y = yNew;
00290 if ((stopper) && (deltaY >= (100.0 * cold))) break;
00291 }
00292 parameters = xBest;
00293 return t;
00294 }
00295
00296
00297 Double_t TMVA::SimulatedAnnealing::Minimize( std::vector<Double_t>& parameters )
00298 {
00299
00300 std::vector<Double_t> bestParameters(fRanges.size());
00301 std::vector<Double_t> oldParameters (fRanges.size());
00302
00303 Double_t currentTemperature, bestFit, currentFit;
00304 Int_t optimizeCalls, generalCalls, equals;
00305
00306 equals = 0;
00307
00308 if (fUseDefaultTemperature) {
00309 if (fKernelTemperature == kIncreasingAdaptive) {
00310 fMinTemperature = currentTemperature = 1e-06;
00311 FillWithRandomValues( parameters );
00312 }
00313 else fInitialTemperature = currentTemperature = GenerateMaxTemperature( parameters );
00314 }
00315 else {
00316 if (fKernelTemperature == kIncreasingAdaptive)
00317 currentTemperature = fMinTemperature;
00318 else
00319 currentTemperature = fInitialTemperature;
00320 FillWithRandomValues( parameters );
00321 }
00322
00323 if (fUseDefaultScale) SetDefaultScale();
00324
00325 Log() << kINFO
00326 << "Temperatur scale = " << fTemperatureScale
00327 << ", current temperature = " << currentTemperature << Endl;
00328
00329 bestParameters = parameters;
00330 bestFit = currentFit = fFitterTarget.EstimatorFunction( bestParameters );
00331
00332 optimizeCalls = fMaxCalls/100;
00333 generalCalls = fMaxCalls - optimizeCalls;
00334 fProgress = 0.0;
00335
00336 Timer timer( fMaxCalls, fLogger->GetSource().c_str() );
00337
00338 for (Int_t sample = 0; sample < generalCalls; sample++) {
00339 GenerateNeighbour( parameters, oldParameters, currentTemperature );
00340 Double_t localFit = fFitterTarget.EstimatorFunction( parameters );
00341
00342 if (localFit < currentFit || TMath::Abs(currentFit-localFit) < fEps) {
00343 if (TMath::Abs(currentFit-localFit) < fEps) {
00344 equals++;
00345 if (equals >= 3)
00346 fProgress+=1.0;
00347 }
00348 else {
00349 fProgress = 0.0;
00350 equals = 0;
00351 }
00352
00353 currentFit = localFit;
00354
00355 if (currentFit < bestFit) {
00356 ReWriteParameters( parameters, bestParameters );
00357 bestFit = currentFit;
00358 }
00359 }
00360 else {
00361 if (!ShouldGoIn(localFit, currentFit, currentTemperature))
00362 ReWriteParameters( oldParameters, parameters );
00363 else
00364 currentFit = localFit;
00365
00366 fProgress+=1.0;
00367 equals = 0;
00368 }
00369
00370 GenerateNewTemperature( currentTemperature, sample );
00371
00372 if ((fMaxCalls<100) || sample%Int_t(fMaxCalls/100.0) == 0) timer.DrawProgressBar( sample );
00373 }
00374
00375
00376 Log() << kINFO << "Elapsed time: " << timer.GetElapsedTime()
00377 << " " << Endl;
00378
00379
00380
00381 Double_t startingTemperature = fMinTemperature*(fRanges.size())*2.0;
00382 currentTemperature = startingTemperature;
00383
00384 Int_t changes = 0;
00385 for (Int_t sample=0;sample<optimizeCalls;sample++) {
00386 GenerateNeighbour( parameters, oldParameters, currentTemperature );
00387 Double_t localFit = fFitterTarget.EstimatorFunction( parameters );
00388
00389 if (localFit < currentFit) {
00390 currentFit = localFit;
00391 changes++;
00392
00393 if (currentFit < bestFit) {
00394 ReWriteParameters( parameters, bestParameters );
00395 bestFit = currentFit;
00396 }
00397 }
00398 else ReWriteParameters( oldParameters, parameters );
00399
00400 currentTemperature-=(startingTemperature - fEps)/optimizeCalls;
00401 }
00402
00403 ReWriteParameters( bestParameters, parameters );
00404
00405 return bestFit;
00406 }
00407
