Analysis Software
Documentation for sPHENIX simulation 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 "Acts/Definitions/Algebra.hpp"
10 #include "Acts/Definitions/Units.hpp"
11 #include "Acts/Surfaces/BoundaryCheck.hpp"
12 #include "Acts/Tests/CommonHelpers/BenchmarkTools.hpp"
13 #include "Acts/Utilities/BoundingBox.hpp"
14 #include "Acts/Utilities/Frustum.hpp"
15 #include "Acts/Utilities/Ray.hpp"
16 
17 #include <algorithm>
18 #include <chrono>
19 #include <functional>
20 #include <iostream>
21 #include <map>
22 #include <random>
23 #include <vector>
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  }
99  double t1 = (box.min()[i] - origin[i]) / d[i];
100  double t2 = (box.max()[i] - origin[i]) / d[i];
101  tmin = std::max(tmin, std::min(t1, t2));
102  tmax = std::min(tmax, std::max(t1, t2));
103  }
104 
105  return tmax > tmin && tmax > 0.0;
106  };
107 
108  rayVariants["Incl. div., min/max alt., unroll"] =
109  [](const Box& box, const Ray<float, 3>& ray) {
110  const VertexType& origin = ray.origin();
111  const vertex_array_type& d = ray.dir();
112 
113  double tx1 = (box.min().x() - origin.x()) / d.x();
114  double tx2 = (box.max().x() - origin.x()) / d.x();
115  double tmin = std::min(tx1, tx2);
116  double tmax = std::max(tx1, tx2);
117 
118  double ty1 = (box.min().y() - origin.y()) / d.y();
119  double ty2 = (box.max().y() - origin.y()) / d.y();
120  tmin = std::max(tmin, std::min(ty1, ty2));
121  tmax = std::min(tmax, std::max(ty1, ty2));
122 
123  double tz1 = (box.min().z() - origin.z()) / d.z();
124  double tz2 = (box.max().z() - origin.z()) / d.z();
125  tmin = std::max(tmin, std::min(tz1, tz2));
126  tmax = std::min(tmax, std::max(tz1, tz2));
127 
128  return tmax > tmin && tmax > 0.0;
129  };
130 
131  rayVariants["No div., min/max alt, unroll"] = [](const Box& box,
132  const Ray<float, 3>& ray) {
133  const VertexType& origin = ray.origin();
134  const vertex_array_type& id = ray.idir();
135 
136  double tx1 = (box.min().x() - origin.x()) * id.x();
137  double tx2 = (box.max().x() - origin.x()) * id.x();
138  double tmin = std::min(tx1, tx2);
139  double tmax = std::max(tx1, tx2);
140 
141  double ty1 = (box.min().y() - origin.y()) * id.y();
142  double ty2 = (box.max().y() - origin.y()) * id.y();
143  tmin = std::max(tmin, std::min(ty1, ty2));
144  tmax = std::min(tmax, std::max(ty1, ty2));
145 
146  double tz1 = (box.min().z() - origin.z()) * id.z();
147  double tz2 = (box.max().z() - origin.z()) * id.z();
148  tmin = std::max(tmin, std::min(tz1, tz2));
149  tmax = std::min(tmax, std::max(tz1, tz2));
150 
151  return tmax > tmin && tmax > 0.0;
152  };
153 
154  rayVariants["No div., min/max orig, loop"] = [](const Box& box,
155  const Ray<float, 3>& ray) {
156  const VertexType& origin = ray.origin();
157  const vertex_array_type& id = ray.idir();
158  double tmin = -INFINITY, tmax = INFINITY;
159 
160  for (size_t i = 0; i < 3; i++) {
161  double t1 = (box.min()[i] - origin[i]) * id[i];
162  double t2 = (box.max()[i] - origin[i]) * id[i];
163  tmin = std::max(tmin, std::min(t1, t2));
164  tmax = std::min(tmax, std::max(t1, t2));
165  }
166 
167  return tmax > tmin && tmax > 0.0;
168  };
169 
170  using Vector3F = Eigen::Matrix<float, 3, 1>;
171 
172  std::vector<Ray3> rays{n, Ray3{Vector3F{0, 0, 0}, Vector3F{1, 0, 0}}};
173  std::generate(rays.begin(), rays.end(), [&]() {
174  const Vector3F d{dir(rng), dir(rng), dir(rng)};
175  const Vector3F l{loc(rng), loc(rng), loc(rng)};
176  return Ray3{l, d.normalized()};
177  });
178 
179  std::cout << "Make sure ray implementations are identical" << std::endl;
180  for (const auto& ray : rays) {
181  std::vector<std::pair<std::string, bool>> results;
182 
183  std::transform(rayVariants.begin(), rayVariants.end(),
184  std::back_inserter(results),
185  [&](const auto& p) -> decltype(results)::value_type {
186  const auto& [name, func] = p;
187  return {name, func(testBox, ray)};
188  });
189 
190  bool all = std::all_of(results.begin(), results.end(),
191  [](const auto& r) { return r.second; });
192  bool none = std::none_of(results.begin(), results.end(),
193  [](const auto& r) { return r.second; });
194 
195  if (!all && !none) {
196  std::cerr << "Discrepancy: " << std::endl;
197  for (const auto& [name, result] : results) {
198  std::cerr << " - " << name << ": " << result << std::endl;
199  }
200 
201  testBox.toStream(std::cerr);
202  std::cerr << std::endl;
203  std::cerr << "Ray: [" << ray.origin().transpose() << "], ["
204  << ray.dir().transpose() << "]" << std::endl;
205  return -1;
206  }
207  }
208  std::cout << "Seems ok" << std::endl;
209 
210  std::cout << "Run benchmarks: " << std::endl;
211  for (const auto& p : rayVariants) {
212  // can't capture structured binding, so pair access it is.
213  std::cout << "- Benchmarking variant: '" << p.first << "'" << std::endl;
214  auto bench_result = Acts::Test::microBenchmark(
215  [&](const auto& ray) { return p.second(testBox, ray); }, rays);
216  std::cout << " " << bench_result << std::endl;
217  }
218 
219  std::cout << "\n==== FRUSTUM ====\n" << std::endl;
220 
221  std::map<std::string, bool (*)(const Box&, const Frustum3&)> frustumVariants;
222 
223  frustumVariants["Nominal"] = [](const auto& box,
224  const auto& frustum) -> bool {
225  return box.intersect(frustum);
226  };
227 
228  frustumVariants["Manual constexpr loop unroll, early ret."] =
229  [](const Box& box, const Frustum3& fr) {
230  constexpr size_t sides = 4; // yes this is pointless, I just want to
231  // kind of match the other impl
232 
233  const auto& normals = fr.normals();
234  const vertex_array_type fr_vmin = box.min() - fr.origin();
235  const vertex_array_type fr_vmax = box.max() - fr.origin();
236 
237  auto calc = [&](const auto& normal) {
238  return (normal.array() < 0).template cast<value_type>() * fr_vmin +
239  (normal.array() >= 0).template cast<value_type>() * fr_vmax;
240  };
241 
242  VertexType p_vtx;
243 
244  p_vtx = calc(normals[0]);
245  if (p_vtx.dot(normals[0]) < 0) {
246  return false;
247  }
248 
249  p_vtx = calc(normals[1]);
250  if (p_vtx.dot(normals[1]) < 0) {
251  return false;
252  }
253 
254  p_vtx = calc(normals[2]);
255  if (p_vtx.dot(normals[2]) < 0) {
256  return false;
257  }
258 
259  if constexpr (sides > 2) {
260  p_vtx = calc(normals[3]);
261  if (p_vtx.dot(normals[3]) < 0) {
262  return false;
263  }
264  }
265 
266  if constexpr (sides > 3) {
267  p_vtx = calc(normals[4]);
268  if (p_vtx.dot(normals[4]) < 0) {
269  return false;
270  }
271  }
272 
273  if constexpr (sides > 4) {
274  for (size_t i = 5; i <= fr.sides; i++) {
275  const VertexType& normal = normals[i];
276 
277  p_vtx = calc(normal);
278  if (p_vtx.dot(normal) < 0) {
279  return false;
280  }
281  }
282  }
283 
284  return true;
285  };
286 
287  frustumVariants["Nominal, no early ret."] = [](const Box& box,
288  const Frustum3& fr) {
289  const auto& normals = fr.normals();
290  const vertex_array_type fr_vmin = box.min() - fr.origin();
291  const vertex_array_type fr_vmax = box.max() - fr.origin();
292 
293  VertexType p_vtx;
294  bool result = true;
295  for (size_t i = 0; i < fr.sides + 1; i++) {
296  const VertexType& normal = normals[i];
297 
298  p_vtx = (normal.array() < 0).template cast<value_type>() * fr_vmin +
299  (normal.array() >= 0).template cast<value_type>() * fr_vmax;
300 
301  result = result && (p_vtx.dot(normal) >= 0);
302  }
303  return result;
304  };
305 
306  frustumVariants["Manual constexpr unroll, early ret."] =
307  [](const Box& box, const Frustum3& fr) {
308  constexpr size_t sides = 4; // yes this is pointless, I just want to
309  // kind of match the other impl
310 
311  const auto& normals = fr.normals();
312  const vertex_array_type fr_vmin = box.min() - fr.origin();
313  const vertex_array_type fr_vmax = box.max() - fr.origin();
314 
315  auto calc = [&](const auto& normal) {
316  return (normal.array() < 0).template cast<value_type>() * fr_vmin +
317  (normal.array() >= 0).template cast<value_type>() * fr_vmax;
318  };
319 
320  VertexType p_vtx;
321  bool result = true;
322 
323  p_vtx = calc(normals[0]);
324  result = result && (p_vtx.dot(normals[0]) >= 0);
325 
326  p_vtx = calc(normals[1]);
327  result = result && (p_vtx.dot(normals[1]) >= 0);
328 
329  p_vtx = calc(normals[2]);
330  result = result && (p_vtx.dot(normals[2]) >= 0);
331 
332  if constexpr (sides > 2) {
333  p_vtx = calc(normals[3]);
334  result = result && (p_vtx.dot(normals[3]) >= 0);
335  }
336 
337  if constexpr (sides > 3) {
338  p_vtx = calc(normals[4]);
339  result = result && (p_vtx.dot(normals[4]) >= 0);
340  }
341 
342  if constexpr (sides > 4) {
343  for (size_t i = 5; i <= fr.sides; i++) {
344  const VertexType& normal = normals[i];
345 
346  p_vtx = calc(normal);
347  result = result && (p_vtx.dot(normal) >= 0);
348  }
349  }
350 
351  return result;
352  };
353 
354  std::vector<Frustum3> frustums{n, Frustum3{{0, 0, 0}, {1, 0, 0}, M_PI / 2.}};
355  std::generate(frustums.begin(), frustums.end(), [&]() {
356  const Vector3F d{dir(rng), dir(rng), dir(rng)};
357  const Vector3F l{loc(rng), loc(rng), loc(rng)};
358  return Frustum3{l, d.normalized(), ang(rng)};
359  });
360 
361  std::cout << "Make sure frustum implementations are identical" << std::endl;
362  for (const auto& fr : frustums) {
363  std::vector<std::pair<std::string, bool>> results;
364 
365  std::transform(frustumVariants.begin(), frustumVariants.end(),
366  std::back_inserter(results),
367  [&](const auto& p) -> decltype(results)::value_type {
368  const auto& [name, func] = p;
369  return {name, func(testBox, fr)};
370  });
371 
372  bool all = std::all_of(results.begin(), results.end(),
373  [](const auto& r) { return r.second; });
374  bool none = std::none_of(results.begin(), results.end(),
375  [](const auto& r) { return r.second; });
376 
377  if (!all && !none) {
378  std::cerr << "Discrepancy: " << std::endl;
379  for (const auto& [name, result] : results) {
380  std::cerr << " - " << name << ": " << result << std::endl;
381  }
382 
383  testBox.toStream(std::cerr);
384  std::cerr << std::endl;
385  std::cerr << "Frustum: [" << fr.origin().transpose() << "], ["
386  << fr.dir().transpose() << "]" << std::endl;
387  return -1;
388  }
389  }
390  std::cout << "Seems ok" << std::endl;
391 
392  size_t iters_per_run = 1000;
393 
394  std::vector<std::pair<std::string, Frustum3>> testFrusts = {
395  {"away", Frustum3{{0, 0, -10}, {0, 0, -1}, M_PI / 4.}},
396  {"towards", Frustum3{{0, 0, -10}, {0, 0, 1}, M_PI / 4.}},
397  {"left", Frustum3{{0, 0, -10}, {0, 1, 0}, M_PI / 4.}},
398  {"right", Frustum3{{0, 0, -10}, {0, -1, 0}, M_PI / 4.}},
399  {"up", Frustum3{{0, 0, -10}, {1, 0, 0}, M_PI / 4.}},
400  {"down", Frustum3{{0, 0, -10}, {-1, 0, 0}, M_PI / 4.}},
401  };
402 
403  std::cout << "Run benchmarks: " << std::endl;
404 
405  for (const auto& fr_pair : testFrusts) {
406  std::cout << "Frustum '" << fr_pair.first << "'" << std::endl;
407 
408  for (const auto& p : frustumVariants) {
409  // can't capture structured binding, so pair access it is.
410  std::cout << "- Benchmarking variant: '" << p.first << "'" << std::endl;
411  auto bench_result = Acts::Test::microBenchmark(
412  [&]() { return p.second(testBox, fr_pair.second); }, iters_per_run);
413  std::cout << " " << bench_result << std::endl;
414  }
415 
416  std::cout << std::endl;
417  }
418  return 0;
419 }