VertexTools/src/DeterministicAnnealing.cc

0001 #include "RecoVertex/VertexTools/interface/DeterministicAnnealing.h"
0002 #include "FWCore/Utilities/interface/isFinite.h"
0003 #include <cmath>
0004 #include <vector>
0005 #include <iostream>
0006 #include <limits>
0007
0008 using namespace std;
0009
0010 DeterministicAnnealing::DeterministicAnnealing(float cutoff)
0011     : theTemperatures({256, 64, 16, 4, 2, 1}), theIndex(0), theChi2cut(cutoff * cutoff), theIsAnnealed(false) {}
0012
0013 DeterministicAnnealing::DeterministicAnnealing(const vector<float>& sched, float cutoff)
0014     : theTemperatures(sched), theIndex(0), theChi2cut(cutoff * cutoff), theIsAnnealed(false) {}
0015
0016 void DeterministicAnnealing::anneal() {
0017   if (theIndex < (theTemperatures.size() - 1)) {
0018     theIndex++;
0019   } else {
0020     theIsAnnealed = true;
0021   };
0022 }
0023
0024 double DeterministicAnnealing::weight(double chi2) const {
0025   long double mphi = phi(chi2);
0026   /*
0027   if ( mphi < std::numeric_limits<double>::epsilon() ) return 0.;
0028   return 1. / ( 1. + phi ( theChi2cut * theChi2cut ) / mphi );
0029   */
0030   // return mphi / ( mphi + phi ( theChi2cut ) );
0031   long double newtmp = mphi / (mphi + phi(theChi2cut));
0032   if (edm::isNotFinite(newtmp)) {
0033     if (chi2 < theChi2cut)
0034       newtmp = 1.;
0035     else
0036       newtmp = 0.;
0037   }
0038   return newtmp;
0039 }
0040
0041 void DeterministicAnnealing::resetAnnealing() {
0042   theIndex = 0;
0043   theIsAnnealed = false;
0044 }
0045
0046 inline double DeterministicAnnealing::phi(double chi2) const { return exp(-.5 * chi2 / theTemperatures[theIndex]); }
0047
0048 double DeterministicAnnealing::cutoff() const { return sqrt(theChi2cut); }
0049
0050 double DeterministicAnnealing::currentTemp() const { return theTemperatures[theIndex]; }
0051
0052 double DeterministicAnnealing::initialTemp() const { return theTemperatures[0]; }
0053
0054 bool DeterministicAnnealing::isAnnealed() const { return theIsAnnealed; }
0055
0056 void DeterministicAnnealing::debug() const {
0057   cout << "[DeterministicAnnealing] schedule=";
0058   for (vector<float>::const_iterator i = theTemperatures.begin(); i != theTemperatures.end(); ++i) {
0059     cout << *i << " ";
0060   };
0061   cout << endl;
0062 }