00001 #include "RecoVertex/VertexTools/interface/DeterministicAnnealing.h" 00002 #include <cmath> 00003 #include <vector> 00004 #include <iostream> 00005 #include <limits> 00006 00007 using namespace std; 00008 00009 namespace { 00010 vector < float > temperatures; 00011 } 00012 00013 DeterministicAnnealing::DeterministicAnnealing ( float cutoff ) : 00014 theIndex(0), theChi2cut ( cutoff*cutoff ), theIsAnnealed ( false ) 00015 { 00016 temperatures.push_back(256); 00017 temperatures.push_back(64); 00018 temperatures.push_back(16); 00019 temperatures.push_back(4); 00020 temperatures.push_back(2); 00021 temperatures.push_back(1); 00022 } 00023 00024 DeterministicAnnealing::DeterministicAnnealing( const vector < float > & sched, 00025 float cutoff ) : theIndex(0), theChi2cut ( cutoff*cutoff ), theIsAnnealed ( false ) 00026 { 00027 temperatures = sched; 00028 } 00029 00030 void DeterministicAnnealing::anneal() 00031 { 00032 if ( theIndex < ( temperatures.size() - 1 ) ) 00033 { 00034 theIndex++; 00035 } else { 00036 theIsAnnealed = true; 00037 }; 00038 } 00039 00040 double DeterministicAnnealing::weight ( double chi2 ) const 00041 { 00042 long double mphi = phi ( chi2 ); 00043 /* 00044 if ( mphi < std::numeric_limits<double>::epsilon() ) return 0.; 00045 return 1. / ( 1. + phi ( theChi2cut * theChi2cut ) / mphi ); 00046 */ 00047 // return mphi / ( mphi + phi ( theChi2cut ) ); 00048 long double newtmp = mphi / ( mphi + phi ( theChi2cut ) ); 00049 if ( !finite(newtmp ) ) 00050 { 00051 if ( chi2 < theChi2cut ) newtmp=1.; 00052 else newtmp=0.; 00053 } 00054 return newtmp; 00055 } 00056 00057 void DeterministicAnnealing::resetAnnealing() 00058 { 00059 theIndex=0; 00060 theIsAnnealed = false; 00061 } 00062 00063 inline double DeterministicAnnealing::phi( double chi2 ) const 00064 { 00065 return exp ( -.5 * chi2 / temperatures[theIndex] ); 00066 } 00067 00068 double DeterministicAnnealing::cutoff() const 00069 { 00070 return sqrt(theChi2cut); 00071 } 00072 00073 double DeterministicAnnealing::currentTemp() const 00074 { 00075 return temperatures[theIndex]; 00076 } 00077 00078 double DeterministicAnnealing::initialTemp() const 00079 { 00080 return temperatures[0]; 00081 } 00082 00083 bool DeterministicAnnealing::isAnnealed() const 00084 { 00085 return theIsAnnealed; 00086 } 00087 00088 void DeterministicAnnealing::debug() const 00089 { 00090 cout << "[DeterministicAnnealing] schedule="; 00091 for ( vector< float >::const_iterator i=temperatures.begin(); 00092 i!=temperatures.end() ; ++i ) 00093 { 00094 cout << *i << " "; 00095 }; 00096 cout << endl; 00097 }