EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RayFrustumBenchmark.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RayFrustumBenchmark.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2017-2018 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 #include <algorithm>
10 #include <chrono>
11 #include <functional>
12 #include <iostream>
13 #include <map>
14 #include <random>
15 #include <vector>
16 
17 #include "Acts/Surfaces/BoundaryCheck.hpp"
18 #include "Acts/Tests/CommonHelpers/BenchmarkTools.hpp"
19 #include "Acts/Utilities/BoundingBox.hpp"
20 #include "Acts/Utilities/Definitions.hpp"
21 #include "Acts/Utilities/Frustum.hpp"
22 #include "Acts/Utilities/Ray.hpp"
23 #include "Acts/Utilities/Units.hpp"
24 
25 using namespace Acts;
26 
27 struct O {};
28 using Box = Acts::AxisAlignedBoundingBox<O, float, 3>;
29 using VertexType = Box::VertexType;
30 using vertex_array_type = Box::vertex_array_type;
31 using value_type = Box::value_type;
32 using Frustum3 = Frustum<float, 3, 4>;
33 using Ray3 = Ray<float, 3>;
34 
35 int main(int /*argc*/, char** /*argv[]*/) {
36  size_t n = 1000;
37 
38  std::mt19937 rng(42);
39  std::uniform_real_distribution<float> dir(0, 1);
40  std::uniform_real_distribution<float> loc(-10, 10);
41  std::uniform_real_distribution<float> ang(M_PI / 10., M_PI / 4.);
42 
43  Box testBox{nullptr, {0, 0, 0}, Box::Size{{1, 2, 3}}};
44 
45  std::cout << "\n==== RAY ====\n" << std::endl;
46 
47  std::map<std::string, bool (*)(const Box&, const Ray3&)> rayVariants;
48 
49  rayVariants["Nominal"] = [](const auto& box, const auto& ray) -> bool {
50  return box.intersect(ray);
51  };
52 
53  rayVariants["Incl. div., unroll"] = [](const Box& box,
54  const Ray<float, 3>& ray) {
55  const VertexType& origin = ray.origin();
56 
57  const vertex_array_type& d = ray.dir();
58 
59  double tmin = -INFINITY, tmax = INFINITY;
60  if (d.x() != 0.0) {
61  double tx1 = (box.min().x() - origin.x()) / d.x();
62  double tx2 = (box.max().x() - origin.x()) / d.x();
63 
64  tmin = std::max(tmin, std::min(tx1, tx2));
65  tmax = std::min(tmax, std::max(tx1, tx2));
66  }
67 
68  if (d.y() != 0.0) {
69  double ty1 = (box.min().y() - origin.y()) / d.y();
70  double ty2 = (box.max().y() - origin.y()) / d.y();
71 
72  tmin = std::max(tmin, std::min(ty1, ty2));
73  tmax = std::min(tmax, std::max(ty1, ty2));
74  }
75 
76  if (d.z() != 0.0) {
77  double tz1 = (box.min().z() - origin.z()) / d.z();
78  double tz2 = (box.max().z() - origin.z()) / d.z();
79 
80  tmin = std::max(tmin, std::min(tz1, tz2));
81  tmax = std::min(tmax, std::max(tz1, tz2));
82  }
83 
84  return tmax > tmin && tmax > 0.0;
85  };
86 
87  rayVariants["Incl. div., loop"] = [](const Box& box,
88  const Ray<float, 3>& ray) {
89  const VertexType& origin = ray.origin();
90 
91  const vertex_array_type& d = ray.dir();
92 
93  double tmin = -INFINITY, tmax = INFINITY;
94 
95  for (size_t i = 0; i < 3; i++) {
96  if (d[i] == 0.0)
97  continue;
98  double t1 = (box.min()[i] - origin[i]) / d[i];
99  double t2 = (box.max()[i] - origin[i]) / d[i];
100  tmin = std::max(tmin, std::min(t1, t2));
101  tmax = std::min(tmax, std::max(t1, t2));
102  }
103 
104  return tmax > tmin && tmax > 0.0;
105  };
106 
107  rayVariants["Incl. div., min/max alt., unroll"] =
108  [](const Box& box, const Ray<float, 3>& ray) {
109  const VertexType& origin = ray.origin();
110  const vertex_array_type& d = ray.dir();
111 
112  double tx1 = (box.min().x() - origin.x()) / d.x();
113  double tx2 = (box.max().x() - origin.x()) / d.x();
114  double tmin = std::min(tx1, tx2);
115  double tmax = std::max(tx1, tx2);
116 
117  double ty1 = (box.min().y() - origin.y()) / d.y();
118  double ty2 = (box.max().y() - origin.y()) / d.y();
119  tmin = std::max(tmin, std::min(ty1, ty2));
120  tmax = std::min(tmax, std::max(ty1, ty2));
121 
122  double tz1 = (box.min().z() - origin.z()) / d.z();
123  double tz2 = (box.max().z() - origin.z()) / d.z();
124  tmin = std::max(tmin, std::min(tz1, tz2));
125  tmax = std::min(tmax, std::max(tz1, tz2));
126 
127  return tmax > tmin && tmax > 0.0;
128  };
129 
130  rayVariants["No div., min/max alt, unroll"] = [](const Box& box,
131  const Ray<float, 3>& ray) {
132  const VertexType& origin = ray.origin();
133  const vertex_array_type& id = ray.idir();
134 
135  double tx1 = (box.min().x() - origin.x()) * id.x();
136  double tx2 = (box.max().x() - origin.x()) * id.x();
137  double tmin = std::min(tx1, tx2);
138  double tmax = std::max(tx1, tx2);
139 
140  double ty1 = (box.min().y() - origin.y()) * id.y();
141  double ty2 = (box.max().y() - origin.y()) * id.y();
142  tmin = std::max(tmin, std::min(ty1, ty2));
143  tmax = std::min(tmax, std::max(ty1, ty2));
144 
145  double tz1 = (box.min().z() - origin.z()) * id.z();
146  double tz2 = (box.max().z() - origin.z()) * id.z();
147  tmin = std::max(tmin, std::min(tz1, tz2));
148  tmax = std::min(tmax, std::max(tz1, tz2));
149 
150  return tmax > tmin && tmax > 0.0;
151  };
152 
153  rayVariants["No div., min/max orig, loop"] = [](const Box& box,
154  const Ray<float, 3>& ray) {
155  const VertexType& origin = ray.origin();
156  const vertex_array_type& id = ray.idir();
157  double tmin = -INFINITY, tmax = INFINITY;
158 
159  for (size_t i = 0; i < 3; i++) {
160  double t1 = (box.min()[i] - origin[i]) * id[i];
161  double t2 = (box.max()[i] - origin[i]) * id[i];
162  tmin = std::max(tmin, std::min(t1, t2));
163  tmax = std::min(tmax, std::max(t1, t2));
164  }
165 
166  return tmax > tmin && tmax > 0.0;
167  };
168 
169  std::vector<Ray3> rays{n, Ray3{Vector3F{0, 0, 0}, Vector3F{1, 0, 0}}};
170  std::generate(rays.begin(), rays.end(), [&]() {
171  const Vector3F d{dir(rng), dir(rng), dir(rng)};
172  const Vector3F l{loc(rng), loc(rng), loc(rng)};
173  return Ray3{l, d.normalized()};
174  });
175 
176  std::cout << "Make sure ray implementations are identical" << std::endl;
177  for (const auto& ray : rays) {
178  std::vector<std::pair<std::string, bool>> results;
179 
180  std::transform(rayVariants.begin(), rayVariants.end(),
181  std::back_inserter(results),
182  [&](const auto& p) -> decltype(results)::value_type {
183  const auto& [name, func] = p;
184  return {name, func(testBox, ray)};
185  });
186 
187  bool all = std::all_of(results.begin(), results.end(),
188  [](const auto& r) { return r.second; });
189  bool none = std::none_of(results.begin(), results.end(),
190  [](const auto& r) { return r.second; });
191 
192  if (!all && !none) {
193  std::cerr << "Discrepancy: " << std::endl;
194  for (const auto& [name, result] : results) {
195  std::cerr << " - " << name << ": " << result << std::endl;
196  }
197 
198  testBox.toStream(std::cerr);
199  std::cerr << std::endl;
200  std::cerr << "Ray: [" << ray.origin().transpose() << "], ["
201  << ray.dir().transpose() << "]" << std::endl;
202  return -1;
203  }
204  }
205  std::cout << "Seems ok" << std::endl;
206 
207  std::cout << "Run benchmarks: " << std::endl;
208  for (const auto& p : rayVariants) {
209  // can't capture structured binding, so pair access it is.
210  std::cout << "- Benchmarking variant: '" << p.first << "'" << std::endl;
211  auto bench_result = Acts::Test::microBenchmark(
212  [&](const auto& ray) { return p.second(testBox, ray); }, rays);
213  std::cout << " " << bench_result << std::endl;
214  }
215 
216  std::cout << "\n==== FRUSTUM ====\n" << std::endl;
217 
218  std::map<std::string, bool (*)(const Box&, const Frustum3&)> frustumVariants;
219 
220  frustumVariants["Nominal"] = [](const auto& box,
221  const auto& frustum) -> bool {
222  return box.intersect(frustum);
223  };
224 
225  frustumVariants["Manual constexpr loop unroll, early ret."] =
226  [](const Box& box, const Frustum3& fr) {
227  constexpr size_t sides = 4; // yes this is pointless, i just want to
228  // kind of match the other impl
229 
230  const auto& normals = fr.normals();
231  const vertex_array_type fr_vmin = box.min() - fr.origin();
232  const vertex_array_type fr_vmax = box.max() - fr.origin();
233 
234  auto calc = [&](const auto& normal) {
235  return (normal.array() < 0).template cast<value_type>() * fr_vmin +
236  (normal.array() >= 0).template cast<value_type>() * fr_vmax;
237  };
238 
239  VertexType p_vtx;
240 
241  p_vtx = calc(normals[0]);
242  if (p_vtx.dot(normals[0]) < 0) {
243  return false;
244  }
245 
246  p_vtx = calc(normals[1]);
247  if (p_vtx.dot(normals[1]) < 0) {
248  return false;
249  }
250 
251  p_vtx = calc(normals[2]);
252  if (p_vtx.dot(normals[2]) < 0) {
253  return false;
254  }
255 
256  if constexpr (sides > 2) {
257  p_vtx = calc(normals[3]);
258  if (p_vtx.dot(normals[3]) < 0) {
259  return false;
260  }
261  }
262 
263  if constexpr (sides > 3) {
264  p_vtx = calc(normals[4]);
265  if (p_vtx.dot(normals[4]) < 0) {
266  return false;
267  }
268  }
269 
270  if constexpr (sides > 4) {
271  for (size_t i = 5; i <= fr.sides; i++) {
272  const VertexType& normal = normals[i];
273 
274  p_vtx = calc(normal);
275  if (p_vtx.dot(normal) < 0) {
276  return false;
277  }
278  }
279  }
280 
281  return true;
282  };
283 
284  frustumVariants["Nominal, no early ret."] = [](const Box& box,
285  const Frustum3& fr) {
286  const auto& normals = fr.normals();
287  const vertex_array_type fr_vmin = box.min() - fr.origin();
288  const vertex_array_type fr_vmax = box.max() - fr.origin();
289 
290  VertexType p_vtx;
291  bool result = true;
292  for (size_t i = 0; i < fr.sides + 1; i++) {
293  const VertexType& normal = normals[i];
294 
295  p_vtx = (normal.array() < 0).template cast<value_type>() * fr_vmin +
296  (normal.array() >= 0).template cast<value_type>() * fr_vmax;
297 
298  result = result && (p_vtx.dot(normal) >= 0);
299  }
300  return result;
301  };
302 
303  frustumVariants["Manual constexpr unroll, early ret."] =
304  [](const Box& box, const Frustum3& fr) {
305  constexpr size_t sides = 4; // yes this is pointless, i just want to
306  // kind of match the other impl
307 
308  const auto& normals = fr.normals();
309  const vertex_array_type fr_vmin = box.min() - fr.origin();
310  const vertex_array_type fr_vmax = box.max() - fr.origin();
311 
312  auto calc = [&](const auto& normal) {
313  return (normal.array() < 0).template cast<value_type>() * fr_vmin +
314  (normal.array() >= 0).template cast<value_type>() * fr_vmax;
315  };
316 
317  VertexType p_vtx;
318  bool result = true;
319 
320  p_vtx = calc(normals[0]);
321  result = result && (p_vtx.dot(normals[0]) >= 0);
322 
323  p_vtx = calc(normals[1]);
324  result = result && (p_vtx.dot(normals[1]) >= 0);
325 
326  p_vtx = calc(normals[2]);
327  result = result && (p_vtx.dot(normals[2]) >= 0);
328 
329  if constexpr (sides > 2) {
330  p_vtx = calc(normals[3]);
331  result = result && (p_vtx.dot(normals[3]) >= 0);
332  }
333 
334  if constexpr (sides > 3) {
335  p_vtx = calc(normals[4]);
336  result = result && (p_vtx.dot(normals[4]) >= 0);
337  }
338 
339  if constexpr (sides > 4) {
340  for (size_t i = 5; i <= fr.sides; i++) {
341  const VertexType& normal = normals[i];
342 
343  p_vtx = calc(normal);
344  result = result && (p_vtx.dot(normal) >= 0);
345  }
346  }
347 
348  return result;
349  };
350 
351  std::vector<Frustum3> frustums{n, Frustum3{{0, 0, 0}, {1, 0, 0}, M_PI / 2.}};
352  std::generate(frustums.begin(), frustums.end(), [&]() {
353  const Vector3F d{dir(rng), dir(rng), dir(rng)};
354  const Vector3F l{loc(rng), loc(rng), loc(rng)};
355  return Frustum3{l, d.normalized(), ang(rng)};
356  });
357 
358  std::cout << "Make sure frustum implementations are identical" << std::endl;
359  for (const auto& fr : frustums) {
360  std::vector<std::pair<std::string, bool>> results;
361 
362  std::transform(frustumVariants.begin(), frustumVariants.end(),
363  std::back_inserter(results),
364  [&](const auto& p) -> decltype(results)::value_type {
365  const auto& [name, func] = p;
366  return {name, func(testBox, fr)};
367  });
368 
369  bool all = std::all_of(results.begin(), results.end(),
370  [](const auto& r) { return r.second; });
371  bool none = std::none_of(results.begin(), results.end(),
372  [](const auto& r) { return r.second; });
373 
374  if (!all && !none) {
375  std::cerr << "Discrepancy: " << std::endl;
376  for (const auto& [name, result] : results) {
377  std::cerr << " - " << name << ": " << result << std::endl;
378  }
379 
380  testBox.toStream(std::cerr);
381  std::cerr << std::endl;
382  std::cerr << "Frustum: [" << fr.origin().transpose() << "], ["
383  << fr.dir().transpose() << "]" << std::endl;
384  return -1;
385  }
386  }
387  std::cout << "Seems ok" << std::endl;
388 
389  size_t iters_per_run = 1000;
390 
391  std::vector<std::pair<std::string, Frustum3>> testFrusts = {
392  {"away", Frustum3{{0, 0, -10}, {0, 0, -1}, M_PI / 4.}},
393  {"towards", Frustum3{{0, 0, -10}, {0, 0, 1}, M_PI / 4.}},
394  {"left", Frustum3{{0, 0, -10}, {0, 1, 0}, M_PI / 4.}},
395  {"right", Frustum3{{0, 0, -10}, {0, -1, 0}, M_PI / 4.}},
396  {"up", Frustum3{{0, 0, -10}, {1, 0, 0}, M_PI / 4.}},
397  {"down", Frustum3{{0, 0, -10}, {-1, 0, 0}, M_PI / 4.}},
398  };
399 
400  std::cout << "Run benchmarks: " << std::endl;
401 
402  for (const auto& fr_pair : testFrusts) {
403  std::cout << "Frustum '" << fr_pair.first << "'" << std::endl;
404 
405  for (const auto& p : frustumVariants) {
406  // can't capture structured binding, so pair access it is.
407  std::cout << "- Benchmarking variant: '" << p.first << "'" << std::endl;
408  auto bench_result = Acts::Test::microBenchmark(
409  [&]() { return p.second(testBox, fr_pair.second); }, iters_per_run);
410  std::cout << " " << bench_result << std::endl;
411  }
412 
413  std::cout << std::endl;
414  }
415  return 0;
416 }