-
Notifications
You must be signed in to change notification settings - Fork 35
/
ArborXBenchmark_PointClouds.hpp
223 lines (187 loc) · 6.39 KB
/
ArborXBenchmark_PointClouds.hpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
/****************************************************************************
* Copyright (c) 2017-2022 by the ArborX authors *
* All rights reserved. *
* *
* This file is part of the ArborX library. ArborX is *
* distributed under a BSD 3-clause license. For the licensing terms see *
* the LICENSE file in the top-level directory. *
* *
* SPDX-License-Identifier: BSD-3-Clause *
****************************************************************************/
#ifndef ARBORX_BENCHMARK_POINT_CLOUDS_HPP
#define ARBORX_BENCHMARK_POINT_CLOUDS_HPP
#include <ArborX_GeometryTraits.hpp>
#include <Kokkos_Core.hpp>
#include <Kokkos_Random.hpp>
#include <random>
namespace ArborXBenchmark
{
enum class PointCloudType
{
filled_box,
hollow_box,
filled_sphere,
hollow_sphere
};
inline PointCloudType to_point_cloud_enum(std::string const &str)
{
if (str == "filled_box")
return PointCloudType::filled_box;
if (str == "hollow_box")
return PointCloudType::hollow_box;
if (str == "filled_sphere")
return PointCloudType::filled_sphere;
if (str == "hollow_sphere")
return PointCloudType::hollow_sphere;
throw std::runtime_error(str +
" doesn't correspond to any known PointCloudType!");
}
namespace Details
{
template <typename Point>
class PointGenerationFunctor
{
private:
static_assert(ArborX::GeometryTraits::is_point_v<Point>);
static constexpr int DIM = ArborX::GeometryTraits::dimension_v<Point>;
using Coordinate = ArborX::GeometryTraits::coordinate_type_t<Point>;
PointCloudType _point_cloud_type;
public:
PointGenerationFunctor(PointCloudType point_cloud_type)
: _point_cloud_type(point_cloud_type)
{}
// Assume random generates in (-1, 1) interval
template <typename Random>
KOKKOS_FUNCTION auto operator()(int i, Random const &random) const
{
switch (_point_cloud_type)
{
case PointCloudType::filled_box:
return filledBoxPoint(random);
case PointCloudType::hollow_box:
return hollowBoxPoint(i, random);
case PointCloudType::filled_sphere:
return filledSpherePoint(random);
case PointCloudType::hollow_sphere:
return hollowSpherePoint(random);
default:
Kokkos::abort("ArborX: implementation bug");
}
}
private:
template <typename Random>
KOKKOS_FUNCTION auto filledBoxPoint(Random const &random) const
{
Point p;
for (int d = 0; d < DIM; ++d)
p[d] = random();
return p;
}
template <typename Random>
KOKKOS_FUNCTION auto hollowBoxPoint(int i, Random const &random) const
{
Point p;
// For 3D, the order is
// Indices: 0 1 2 3 4 5 6 7 8 9
// Axes : 0- 0+ 1- 1+ 2- 2+ 0- 0+ 1- 1+
int axis = (i / 2) % DIM;
for (int d = 0; d < DIM; ++d)
{
if (d != axis)
p[d] = random();
else if (i % 2 == 0)
p[d] = -1;
else
p[d] = 1;
}
return p;
}
template <typename Random>
KOKKOS_FUNCTION auto filledSpherePoint(Random const &random) const
{
Point p;
auto const radius_squared = 1;
Coordinate norm_squared;
do
{
norm_squared = 0;
for (int d = 0; d < DIM; ++d)
{
p[d] = random();
norm_squared += p[d] * p[d];
}
// Only accept points that are in the sphere
} while (norm_squared > radius_squared);
return p;
}
template <typename Random>
KOKKOS_FUNCTION auto hollowSpherePoint(Random const &random) const
{
Point p;
// We use "Alternative method 2" from
// https://corysimon.github.io/articles/uniformdistn-on-sphere/.
// Note that once the dimensions go high, more and more points are going to
// be rejected, and the algorithm based on normal distribution would be
// more efficient.
auto const radius_squared = 1;
Coordinate norm_squared;
do
{
norm_squared = 0;
for (int d = 0; d < DIM; ++d)
{
p[d] = random();
norm_squared += p[d] * p[d];
}
// Only accept points that are in the sphere
} while (norm_squared > radius_squared || norm_squared == 0);
auto norm = Kokkos::sqrt(norm_squared);
for (int d = 0; d < DIM; ++d)
p[d] /= norm;
return p;
}
};
} // namespace Details
template <typename ExecutionSpace, typename Points>
void generatePointCloud(ExecutionSpace const &exec,
PointCloudType const point_cloud_type,
double const length, Points &random_points)
{
static_assert(Kokkos::is_view_v<Points>);
using Point = typename Points::value_type;
using namespace ArborX::GeometryTraits;
check_valid_geometry_traits(Point{});
static_assert(is_point_v<Point>, "ArborX: View must contain point values");
KOKKOS_ASSERT(point_cloud_type == PointCloudType::filled_box ||
point_cloud_type == PointCloudType::hollow_box ||
point_cloud_type == PointCloudType::filled_sphere ||
point_cloud_type == PointCloudType::hollow_sphere);
static constexpr int DIM = ArborX::GeometryTraits::dimension_v<Point>;
using Coordinate = ArborX::GeometryTraits::coordinate_type_t<Point>;
using GeneratorPool = Kokkos::Random_XorShift1024_Pool<ExecutionSpace>;
using Generator = typename GeneratorPool::generator_type;
constexpr unsigned int batch_size = 8;
GeneratorPool random_pool(0);
Details::PointGenerationFunctor<Point> functor(point_cloud_type);
unsigned int const n = random_points.extent(0);
Kokkos::parallel_for(
"ArborXBenchmark::generatePointCloud::generate",
Kokkos::RangePolicy<ExecutionSpace>(exec, 0, n / batch_size),
KOKKOS_LAMBDA(int i) {
auto generator = random_pool.get_state();
auto random = [&generator]() {
return Kokkos::rand<Generator, Coordinate>::draw(generator, -1, 1);
};
auto begin = i * batch_size;
auto end = Kokkos::min((i + 1) * batch_size, n);
for (unsigned int k = begin; k < end; ++k)
{
random_points(k) = functor(k, random);
for (int d = 0; d < DIM; ++d)
random_points(k)[d] *= length;
}
random_pool.free_state(generator);
});
}
} // namespace ArborXBenchmark
#endif