LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
CMatchManager.cxx
Go to the documentation of this file.
2 
3 #include <algorithm>
4 #include <iostream>
5 #include <set>
6 #include <stdlib.h>
7 #include <string>
8 #include <utility>
9 #include <vector>
10 
11 #include "TStopwatch.h"
18 
19 namespace cmtool {
20 
22  {
23  throw CMTException("Default ctor needs # planes as an argument!");
24  }
25 
27  {
28  _match_algo = nullptr;
29  _nplanes = nplanes;
30  Reset();
31  }
32 
34  {
39  }
40 
42  {
43  if (_debug_mode <= kPerMerging) {
46  }
47 
50  }
51 
53  {
56  }
57 
59  {
62  }
63 
65  {
68  }
69 
70  unsigned int CMFactorial(unsigned int n)
71  {
72  return (n == 1 || n == 0) ? 1 : CMFactorial(n - 1) * n;
73  }
74 
75  std::vector<std::vector<size_t>> SimpleCombination(size_t n, size_t r)
76  {
77 
78  if (!n || !r) exit(1);
79  if (r > n) exit(1);
80 
81  std::vector<bool> v(n, false);
82  std::fill(v.begin() + n - r, v.end(), true);
83  std::vector<std::vector<size_t>> res;
84  res.reserve(CMFactorial(n) / CMFactorial(n - r) / CMFactorial(r));
85 
86  do {
87  std::vector<size_t> tmp;
88  tmp.reserve(r);
89 
90  for (size_t i = 0; i < n; ++i) {
91  if (v[i]) tmp.push_back(i);
92  }
93  res.push_back(tmp);
94  } while (std::next_permutation(v.begin(), v.end()));
95 
96  return res;
97  }
98 
99  std::vector<std::vector<size_t>> ClusterCombinations(const std::vector<size_t>& seed)
100  {
101 
102  std::vector<size_t> ctr(seed.size(), 0);
103 
104  std::vector<std::vector<size_t>> res;
105 
106  while (1) {
107 
108  res.push_back(std::vector<size_t>(seed.size(), 0));
109  for (size_t index = 0; index < ctr.size(); ++index)
110 
111  (*res.rbegin())[index] = ctr.at(index);
112 
113  for (size_t i = 0; i < ctr.size(); ++i) {
114 
115  size_t index = (size_t)(ctr.size() - i - 1);
116 
117  ctr.at(index) += 1;
118 
119  if (ctr.at(index) < seed.at(index)) break;
120 
121  ctr.at(index) = 0;
122  }
123 
124  bool abort = true;
125  for (auto const& value : ctr)
126 
127  abort = abort && (!(value));
128 
129  if (abort) break;
130  }
131  return res;
132  }
133 
134  std::vector<std::vector<std::pair<size_t, size_t>>> PlaneClusterCombinations(
135  const std::vector<size_t>& seed)
136  {
137  // Result container
138  std::vector<std::vector<std::pair<size_t, size_t>>> result;
139 
140  // Loop over N-planes: start from max number of planes => down to 2 planes
141  for (size_t i = 0; i < seed.size(); ++i) {
142 
143  // If finish making clusters down to 2 palnes, break
144  if (seed.size() < 2 + i) break;
145 
146  // Compute possible N-plane combinations
147  auto const& plane_comb_v = SimpleCombination(seed.size(), seed.size() - i);
148 
149  // Loop over possible N-plane combinations
150  for (auto const& plane_comb : plane_comb_v) {
151 
152  // Make a seed for cluster combinations
153  std::vector<size_t> cluster_seed_v;
154  cluster_seed_v.reserve(plane_comb.size());
155  for (auto const& index : plane_comb)
156  cluster_seed_v.push_back(seed[index]);
157 
158  // Compute cluster combinations
159  for (auto const& cluster_comb : ClusterCombinations(cluster_seed_v)) {
160 
161  // Store result
162  result.push_back(std::vector<std::pair<size_t, size_t>>());
163  for (size_t i = 0; i < cluster_comb.size(); ++i)
164 
165  (*result.rbegin()).push_back(std::make_pair(plane_comb.at(i), cluster_comb.at(i)));
166  }
167  }
168  }
169  return result;
170  }
171 
173  {
174  TStopwatch localWatch;
175 
176  //
177  // Create plane-by-plane vectors
178  //
180 
181  if (_planes.size() < 2) return false;
182 
183  if (_planes.size() > _nplanes)
184  throw CMTException("Found more plane IDs than specified number of planes!");
185 
186  // Index array of clusters per plane
187  std::vector<std::vector<size_t>> cluster_array;
188 
189  // Resize to # of planes w/ clusters
190  cluster_array.reserve(_planes.size());
191 
192  // plane-to-index map
193  std::vector<size_t> plane_to_index(_nplanes, _nplanes);
194 
195  // Fill plane-to-index map
196  for (size_t plane = 0; plane < _nplanes; ++plane) {
197  if (_planes.find(plane) != _planes.end()) {
198  plane_to_index[plane] = cluster_array.size();
199 
200  cluster_array.push_back(std::vector<size_t>());
201  }
202  }
203 
204  // Fill cluster_array
205  for (auto riter = _priority.rbegin(); riter != _priority.rend(); ++riter)
206  cluster_array.at(plane_to_index.at(_in_clusters.at((*riter).second).Plane()))
207  .push_back((*riter).second);
208 
209  // Find combinations
210  std::vector<size_t> seed;
211  seed.reserve(cluster_array.size());
212  for (auto const& clusters_per_plane : cluster_array)
213  seed.push_back(clusters_per_plane.size());
214 
215  auto const& combinations = PlaneClusterCombinations(seed);
216 
217  // Loop over combinations and call algorithm
218  for (auto const& comb : combinations) {
219 
220  std::vector<const cluster::ClusterParamsAlg*> ptr_v;
221 
222  std::vector<unsigned int> tmp_index_v;
223  tmp_index_v.reserve(comb.size());
224 
225  ptr_v.reserve(comb.size());
226 
227  for (auto const& plane_cluster : comb) {
228  auto const& in_cluster_index =
229  cluster_array.at(plane_cluster.first).at(plane_cluster.second);
230 
231  tmp_index_v.push_back(in_cluster_index);
232 
233  ptr_v.push_back(&(_in_clusters.at(in_cluster_index)));
234  }
235 
236  if (_debug_mode <= kPerMerging) {
237  std::cout << " \033[93m"
238  << "Inspecting a pair (";
239  for (auto const& index : tmp_index_v)
240  std::cout << index << " ";
241  std::cout << ") \033[00m" << std::flush;
242 
243  localWatch.Start();
244  }
245 
246  auto const& score = _match_algo->Float(gser, ptr_v);
247 
248  if (_debug_mode <= kPerMerging)
249  std::cout << " ... Time taken = " << localWatch.RealTime() << " [s]" << std::endl;
250 
251  if (score > 0) _book_keeper.Match(tmp_index_v, score);
252  }
253 
254  if (_debug_mode <= kPerIteration) {
257  }
258 
259  return false;
260  }
261 
262 }
unsigned int CMFactorial(unsigned int n)
TRandom r
Definition: spectrum.C:23
virtual void EventEnd()
FMWK function called @ end of Process()
Somewhat verbose (cout per merging iteration)
Definition: CMManagerBase.h:48
Class def header for a class CMatchBookKeeper.
virtual void Reset()
Method to reset itself.
virtual void Report()
Definition: CMAlgoBase.h:69
Class def header for a class CFloatAlgoBase.
std::vector< std::vector< std::pair< size_t, size_t > > > PlaneClusterCombinations(const std::vector< size_t > &seed)
::cmtool::CFloatAlgoBase * _match_algo
Merging algorithm.
Definition: CMatchManager.h:76
Float_t tmp
Definition: plot.C:35
CMatchManager()
Default constructor is private because we need an argument to configure w/ # planes in the detector...
virtual float Float(util::GeometryUtilities const &, const std::vector< const cluster::ClusterParamsAlg * > &clusters)
std::vector< std::vector< size_t > > ClusterCombinations(const std::vector< size_t > &seed)
CMatchBookKeeper _book_keeper
Book keeper instance.
Definition: CMatchManager.h:73
void ComputePriority(const std::vector< cluster::ClusterParamsAlg > &clusters)
Function to compute priority.
virtual void EventEnd()
Definition: CMAlgoBase.h:50
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
virtual bool IterationProcess(util::GeometryUtilities const &gser)
FMWK function called @ iterative loop inside Process()
Class def header for a class CPriorityAlgoBase.
virtual void EventBegin()
FMWK function called @ beginning of Process()
void Match(const std::vector< unsigned int > &matched_indexes, const float &score)
Method to register matched clusters.
std::multimap< float, size_t > _priority
Priority record.
virtual void Reset()
Function to reset the algorithm instance called within CMergeManager/CMatchManager&#39;s Reset() ...
Definition: CMAlgoBase.h:40
long seed
Definition: chem4.cc:67
Class def header for exception classes in CMTException.
void fill(const art::PtrVector< recob::Hit > &hits, int only_plane)
CMMSGLevel_t _debug_mode
Debug mode switch.
double value
Definition: spectrum.C:18
Class def header for a class CMatchManager.
std::vector< cluster::ClusterParamsAlg > _in_clusters
Input clusters.
virtual void IterationEnd()
FMWK function called @ end of iterative loop inside Process()
::cmtool::CPriorityAlgoBase * _priority_algo
Priority algorithm.
std::set< UChar_t > _planes
A holder for # of unique planes in the clusters, computed in ComputePriority() function.
void Reset()
Reset method.
virtual void SetVerbose(bool doit=true)
Setter function for verbosity.
Definition: CMAlgoBase.h:75
Char_t n[5]
virtual void IterationBegin(const std::vector< cluster::ClusterParamsAlg > &)
Definition: CMAlgoBase.h:57
Extremely verbose (cout per individual algorithm execution)
Definition: CMManagerBase.h:46
virtual void IterationBegin()
FMWK function called @ beginning of iterative loop inside Process()
Class def header for a class CMManagerBase.
std::vector< std::vector< size_t > > SimpleCombination(size_t n, size_t r)
size_t _nplanes
Number of planes.
Definition: CMatchManager.h:79
virtual void IterationEnd()
Definition: CMAlgoBase.h:62
void Reset()
Method to reset itself.