PACMAN 0.1.0
Portable Algorithms for Coupling, Mapping, and Adaptive iNterpolation
Loading...
Searching...
No Matches
clustering.hpp
Go to the documentation of this file.
1//
2// This file is subject to the terms and conditions defined in
3// file 'LICENSE', which is part of this source code package.
4//
5
6#pragma once
7
8#include "callbacks.hpp"
9#include "common/types.hpp"
10#include "interpolator.hxx"
11#include "utils/operators.hpp"
12
13namespace PACMAN {
14namespace RbfPum {
15
20template <int_t X, int_t N> constexpr index_t CompileTimePow(void) {
21 index_t prod = 1;
22 for (index_t i = 0; i < N; ++i) {
23 prod *= X;
24 }
25 return prod;
26}
27
34template <int_t Dim> constexpr auto ZCurveNeighborOffsets(int_t shape[Dim]) {
35 constexpr index_t N = CompileTimePow<3, Dim>() - 1;
36 std::array<int_t, N> result{};
37
38 std::array<int_t, Dim> offset{};
39 index_t out = 0;
40
41 const index_t total = N + 1;
42
43 for (index_t i = 0; i < total; ++i) {
44 index_t t = i;
45 bool all_zero = true;
46
47 for (int_t d = 0; d < Dim; ++d) {
48 int_t digit = static_cast<int_t>(t % 3);
49 t /= 3;
50 digit -= 1;
51 offset[d] = digit;
52 if (digit != 0) {
53 all_zero = false;
54 }
55 }
56
57 if (!all_zero) {
58 // Indice zCurve du offset calculé
59 int_t index = 0;
60 index_t stride = 1;
61 for (int_t k = Dim - 1; k >= 0; --k) {
62 index += offset[k] * stride;
63 stride *= shape[k];
64 }
65 result[out++] = index;
66 }
67 }
68 return result;
69}
70
72template <> constexpr auto ZCurveNeighborOffsets<2>(int_t shape[2]) {
73 std::array<int_t, 8> result{};
74
75 const std::array<std::array<int_t, 2>, 8> offsets = {
76 {{-1, -1}, {-1, 0}, {-1, 1}, {0, -1}, {0, 1}, {1, -1}, {1, 0}, {1, 1}}};
77
78 index_t out = 0;
79 for (const auto offset : offsets) {
80 int_t index = 0;
81 index_t stride = 1;
82 for (int_t k = 1; k >= 0; --k) {
83 index += offset[k] * stride;
84 stride *= shape[k];
85 }
86 result[out++] = index;
87 }
88
89 return result;
90}
91
93template <> constexpr auto ZCurveNeighborOffsets<3>(int_t shape[3]) {
94 std::array<int_t, 26> result{};
95
96 const std::array<std::array<int_t, 3>, 26> offsets = {
97 {{-1, -1, -1}, {-1, -1, 0}, {-1, -1, 1}, {-1, 0, 0}, {-1, 0, 1},
98 {-1, 0, -1}, {-1, 1, 0}, {-1, 1, 1}, {-1, 1, -1}, {0, -1, -1},
99 {0, -1, 0}, {0, -1, 1}, {0, 0, 1}, {0, 0, -1}, {0, 1, 0},
100 {0, 1, 1}, {0, 1, -1}, {1, -1, -1}, {1, -1, 0}, {1, -1, 1},
101 {1, 0, 0}, {1, 0, 1}, {1, 0, -1}, {1, 1, 0}, {1, 1, 1},
102 {1, 1, -1}}};
103
104 index_t out = 0;
105 for (const auto offset : offsets) {
106 int_t index = 0;
107 index_t stride = 1;
108 for (int_t k = 2; k >= 0; --k) {
109 index += offset[k] * stride;
110 stride *= shape[k];
111 }
112 result[out++] = index;
113 }
114
115 return result;
116}
117
120template <typename ValueType> struct RemoveTaggedCenters {
121 KOKKOS_INLINE_FUNCTION
122 bool operator()(const ValueType &v) const { return v[0] != v[0]; }
123};
124
127FULL_TEMPLATE
128void TEMPLATED_CLASSNAME::CreateClusters(void) {
129 assert(this->mRadius > 0);
130 assert(this->mRelativeOverlap > 0 && this->mRelativeOverlap < 1);
131
132 const std::string _region_name = "RbfPumInterpolator::create_clusters";
133 Kokkos::Profiling::ScopedRegion region(_region_name);
134 const ExecSpace execspace{};
135
136 // 1. computation of the source points bounding box
137 const Point lower = this->mSourceBvh.bounds().minCorner();
138 const Point upper = this->mSourceBvh.bounds().maxCorner();
139
140 // 2. compatation of the max distance between two points to guarantee the
141 // given overlap ratio
142 const fp_t spacing = std::sqrt(4.0 / Dim) * this->mRadius *
143 (fp_consts::one() - this->mRelativeOverlap);
144
145 int_t shape[Dim];
146 fp_t distances[Dim];
147 fp_t min_distance = fp_consts::max();
148 int_t nb_centers = 1;
149 for (int_t axis = 0; axis < Dim; ++axis) {
150 fp_t edge_length = upper[axis] - lower[axis];
151 edge_length = (edge_length < 0) ? (-edge_length) : (edge_length);
152
153 // 3-1. computation of the number of centers we can fit in each direction
154 shape[axis] = std::ceil(std::max(fp_consts::one(), edge_length / spacing));
155
156 // 4. computation of the distances dx, dy, dz, the space between centers in
157 // each direction
158 distances[axis] = edge_length / shape[axis];
159 if (distances[axis] < min_distance) {
160 min_distance = distances[axis];
161 }
162
163 // we begin fitting centers at the beginning of the bounding box
164 nb_centers *= ++shape[axis];
165 }
166
167 // 5. we create a view with all of the clusters centers we want to try
168 this->mClusters = decltype(this->mClusters)(
169 Kokkos::view_alloc(execspace, Kokkos::WithoutInitializing,
170 "this->mClusters"),
171 nb_centers);
172 auto centers_candidates = this->mClusters;
173
174 Kokkos::parallel_for(
175 _region_name + "::p_for fill centers candidates",
176 Kokkos::RangePolicy(execspace, 0, nb_centers),
177 KOKKOS_LAMBDA(const index_t &i) {
178 index_t i_copy = i;
179 // a. we compute the coordinates of the point on the centers grid
180 Kokkos::Array<int_t, Dim> coords{};
181 for (int_t axis = Dim - 1; axis >= 0; --axis) {
182 int_t prod = 1;
183 for (int_t d = axis - 1; d >= 0; --d) {
184 prod *= shape[d];
185 }
186 coords[axis] = i_copy / prod;
187 i_copy -= coords[axis] * prod;
188 }
189
190 // b. computation of the zCurve index
191 int_t index = 0;
192 int_t stride = 1;
193 for (int_t k = Dim - 1; k >= 0; --k) {
194 index += coords[k] * stride;
195 stride *= shape[k];
196 }
197
198 // c. computation of the spatial coordinates of the point in space
199 centers_candidates(index) = Point{lower};
200 for (int_t axis = 0; axis < Dim; ++axis) {
201 centers_candidates(index)[axis] += distances[axis] * coords[axis];
202 }
203 });
204 Kokkos::fence();
205
206 // 6. we tag for removal clusters which do not contain source points
207 // = all clusters with norm2(nearest(center), center) > radius
208 TagEmptyCenters tag_empty_centers{centers_candidates};
209 TagEmptyCentersCallback tag_empty_centers_callback{
210 centers_candidates, this->mRadius * this->mRadius};
211 this->mSourceBvh.query(execspace, tag_empty_centers,
212 tag_empty_centers_callback);
213
214 // 7. same step but checking for target points in clusters
215 this->mTargetBvh.query(execspace, tag_empty_centers,
216 tag_empty_centers_callback);
217
218 // 8. we transform each center candidate to its nearest point on the source
219 // mesh/source points cloud
220 TransformToNearest transform_to_nearest{centers_candidates};
221 TransformToNearestCallback transform_to_nearest_callback{centers_candidates};
222
223 this->mSourceBvh.query(execspace, transform_to_nearest,
224 transform_to_nearest_callback);
225
226 // 9. we remove centers that are too close from each others (can be improved
227 // and taken to parallel computations on GPU) threshold = 0.4 * min(distances)
228 fp_t threshold = 0.4 * min_distance;
229 threshold *= threshold;
230 auto centers_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{},
231 centers_candidates);
232 auto neighbors_offsets = ZCurveNeighborOffsets<Dim>(shape);
233 for (index_t i = 0; i < centers_host.extent(0); ++i) {
234 // is not NaN = is not tagged
235 if (centers_host(i)[0] == centers_host(i)[0]) {
236 for (auto off : neighbors_offsets) {
237 if (!off) {
238 continue;
239 }
240 int_t neighbor_id = i + off;
241 if (neighbor_id >= 0 && neighbor_id < centers_host.extent_int(0)) {
242 auto center = centers_host(i);
243 auto neighbor = centers_host(neighbor_id);
244 // is not NaN = is not tagged AND is too close
245 if (neighbor[0] == neighbor[0] &&
246 SquaredDifference(center, neighbor) < threshold) {
247 centers_host(i)[0] = NAN;
248 break;
249 }
250 }
251 }
252 }
253 }
254 Kokkos::deep_copy(centers_candidates, centers_host);
255
256 // 10. we tag for removal the clusters which do not contain any target point:
257 // = all clusters with norm2(nearest(center), center) > radius
258 tag_empty_centers = decltype(tag_empty_centers){centers_candidates};
259 tag_empty_centers_callback = decltype(tag_empty_centers_callback){
260 centers_candidates, this->mRadius * this->mRadius};
261 this->mTargetBvh.query(execspace, tag_empty_centers,
262 tag_empty_centers_callback);
263
264 // 11. we remove tagged centers
265 const auto end = Kokkos::Experimental::remove_if(
266 _region_name + "::remove_if", execspace, centers_candidates,
267 RemoveTaggedCenters<Point>{});
268 const int_t dist = Kokkos::Experimental::distance(
269 Kokkos::Experimental::begin(centers_candidates), end);
270 Kokkos::resize(this->mClusters, dist);
271}
272} // namespace RbfPum
273
274} // namespace PACMAN
PACMAN::RbfPum::ZCurveNeighborOffsets< 2 >
constexpr auto ZCurveNeighborOffsets< 2 >(int_t shape[2])
: 2D specialization of ZCurveNeighborOffsets for faster computations
Definition clustering.hpp:72
PACMAN::RbfPum::CompileTimePow
constexpr index_t CompileTimePow(void)
A pow function computed at compile time.
Definition clustering.hpp:20
PACMAN::RbfPum::SquaredDifference
KOKKOS_INLINE_FUNCTION fp_t SquaredDifference(const ::ArborX::Point< Dim, coordinates_t > &lhs, const ::ArborX::Point< Dim, coordinates_t > &rhs) noexcept
Returns the squared distance between two points, and check for NaN values.
Definition operators.hpp:127
PACMAN::RbfPum::ZCurveNeighborOffsets
constexpr auto ZCurveNeighborOffsets(int_t shape[Dim])
Return the indices offsets of the spatial neighbors given the shape of the space.
Definition clustering.hpp:34
PACMAN::RbfPum::ZCurveNeighborOffsets< 3 >
constexpr auto ZCurveNeighborOffsets< 3 >(int_t shape[3])
: 3D specialization of ZCurveNeighborOffsets for faster computations
Definition clustering.hpp:93
PACMAN::fp_consts::max
KOKKOS_INLINE_FUNCTION fp_t max(void)
Definition types.hpp:89
PACMAN::fp_consts::one
KOKKOS_INLINE_FUNCTION consteval fp_t one(void)
Definition types.hpp:95
PACMAN::fp_t
double fp_t
Definition types.hpp:15
PACMAN::index_t
uint32_t index_t
Definition types.hpp:20
PACMAN::int_t
int32_t int_t
Definition types.hpp:18
PACMAN::RbfPum::RemoveTaggedCenters
A functor which returns true is a given point is tagged for removal (=> if its first dimension coordi...
Definition clustering.hpp:120
PACMAN::RbfPum::RemoveTaggedCenters::operator()
KOKKOS_INLINE_FUNCTION bool operator()(const ValueType &v) const
Definition clustering.hpp:122
PACMAN::is_kokkos_view
Definition concepts.hpp:17