GsfTools/interface/CloseComponentsMerger.icc

0001 #include "TrackingTools/GsfTools/interface/MultiGaussianState.h"
0002 #include "TrackingTools/GsfTools/interface/MultiGaussianStateAssembler.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 #include <algorithm>
0006 #include <queue>
0007 #include <cfloat>
0008 #include "CommonTools/Utils/interface/DynArray.h"
0009
0010 template <unsigned int N>
0011 CloseComponentsMerger<N>::CloseComponentsMerger(int maxNumberOfComponents, const DistanceBetweenComponents<N>* distance)
0012     : theMaxNumberOfComponents(maxNumberOfComponents), theDistance(distance->clone()) {}
0013
0014 template <unsigned int N>
0015 MultiGaussianState<N> CloseComponentsMerger<N>::merge(const MultiState& mgs) const {
0016   //auto oldRes = mergeOld(mgs);
0017
0018   typedef std::vector<SingleStatePtr> SingleStateVector;
0019
0020   auto ori = mgs.components();
0021   SingleStateVector finalComponents;
0022
0023   int noComp = ori.size();
0024   if (noComp <= theMaxNumberOfComponents)
0025     return mgs;
0026
0027   while (true) {  // termitates when the nunmber of components becomes less than allowed maximum
0028     SingleStateVector merged;
0029     merged.reserve(noComp);
0030
0031     declareDynArray(float, noComp, weights);
0032     initDynArray(bool, noComp, active, true);
0033     for (int i = 0; i < noComp; ++i) {
0034       weights[i] = ori[i]->weight();
0035     }
0036
0037     auto cmp = [&](int i, int j) { return weights[i] > weights[j]; };
0038     unInitDynArray(int, noComp, qst);  // queue storage
0039     std::priority_queue<int, DynArray<int>, decltype(cmp)> toMerge(cmp, std::move(qst));
0040     for (int i = 0; i < noComp; ++i)
0041       toMerge.push(i);
0042
0043     auto minDistToMax = [&]() -> int {
0044       auto mind = std::numeric_limits<double>::max();
0045       int im = 0;
0046       auto topI = toMerge.top();
0047       auto const& tc = *ori[topI];
0048       active[topI] = false;
0049       for (int i = 0; i < noComp; ++i) {
0050         if (!active[i])
0051           continue;
0052         // assert(weights[topI]<=weights[i]);
0053         auto dist = (*theDistance)(tc, *ori[i]);
0054         if (dist < mind) {
0055           mind = dist;
0056           im = i;
0057         }
0058       }
0059       return im;
0060     };
0061
0062     auto nComp = noComp;
0063     auto nAct = nComp;
0064     while ((nAct > 0) & (nComp > theMaxNumberOfComponents)) {
0065       if (nAct == 1) {
0066         merged.push_back(ori[toMerge.top()]);
0067         nAct = 0;
0068         break;
0069       }
0070
0071       auto ii = minDistToMax();
0072       SingleStateVector comp;
0073       comp.push_back(ori[toMerge.top()]);
0074       comp.push_back(ori[ii]);
0075       active[ii] = false;
0076       while ((!toMerge.empty()) && (!active[toMerge.top()])) {
0077         toMerge.pop();
0078       }
0079       --nComp;
0080       nAct -= 2;
0081       merged.push_back(MultiGaussianStateCombiner<N>().combine(comp));
0082       // assert(toMerge.size()>=0);
0083       // assert(int(toMerge.size())>=nAct);
0084     }
0085
0086     // std::cout << theMaxNumberOfComponents << ' ' << noComp << ' ' << nComp << ' ' << nAct << ' ' << toMerge.size() << std::endl;
0087
0088     if (nComp <= theMaxNumberOfComponents) {  // end game
0089
0090       /*
0091         {
0092           int kk=merged.size();
0093           for (int i=0; i<noComp; ++i) { if (active[i]) ++kk;}
0094           assert(kk==nComp);
0095         }
0096         */
0097
0098       MultiGaussianStateAssembler<N> result;
0099
0100       for (int i = 0; i < noComp; ++i) {
0101         if (active[i])
0102           result.addState(ori[i]);
0103       }
0104       for (auto&& s : merged)
0105         result.addState(s);
0106
0107       auto res = result.combinedState();
0108
0109       /*
0110       assert(res.components().size()==oldRes.components().size());
0111       std::cout << res.weight() << ' ' << oldRes.weight() << std::endl;
0112       for(auto const & c : res.components()) std::cout << c->weight() << '/';
0113       std::cout << std::endl;
0114       for(auto const & c : oldRes.components()) std::cout << c->weight() << '/';
0115       std::cout << std::endl;
0116       */
0117
0118       return res;
0119     }
0120
0121     // assert(nAct==0);
0122     std::swap(ori, merged);
0123     noComp = ori.size();
0124   }
0125
0126   //
0127 }
0128
0129 template <unsigned int N>
0130 MultiGaussianState<N> CloseComponentsMerger<N>::mergeOld(const MultiState& mgs) const {
0131   typedef std::vector<SingleStatePtr> SingleStateVector;
0132   SingleStateVector unmergedComponents(mgs.components());
0133   SingleStateVector finalComponents;
0134   int nComp = unmergedComponents.size();
0135
0136   if (unmergedComponents.empty()) {
0137     edm::LogError("CloseComponentsMerger") << "Trying to merge trajectory state with zero components!";
0138     return mgs;  // ThS: TSOS();
0139   }
0140
0141   // ThS: Don't you want to throw an exception at construction of the class?
0142   if (theMaxNumberOfComponents <= 0) {
0143     edm::LogError("CloseComponentsMerger")
0144         << "Trying to merge state into zero (or less!) components, returning invalid state!";
0145     return mgs;  // ThS: TSOS();
0146   }
0147
0148   // ThS: Of course, here the TSOS will not be invalid. But it will have 0 components
0149   if (mgs.weight() == 0) {
0150     edm::LogError("CloseComponentsMerger") << "Trying to merge mixture with sum of weights equal to zero!";
0151     //     return mgs.createNewState(finalComponents);
0152     return MultiState();
0153   }
0154
0155   if (nComp < theMaxNumberOfComponents + 1)
0156     return mgs;
0157   // ThS: Why not the initial object, as above?
0158   //      return TSOS(new BasicMultiTrajectoryState(unmergedComponents));
0159
0160   SingleStateMap mapUnmergedComp;
0161   SingleStateMap mapMergedComp;
0162
0163   for (typename SingleStateVector::const_iterator it = unmergedComponents.begin(); it != unmergedComponents.end();
0164        it++) {
0165     mapUnmergedComp.insert(std::make_pair((**it).weight(), *it));
0166   }
0167
0168   while (nComp > theMaxNumberOfComponents) {
0169     mapMergedComp.clear();
0170     while (nComp > theMaxNumberOfComponents && !mapUnmergedComp.empty()) {
0171       if (mapUnmergedComp.size() > 1) {
0172         MinDistResult pairMinDist = compWithMinDistToLargestWeight(mapUnmergedComp);
0173         SingleStateVector comp;
0174         comp.push_back(mapUnmergedComp.begin()->second);
0175         comp.push_back(pairMinDist.first);
0176         mapUnmergedComp.erase(pairMinDist.second);
0177         mapUnmergedComp.erase(mapUnmergedComp.begin());
0178         SingleStatePtr mergedComp = MultiGaussianStateCombiner<N>().combine(comp);
0179         mapMergedComp.insert(std::make_pair(mergedComp->weight(), mergedComp));
0180         nComp--;
0181       } else {
0182         mapMergedComp.insert(std::make_pair(mapUnmergedComp.begin()->first, mapUnmergedComp.begin()->second));
0183         mapUnmergedComp.erase(mapUnmergedComp.begin());
0184       }
0185     }
0186     if (mapUnmergedComp.empty() && nComp > theMaxNumberOfComponents) {
0187       mapUnmergedComp = mapMergedComp;
0188     }
0189   }
0190
0191   MultiGaussianStateAssembler<N> result;
0192
0193   for (typename SingleStateMap::const_iterator it = mapUnmergedComp.begin(); it != mapUnmergedComp.end(); it++) {
0194     result.addState(it->second);
0195   }
0196
0197   for (typename SingleStateMap::const_iterator it = mapMergedComp.begin(); it != mapMergedComp.end(); it++) {
0198     result.addState(it->second);
0199   }
0200
0201   return result.combinedState();
0202 }
0203
0204 template <unsigned int N>
0205 typename CloseComponentsMerger<N>::MinDistResult CloseComponentsMerger<N>::compWithMinDistToLargestWeight(
0206     SingleStateMap& unmergedComp) const {
0207   // template <unsigned int N>
0208   // std::pair<SingleGaussianState<N>,
0209   //      typename std::multimap<double, std::shared_ptr< MultiGaussianState<N> > >::iterator>
0210   // CloseComponentsMerger<N>::compWithMinDistToLargestWeight(SingleStateMap& unmergedComp) const {
0211   double large = DBL_MAX;
0212   double minDist = large;
0213   typename SingleStateMap::iterator iterMinDist(0);
0214   for (typename SingleStateMap::iterator it = unmergedComp.begin(); it != unmergedComp.end(); it++) {
0215     if (it != unmergedComp.begin()) {
0216       // assert(unmergedComp.begin()->first >=it->first);
0217       double dist = (*theDistance)(*unmergedComp.begin()->second, *it->second);
0218       if (dist < minDist) {
0219         iterMinDist = it;
0220         minDist = dist;
0221       }
0222     }
0223   }
0224   //   SingleStatePtr minDistComp(iterMinDist->second);
0225   //   return std::make_pair(minDistComp, iterMinDist);
0226   return std::make_pair(iterMinDist->second, iterMinDist);
0227 }