utils/grids/grid.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 // Detray include(s)
0010 #include "detray/utils/grid/grid.hpp"
0011
0012 #include "detray/definitions/indexing.hpp"
0013 #include "detray/geometry/shapes/cuboid3D.hpp"
0014 #include "detray/utils/grid/concepts.hpp"
0015
0016 // Detray test include(s)
0017 #include "detray/test/framework/types.hpp"
0018
0019 // Vecmem include(s)
0020 #include <vecmem/containers/device_vector.hpp>
0021 #include <vecmem/containers/vector.hpp>
0022 #include <vecmem/memory/host_memory_resource.hpp>
0023
0024 // GTest include(s)
0025 #include <gtest/gtest.h>
0026
0027 // System include(s)
0028 #include <algorithm>
0029 #include <limits>
0030 #include <random>
0031
0032 using namespace detray;
0033 using namespace detray::axis;
0034
0035 namespace {
0036
0037 // Algebra definitions
0038 using test_algebra = test::algebra;
0039 using scalar = test::scalar;
0040 using point3 = test::point3;
0041
0042 constexpr scalar inf{std::numeric_limits<scalar>::max()};
0043 constexpr scalar tol{1e-7f};
0044
0045 // Either a data owning or non-owning 3D cartesian multi-axis
0046 template <bool ownership = true, typename containers = host_container_types>
0047 using cartesian_3D =
0048     coordinate_axes<axes<cuboid3D>, test_algebra, ownership, containers>;
0049
0050 // non-owning multi-axis: Takes external containers
0051 bool constexpr is_owning = true;
0052 bool constexpr is_n_owning = false;
0053
0054 // Bin edges for all axes
0055 dvector<scalar> bin_edges = {-10.f, 10.f, -20.f, 20.f, 0.f, 100.f};
0056 // Offsets into edges container and #bins for all axes
0057 dvector<dsized_index_range> edge_ranges = {{0u, 20u}, {2u, 40u}, {4u, 50u}};
0058
0059 // non-owning multi-axis for the non-owning grid
0060 cartesian_3D<is_n_owning, host_container_types> ax_n_own(edge_ranges,
0061                                                          bin_edges);
0062
0063 // Create some bin data for non-owning grid
0064 template <typename populator_t, typename bin_t>
0065 struct bin_content_sequence {
0066   using entry_t = typename bin_t::entry_type;
0067   entry_t entry{0};
0068
0069   auto operator()() {
0070     entry += entry_t{1};
0071     bin_t bin{};
0072     populator_t{}(bin, entry);
0073     return bin;
0074   }
0075 };
0076
0077 /// Test bin content element by element
0078 template <concepts::grid grid_t, typename content_t>
0079 void test_content(const grid_t& g, const point3& p, const content_t& expected) {
0080   dindex i = 0u;
0081   for (const auto& entry : g.bin(p)) {
0082     ASSERT_NEAR(entry, expected[i++], tol) << " at index " << i - 1u;
0083   }
0084 }
0085
0086 }  // anonymous namespace
0087
0088 /// Unittest: Test single-bin grid construction
0089 GTEST_TEST(detray_grid, single_grid) {
0090   // Owning and non-owning, cartesian, 3-dimensional grids
0091   using grid_owning_t =
0092       grid<test_algebra, axes<cuboid3D>, bins::single<scalar>>;
0093
0094   using grid_n_owning_t =
0095       grid<test_algebra, axes<cuboid3D>, bins::single<scalar>,
0096            simple_serializer, host_container_types, false>;
0097
0098   using grid_device_t = grid<test_algebra, axes<cuboid3D>, bins::single<scalar>,
0099                              simple_serializer, device_container_types>;
0100
0101   static_assert(concepts::grid<grid_owning_t>);
0102   static_assert(concepts::grid<grid_n_owning_t>);
0103   static_assert(concepts::grid<grid_device_t>);
0104
0105   // Fill the bin data for every test
0106   // bin test entries
0107   grid_owning_t::bin_container_type bin_data{};
0108   bin_data.resize(40'000u);
0109   std::ranges::generate_n(
0110       bin_data.begin(), 40'000u,
0111       bin_content_sequence<replace<>, bins::single<scalar>>());
0112
0113   // Copy data that will be moved into the data owning types
0114   dvector<scalar> bin_edges_cp(bin_edges);
0115   dvector<dsized_index_range> edge_ranges_cp(edge_ranges);
0116   grid_owning_t::bin_container_type bin_data_cp(bin_data);
0117
0118   // Data-owning axes and grid
0119   cartesian_3D<is_owning, host_container_types> axes_own(
0120       std::move(edge_ranges_cp), std::move(bin_edges_cp));
0121   grid_owning_t grid_own(std::move(bin_data_cp), std::move(axes_own));
0122
0123   // Copy a second time for the comparison
0124   dvector<scalar> bin_edges_cp2(bin_edges);
0125   dvector<dsized_index_range> edge_ranges_cp2(edge_ranges);
0126   grid_owning_t::bin_container_type bin_data_cp2(bin_data);
0127
0128   // Make a second grid
0129   cartesian_3D<is_owning, host_container_types> axes_own2(
0130       std::move(edge_ranges_cp2), std::move(bin_edges_cp2));
0131   grid_owning_t grid_own2(std::move(bin_data_cp2), std::move(axes_own2));
0132
0133   // CHECK equality
0134   EXPECT_TRUE(grid_own == grid_own2);
0135
0136   // Check a few basics
0137   EXPECT_EQ(grid_own.dim, 3u);
0138   EXPECT_EQ(grid_own.nbins(), 40'000u);
0139   auto y_axis = grid_own.get_axis<label::e_y>();
0140   EXPECT_EQ(y_axis.nbins(), 40u);
0141   auto z_axis =
0142       grid_own.get_axis<single_axis<closed<label::e_z>, regular<scalar>>>();
0143   EXPECT_EQ(z_axis.nbins(), 50u);
0144
0145   // Create non-owning grid
0146   grid_n_owning_t grid_n_own(&bin_data, ax_n_own);
0147
0148   // Test for consistency with owning grid
0149   EXPECT_EQ(grid_n_own.dim, grid_own.dim);
0150   y_axis = grid_n_own.get_axis<label::e_y>();
0151   EXPECT_EQ(y_axis.nbins(), grid_own.get_axis<label::e_y>().nbins());
0152   z_axis = grid_n_own.get_axis<label::e_z>();
0153   EXPECT_EQ(z_axis.nbins(), grid_own.get_axis<label::e_z>().nbins());
0154
0155   // Construct a grid from a view
0156   grid_owning_t::view_type grid_view = get_data(grid_own);
0157   grid_device_t device_grid(grid_view);
0158
0159   // Test for consistency with non-owning grid
0160   EXPECT_EQ(device_grid.dim, grid_n_own.dim);
0161   auto y_axis_dev = device_grid.get_axis<label::e_y>();
0162   EXPECT_EQ(y_axis_dev.nbins(), grid_n_own.get_axis<label::e_y>().nbins());
0163   auto z_axis_dev = device_grid.get_axis<label::e_z>();
0164   EXPECT_EQ(z_axis_dev.nbins(), grid_n_own.get_axis<label::e_z>().nbins());
0165
0166   // Test the global bin iteration: owning grid
0167   auto seq = detray::views::iota(1, 40'001);
0168   auto flat_bin_view = grid_own.all();
0169
0170   static_assert(detray::ranges::random_access_range<decltype(flat_bin_view)>);
0171
0172   EXPECT_EQ(seq.size(), 40'000u);
0173   EXPECT_EQ(flat_bin_view.size(), 40'000u);
0174   EXPECT_EQ(flat_bin_view[42], 43u);
0175   EXPECT_TRUE(
0176       std::equal(flat_bin_view.begin(), flat_bin_view.end(), seq.begin()));
0177
0178   // Test the global bin iteration: non-owning grid
0179   auto flat_bin_view2 = grid_n_own.all();
0180
0181   static_assert(detray::ranges::random_access_range<decltype(flat_bin_view2)>);
0182
0183   EXPECT_EQ(seq.size(), 40'000u);
0184   EXPECT_EQ(flat_bin_view2.size(), 40'000u);
0185   EXPECT_EQ(flat_bin_view2[42], 43u);
0186   EXPECT_TRUE(
0187       std::equal(flat_bin_view2.begin(), flat_bin_view2.end(), seq.begin()));
0188
0189   // Test const grid view
0190   /*auto const_grid_view = get_data(const_cast<const
0191   grid_owning_t&>(grid_own));
0192
0193   static_assert(
0194       std::is_same_v<decltype(const_grid_view),
0195                      typename grid<test_algebra,cartesian_3D<>,
0196   bins::single<const scalar>>::view_type>, "Const grid view was not correctly
0197   constructed!");
0198
0199   grid<test_algebra,cartesian_3D<is_owning, device_container_types>,
0200   bins::single<const scalar>> const_device_grid(const_grid_view);
0201
0202   static_assert(
0203       std::is_same_v<typename decltype(const_device_grid)::bin_type,
0204                      typename replace::template bin_type<const scalar>>,
0205       "Const grid was not correctly constructed from view!");*/
0206 }
0207
0208 /// Unittest: Test dynamic array grid
0209 GTEST_TEST(detray_grid, dynamic_array) {
0210   // Owning and non-owning, cartesian, 3-dimensional grids
0211   using grid_owning_t =
0212       grid<test_algebra, axes<cuboid3D>, bins::dynamic_array<scalar>>;
0213
0214   using grid_n_owning_t =
0215       grid<test_algebra, axes<cuboid3D>, bins::dynamic_array<scalar>,
0216            simple_serializer, host_container_types, false>;
0217
0218   using grid_device_t =
0219       grid<test_algebra, axes<cuboid3D>, bins::dynamic_array<scalar>,
0220            simple_serializer, device_container_types>;
0221
0222   static_assert(concepts::grid<grid_owning_t>);
0223   static_assert(concepts::grid<grid_n_owning_t>);
0224   static_assert(concepts::grid<grid_device_t>);
0225
0226   // Fill the bin data for every test
0227   // bin test entries
0228   grid_owning_t::bin_container_type bin_data{};
0229   // 40 000 entries
0230   bin_data.entries.resize(80'000u);
0231   // 20 000 bins
0232   bin_data.bins.resize(40'000u);
0233
0234   int i{0};
0235   dindex offset{0u};
0236   scalar entry{0.f};
0237   complete<> completer{};
0238
0239   // Test data to compare bin content against
0240   std::vector<scalar> seq;
0241   seq.reserve(80'000);
0242
0243   for (auto& data : bin_data.bins) {
0244     data.offset = offset;
0245     // Every second bin holds one element, otherwise three
0246     data.capacity = (i % 2) != 0u ? 1u : 3u;
0247
0248     detray::bins::dynamic_array bin{bin_data.entries.data(), data};
0249
0250     ASSERT_TRUE(bin.capacity() == (i % 2 ? 1u : 3u));
0251     ASSERT_TRUE(bin.size() == 0);
0252
0253     offset += bin.capacity();
0254
0255     // Populate the bin
0256     completer(bin, entry);
0257
0258     for (auto e : bin) {
0259       ASSERT_TRUE(e == entry);
0260       seq.push_back(e);
0261     }
0262     entry += 1.f;
0263     ++i;
0264   }
0265
0266   // Copy data that will be moved into the data owning types
0267   dvector<scalar> bin_edges_cp(bin_edges);
0268   dvector<dsized_index_range> edge_ranges_cp(edge_ranges);
0269   grid_owning_t::bin_container_type bin_data_cp(bin_data);
0270
0271   // Data-owning axes and grid
0272   cartesian_3D<is_owning, host_container_types> axes_own(
0273       std::move(edge_ranges_cp), std::move(bin_edges_cp));
0274   grid_owning_t grid_own(std::move(bin_data_cp), std::move(axes_own));
0275
0276   // Check a few basics
0277   EXPECT_EQ(grid_own.dim, 3u);
0278   EXPECT_EQ(grid_own.nbins(), 40'000u);
0279   auto y_axis = grid_own.get_axis<label::e_y>();
0280   EXPECT_EQ(y_axis.nbins(), 40u);
0281   auto z_axis =
0282       grid_own.get_axis<single_axis<closed<label::e_z>, regular<scalar>>>();
0283   EXPECT_EQ(z_axis.nbins(), 50u);
0284
0285   // Check equality operator:
0286   // - Copy a second time for the comparison
0287   dvector<scalar> bin_edges_cp2(bin_edges);
0288   dvector<dsized_index_range> edge_ranges_cp2(edge_ranges);
0289   grid_owning_t::bin_container_type bin_data_cp2(bin_data);
0290
0291   // Make a second grid
0292   cartesian_3D<is_owning, host_container_types> axes_own2(
0293       std::move(edge_ranges_cp2), std::move(bin_edges_cp2));
0294
0295   grid_owning_t grid_own2(std::move(bin_data_cp2), std::move(axes_own2));
0296
0297   // CHECK equality
0298   EXPECT_TRUE(grid_own == grid_own2);
0299
0300   // Create non-owning grid
0301   grid_n_owning_t grid_n_own(&bin_data, ax_n_own);
0302
0303   // Test for consistency with owning grid
0304   EXPECT_EQ(grid_n_own.dim, grid_own.dim);
0305   y_axis = grid_n_own.get_axis<label::e_y>();
0306   EXPECT_EQ(y_axis.nbins(), grid_own.get_axis<label::e_y>().nbins());
0307   z_axis = grid_n_own.get_axis<label::e_z>();
0308   EXPECT_EQ(z_axis.nbins(), grid_own.get_axis<label::e_z>().nbins());
0309
0310   // Construct a grid from a view
0311   grid_owning_t::view_type grid_view = get_data(grid_own);
0312   grid_device_t device_grid(grid_view);
0313
0314   // Test for consistency with non-owning grid
0315   EXPECT_EQ(device_grid.dim, grid_n_own.dim);
0316   auto y_axis_dev = device_grid.get_axis<label::e_y>();
0317   EXPECT_EQ(y_axis_dev.nbins(), grid_n_own.get_axis<label::e_y>().nbins());
0318   auto z_axis_dev = device_grid.get_axis<label::e_z>();
0319   EXPECT_EQ(z_axis_dev.nbins(), grid_n_own.get_axis<label::e_z>().nbins());
0320
0321   // Test the global bin iteration
0322   auto flat_bin_view = grid_own.all();
0323
0324   static_assert(detray::ranges::bidirectional_range<decltype(flat_bin_view)>);
0325   // TODO: Const-correctness issue
0326   // static_assert(detray::ranges::random_access_range<typename decltype(
0327   //                  flat_bin_view)>);
0328
0329   EXPECT_EQ(seq.size(), 80'000u);
0330   EXPECT_EQ(flat_bin_view.size(), 80'000u);
0331   EXPECT_TRUE(
0332       std::equal(flat_bin_view.begin(), flat_bin_view.end(), seq.begin()));
0333
0334   // Test the global bin iteration: non-owning grid
0335   auto flat_bin_view2 = grid_n_own.all();
0336
0337   static_assert(detray::ranges::bidirectional_range<decltype(flat_bin_view2)>);
0338
0339   EXPECT_EQ(seq.size(), 80'000u);
0340   EXPECT_EQ(flat_bin_view2.size(), 80'000u);
0341   EXPECT_TRUE(
0342       std::equal(flat_bin_view2.begin(), flat_bin_view2.end(), seq.begin()));
0343 }
0344
0345 /// Integration test: Test replace population
0346 GTEST_TEST(detray_grid, replace_population) {
0347   // Non-owning, 3D cartesian  grid
0348   using grid_t = grid<test_algebra, decltype(ax_n_own), bins::single<scalar>,
0349                       simple_serializer, host_container_types, false>;
0350
0351   static_assert(concepts::grid<grid_t>);
0352
0353   // init
0354   using bin_t = grid_t::bin_type;
0355   grid_t::bin_container_type bin_data{};
0356   bin_data.resize(40'000u, bin_t{});
0357
0358   // Create non-owning grid
0359   grid_t g3r(&bin_data, ax_n_own);
0360
0361   // Test the initialization
0362   point3 p = {-10.f, -20.f, 0.f};
0363   for (int ib0 = 0; ib0 < 20; ++ib0) {
0364     for (int ib1 = 0; ib1 < 40; ++ib1) {
0365       for (int ib2 = 0; ib2 < 100; ib2 += 2) {
0366         p = {static_cast<scalar>(-10 + ib0), static_cast<scalar>(-20 + ib1),
0367              static_cast<scalar>(0 + ib2)};
0368         EXPECT_NEAR(g3r.bin(p)[0], std::numeric_limits<scalar>::max(), tol);
0369       }
0370     }
0371   }
0372
0373   p = {-4.5f, -4.5f, 4.5f};
0374   // Fill and read
0375   g3r.template populate<replace<>>(p, 3.f);
0376   EXPECT_NEAR(g3r.bin(p)[0], static_cast<scalar>(3u), tol);
0377
0378   // Fill and read two times, fill first 0-99, then 100-199
0379   for (unsigned int il = 0u; il < 2u; ++il) {
0380     scalar counter{static_cast<scalar>(il * 100u)};
0381     for (int ib0 = 0; ib0 < 20; ++ib0) {
0382       for (int ib1 = 0; ib1 < 40; ++ib1) {
0383         for (int ib2 = 0; ib2 < 100; ib2 += 2) {
0384           p = {static_cast<scalar>(-10 + ib0), static_cast<scalar>(-20 + ib1),
0385                static_cast<scalar>(0 + ib2)};
0386           g3r.template populate<replace<>>(p, counter);
0387           EXPECT_NEAR(g3r.bin(p)[0], counter, tol);
0388           counter += 1.f;
0389         }
0390       }
0391     }
0392   }
0393 }
0394
0395 /// Test bin entry retrieval
0396 GTEST_TEST(detray_grid, complete_population) {
0397   // Non-owning, 3D cartesian, completing grid (4 dims and sort)
0398   using grid_t =
0399       grid<test_algebra, decltype(ax_n_own), bins::static_array<scalar, 4>,
0400            simple_serializer, host_container_types, false>;
0401   using bin_t = grid_t::bin_type;
0402   using bin_content_t = std::array<scalar, 4>;
0403
0404   static_assert(concepts::grid<grid_t>);
0405
0406   // init
0407   grid_t::bin_container_type bin_data{};
0408   bin_data.resize(40'000u, bin_t{});
0409   // Create non-owning grid
0410   grid_t g3c(&bin_data, ax_n_own);
0411
0412   // Test the initialization
0413   point3 p = {-10.f, -20.f, 0.f};
0414   bin_t invalid{};
0415   for (int ib0 = 0; ib0 < 20; ++ib0) {
0416     for (int ib1 = 0; ib1 < 40; ++ib1) {
0417       for (int ib2 = 0; ib2 < 100; ib2 += 2) {
0418         p = {static_cast<scalar>(-10 + ib0), static_cast<scalar>(-20 + ib1),
0419              static_cast<scalar>(0 + ib2)};
0420         test_content(g3c, p, invalid);
0421       }
0422     }
0423   }
0424
0425   p = {-4.5f, -4.5f, 4.5f};
0426   bin_content_t expected{4.f, 4.f, 4.f, 4.f};
0427   // Fill and read
0428   g3c.template populate<complete<>>(p, 4.f);
0429   test_content(g3c, p, expected);
0430
0431   // Test search without search window
0432   for (scalar entry : g3c.bin(p)) {
0433     EXPECT_EQ(entry, 4.f);
0434   }
0435 }
0436
0437 /// Test bin entry retrieval
0438 GTEST_TEST(detray_grid, regular_attach_population) {
0439   // Non-owning, 3D cartesian, completing grid (4 dims and sort)
0440   using grid_t =
0441       grid<test_algebra, decltype(ax_n_own), bins::static_array<scalar, 4>,
0442            simple_serializer, host_container_types, false>;
0443   using bin_t = grid_t::bin_type;
0444   using bin_content_t = std::array<scalar, 4>;
0445
0446   static_assert(concepts::grid<grid_t>);
0447
0448   // init
0449   grid_t::bin_container_type bin_data{};
0450   bin_data.resize(40'000u, bin_t{});
0451
0452   // Create non-owning grid
0453   grid_t g3ra(&bin_data, ax_n_own);
0454
0455   // Test the initialization
0456   point3 p = {-10.f, -20.f, 0.f};
0457   bin_t invalid{};
0458   for (int ib0 = 0; ib0 < 20; ++ib0) {
0459     for (int ib1 = 0; ib1 < 40; ++ib1) {
0460       for (int ib2 = 0; ib2 < 100; ib2 += 2) {
0461         p = {static_cast<scalar>(-10 + ib0), static_cast<scalar>(-20 + ib1),
0462              static_cast<scalar>(0 + ib2)};
0463         test_content(g3ra, p, invalid);
0464       }
0465     }
0466   }
0467
0468   p = {-4.5f, -4.5f, 4.5f};
0469   bin_content_t expected{5.f, inf, inf, inf};
0470   // Fill and read
0471   g3ra.template populate<attach<>>(p, 5.f);
0472   test_content(g3ra, p, expected);
0473
0474   // Test search without search window
0475   for (scalar entry : g3ra.bin(p)) {
0476     EXPECT_EQ(entry, 5.f);
0477   }
0478 }