TrackerTFP/src/MiniHoughTransform.cc

0001 #include "L1Trigger/TrackerTFP/interface/MiniHoughTransform.h"
0002
0003 #include <numeric>
0004 #include <algorithm>
0005 #include <iterator>
0006 #include <deque>
0007 #include <vector>
0008 #include <set>
0009 #include <utility>
0010 #include <cmath>
0011
0012 using namespace std;
0013 using namespace edm;
0014 using namespace tt;
0015
0016 namespace trackerTFP {
0017
0018   MiniHoughTransform::MiniHoughTransform(const ParameterSet& iConfig,
0019                                          const Setup* setup,
0020                                          const DataFormats* dataFormats,
0021                                          int region)
0022       : enableTruncation_(iConfig.getParameter<bool>("EnableTruncation")),
0023         setup_(setup),
0024         dataFormats_(dataFormats),
0025         inv2R_(dataFormats_->format(Variable::inv2R, Process::ht)),
0026         phiT_(dataFormats_->format(Variable::phiT, Process::ht)),
0027         region_(region),
0028         numBinsInv2R_(setup_->htNumBinsInv2R()),
0029         numCells_(setup_->mhtNumCells()),
0030         numNodes_(setup_->mhtNumDLBNodes()),
0031         numChannel_(setup_->mhtNumDLBChannel()),
0032         input_(numBinsInv2R_) {}
0033
0034   // read in and organize input product (fill vector input_)
0035   void MiniHoughTransform::consume(const StreamsStub& streams) {
0036     auto valid = [](int sum, const FrameStub& frame) { return sum + (frame.first.isNonnull() ? 1 : 0); };
0037     int nStubsHT(0);
0038     for (int binInv2R = 0; binInv2R < numBinsInv2R_; binInv2R++) {
0039       const StreamStub& stream = streams[region_ * numBinsInv2R_ + binInv2R];
0040       nStubsHT += accumulate(stream.begin(), stream.end(), 0, valid);
0041     }
0042     stubsHT_.reserve(nStubsHT);
0043     stubsMHT_.reserve(nStubsHT * numCells_);
0044     for (int binInv2R = 0; binInv2R < numBinsInv2R_; binInv2R++) {
0045       const int inv2R = inv2R_.toSigned(binInv2R);
0046       const StreamStub& stream = streams[region_ * numBinsInv2R_ + binInv2R];
0047       vector<StubHT*>& stubs = input_[binInv2R];
0048       stubs.reserve(stream.size());
0049       // Store input stubs in vector, so rest of MHT algo can work with pointers to them (saves CPU)
0050       for (const FrameStub& frame : stream) {
0051         StubHT* stub = nullptr;
0052         if (frame.first.isNonnull()) {
0053           stubsHT_.emplace_back(frame, dataFormats_, inv2R);
0054           stub = &stubsHT_.back();
0055         }
0056         stubs.push_back(stub);
0057       }
0058     }
0059   }
0060
0061   // fill output products
0062   void MiniHoughTransform::produce(StreamsStub& accepted, StreamsStub& lost) {
0063     // fill MHT cells
0064     vector<deque<StubMHT*>> stubsCells(numBinsInv2R_ * numCells_);
0065     for (int channel = 0; channel < numBinsInv2R_; channel++)
0066       fill(channel, input_[channel], stubsCells);
0067     // perform static load balancing
0068     vector<vector<StubMHT*>> streamsSLB(numBinsInv2R_);
0069     for (int channel = 0; channel < numBinsInv2R_; channel++) {
0070       vector<deque<StubMHT*>> tmp(numCells_);
0071       // gather streams to mux together: same MHT cell of 4 adjacent MHT input streams
0072       for (int k = 0; k < numCells_; k++)
0073         swap(tmp[k], stubsCells[(channel / numCells_) * numBinsInv2R_ + channel % numCells_ + k * numCells_]);
0074       slb(tmp, streamsSLB[channel], lost[channel]);
0075     }
0076     // dynamic load balancing stage 1
0077     vector<vector<StubMHT*>> streamsDLB(numBinsInv2R_);
0078     for (int node = 0; node < numNodes_; node++) {
0079       vector<vector<StubMHT*>> tmp(numChannel_);
0080       // gather streams to dynamically balance them
0081       for (int k = 0; k < numChannel_; k++)
0082         swap(tmp[k], streamsSLB[(node / numCells_) * numNodes_ + node % numCells_ + k * numCells_]);
0083       dlb(tmp);
0084       for (int k = 0; k < numChannel_; k++)
0085         swap(tmp[k], streamsDLB[node * numChannel_ + k]);
0086     }
0087     // dynamic load balancing stage 2
0088     vector<vector<StubMHT*>> streamsMHT(numBinsInv2R_);
0089     for (int node = 0; node < numNodes_; node++) {
0090       vector<vector<StubMHT*>> tmp(numChannel_);
0091       // gather streams to dynamically balance them
0092       for (int k = 0; k < numChannel_; k++)
0093         swap(tmp[k], streamsDLB[node + k * numNodes_]);
0094       dlb(tmp);
0095       for (int k = 0; k < numChannel_; k++)
0096         swap(tmp[k], streamsMHT[node * numChannel_ + k]);
0097     }
0098     // fill output product
0099     for (int channel = 0; channel < numBinsInv2R_; channel++) {
0100       const vector<StubMHT*>& stubs = streamsMHT[channel];
0101       StreamStub& stream = accepted[region_ * numBinsInv2R_ + channel];
0102       stream.reserve(stubs.size());
0103       for (StubMHT* stub : stubs)
0104         stream.emplace_back(stub ? stub->frame() : FrameStub());
0105     }
0106   }
0107
0108   // perform finer pattern recognition per track
0109   void MiniHoughTransform::fill(int channel, const vector<StubHT*>& stubs, vector<deque<StubMHT*>>& streams) {
0110     if (stubs.empty())
0111       return;
0112     int id;
0113     auto differentHT = [&id](StubHT* stub) { return id != stub->trackId(); };
0114     auto differentMHT = [&id](StubMHT* stub) { return !stub || id != stub->trackId(); };
0115     for (auto it = stubs.begin(); it != stubs.end();) {
0116       const auto start = it;
0117       id = (*it)->trackId();
0118       it = find_if(it, stubs.end(), differentHT);
0119       const int size = distance(start, it);
0120       // create finer track candidates stub container
0121       vector<vector<StubMHT*>> mhtCells(numCells_);
0122       for (vector<StubMHT*>& mhtCell : mhtCells)
0123         mhtCell.reserve(size);
0124       // fill finer track candidates stub container
0125       for (auto stub = start; stub != it; stub++) {
0126         const double r = (*stub)->r();
0127         const double chi = (*stub)->phi();
0128         // identify finer track candidates for this stub
0129         // 0 and 1 belong to the MHT cells with larger inv2R; 0 and 2 belong to those with smaller track PhiT
0130         vector<int> cells;
0131         cells.reserve(numCells_);
0132         const bool compA = 2. * abs(chi) < phiT_.base();
0133         const bool compB = 2. * abs(chi) < abs(r * inv2R_.base());
0134         const bool compAB = compA && compB;
0135         if (chi >= 0. && r >= 0.) {
0136           cells.push_back(3);
0137           if (compA)
0138             cells.push_back(1);
0139           if (compAB)
0140             cells.push_back(2);
0141         }
0142         if (chi >= 0. && r < 0.) {
0143           cells.push_back(1);
0144           if (compA)
0145             cells.push_back(3);
0146           if (compAB)
0147             cells.push_back(0);
0148         }
0149         if (chi < 0. && r >= 0.) {
0150           cells.push_back(0);
0151           if (compA)
0152             cells.push_back(2);
0153           if (compAB)
0154             cells.push_back(1);
0155         }
0156         if (chi < 0. && r < 0.) {
0157           cells.push_back(2);
0158           if (compA)
0159             cells.push_back(0);
0160           if (compAB)
0161             cells.push_back(3);
0162         }
0163         // organise stubs in finer track candidates
0164         for (int cell : cells) {
0165           const int inv2R = cell / setup_->mhtNumBinsPhiT();
0166           const int phiT = cell % setup_->mhtNumBinsPhiT();
0167           stubsMHT_.emplace_back(**stub, phiT, inv2R);
0168           mhtCells[cell].push_back(&stubsMHT_.back());
0169         }
0170       }
0171       // perform pattern recognition
0172       for (int sel = 0; sel < numCells_; sel++) {
0173         deque<StubMHT*>& stream = streams[channel * numCells_ + sel];
0174         vector<StubMHT*>& mhtCell = mhtCells[sel];
0175         set<int> layers;
0176         auto toLayer = [](StubMHT* stub) { return stub->layer(); };
0177         transform(mhtCell.begin(), mhtCell.end(), inserter(layers, layers.begin()), toLayer);
0178         if ((int)layers.size() < setup_->mhtMinLayers())
0179           mhtCell.clear();
0180         for (StubMHT* stub : mhtCell)
0181           stream.push_back(stub);
0182         stream.insert(stream.end(), size - (int)mhtCell.size(), nullptr);
0183       }
0184     }
0185     for (int sel = 0; sel < numCells_; sel++) {
0186       deque<StubMHT*>& stream = streams[channel * numCells_ + sel];
0187       // remove all gaps between end and last stub
0188       for (auto it = stream.end(); it != stream.begin();)
0189         it = (*--it) ? stream.begin() : stream.erase(it);
0190       // read out fine track cannot start before rough track has read in completely, add gaps to take this into account
0191       int pos(0);
0192       for (auto it = stream.begin(); it != stream.end();) {
0193         if (!(*it)) {
0194           it = stream.erase(it);
0195           continue;
0196         }
0197         id = (*it)->trackId();
0198         const int s = distance(it, find_if(it, stream.end(), differentMHT));
0199         const int d = distance(stream.begin(), it);
0200         pos += s;
0201         if (d < pos) {
0202           const int diff = pos - d;
0203           it = stream.insert(it, diff, nullptr);
0204           it = next(it, diff);
0205         } else
0206           it = stream.erase(remove(next(stream.begin(), pos), it, nullptr), it);
0207         it = next(it, s);
0208       }
0209       // adjust stream start so that first output stub is in first place in case of quickest track
0210       if (!stream.empty())
0211         stream.erase(stream.begin(), next(stream.begin(), setup_->mhtMinLayers()));
0212     }
0213   }
0214
0215   // Static load balancing of inputs: mux 4 streams to 1 stream
0216   void MiniHoughTransform::slb(vector<deque<StubMHT*>>& inputs, vector<StubMHT*>& accepted, StreamStub& lost) const {
0217     if (all_of(inputs.begin(), inputs.end(), [](const deque<StubMHT*>& stubs) { return stubs.empty(); }))
0218       return;
0219     auto size = [](int sum, const deque<StubMHT*>& stubs) { return sum = stubs.size(); };
0220     const int nFrames = accumulate(inputs.begin(), inputs.end(), 0, size);
0221     accepted.reserve(nFrames);
0222     // input fifos
0223     vector<deque<StubMHT*>> stacks(numCells_);
0224     // helper for handshake
0225     TTBV empty(-1, numCells_, true);
0226     TTBV enable(0, numCells_);
0227     // clock accurate firmware emulation, each while trip describes one clock tick, one stub in and one stub out per tick
0228     while (!all_of(inputs.begin(), inputs.end(), [](const deque<StubMHT*>& d) { return d.empty(); }) or
0229            !all_of(stacks.begin(), stacks.end(), [](const deque<StubMHT*>& d) { return d.empty(); })) {
0230       // store stub in fifo
0231       for (int channel = 0; channel < numCells_; channel++) {
0232         StubMHT* stub = pop_front(inputs[channel]);
0233         if (stub)
0234           stacks[channel].push_back(stub);
0235       }
0236       // identify empty fifos
0237       for (int channel = 0; channel < numCells_; channel++)
0238         empty[channel] = stacks[channel].empty();
0239       // chose new fifo to read from if current fifo got empty
0240       const int iEnableOld = enable.plEncode();
0241       if (enable.none() || empty[iEnableOld]) {
0242         enable.reset();
0243         const int iNotEmpty = empty.plEncode(false);
0244         if (iNotEmpty < numCells_)
0245           enable.set(iNotEmpty);
0246       }
0247       // read from chosen fifo
0248       const int iEnable = enable.plEncode();
0249       if (enable.any())
0250         accepted.push_back(pop_front(stacks[iEnable]));
0251       else
0252         // gap if no fifo has been chosen
0253         accepted.push_back(nullptr);
0254     }
0255     // perform truncation if desired
0256     if (enableTruncation_ && (int)accepted.size() > setup_->numFrames()) {
0257       const auto limit = next(accepted.begin(), setup_->numFrames());
0258       auto valid = [](int sum, StubMHT* stub) { return sum + (stub ? 1 : 0); };
0259       const int nLost = accumulate(limit, accepted.end(), 0, valid);
0260       lost.reserve(nLost);
0261       for (auto it = limit; it != accepted.end(); it++)
0262         if (*it)
0263           lost.emplace_back((*it)->frame());
0264       accepted.erase(limit, accepted.end());
0265     }
0266     // cosmetics -- remove gaps at the end of stream
0267     for (auto it = accepted.end(); it != accepted.begin();)
0268       it = (*--it) == nullptr ? accepted.erase(it) : accepted.begin();
0269   }
0270
0271   // Dynamic load balancing of inputs: swapping parts of streams to balance the amount of tracks per stream
0272   void MiniHoughTransform::dlb(vector<vector<StubMHT*>>& streams) const {
0273     if (all_of(streams.begin(), streams.end(), [](const vector<StubMHT*>& stubs) { return stubs.empty(); }))
0274       return;
0275     auto maxSize = [](int size, const vector<StubMHT*>& stream) { return size = max(size, (int)stream.size()); };
0276     const int nMax = accumulate(streams.begin(), streams.end(), 0, maxSize);
0277     for (vector<StubMHT*>& stream : streams)
0278       stream.resize(nMax, nullptr);
0279     vector<int> prevTrks(numChannel_, -1);
0280     bool swapping(false);
0281     vector<int> loads(numChannel_, 0);
0282     for (int i = 0; i < nMax; i++) {
0283       TTBV newTrks(0, numChannel_);
0284       for (int k = 0; k < numChannel_; k++)
0285         if (!streams[numChannel_ - k - 1][i] && streams[k][i] && streams[k][i]->trackId() != prevTrks[k])
0286           newTrks.set(k);
0287       for (int k = 0; k < numChannel_; k++)
0288         if (newTrks[k])
0289           if ((swapping && loads[numChannel_ - k - 1] > loads[k]) ||
0290               (!swapping && loads[k] > loads[numChannel_ - k - 1]))
0291             swapping = !swapping;
0292       for (int k = 0; k < numChannel_; k++) {
0293         if (streams[k][i])
0294           loads[swapping ? numChannel_ - k - 1 : k]++;
0295         prevTrks[k] = streams[k][i] ? streams[k][i]->trackId() : -1;
0296       }
0297       if (swapping)
0298         swap(streams[0][i], streams[1][i]);
0299     }
0300     // remove all gaps between end and last stub
0301     for (vector<StubMHT*>& stream : streams)
0302       for (auto it = stream.end(); it != stream.begin();)
0303         it = (*--it) ? stream.begin() : stream.erase(it);
0304   }
0305
0306   // remove and return first element of deque, returns nullptr if empty
0307   template <class T>
0308   T* MiniHoughTransform::pop_front(deque<T*>& ts) const {
0309     T* t = nullptr;
0310     if (!ts.empty()) {
0311       t = ts.front();
0312       ts.pop_front();
0313     }
0314     return t;
0315   }
0316
0317 }  // namespace trackerTFP