BioCMAMC-ST
particles_container.hpp
1#ifndef __PARTICLES_CONTAINER_HPP__
2#define __PARTICLES_CONTAINER_HPP__
3
4#include "Kokkos_Core_fwd.hpp"
5#include "common/common.hpp"
6#include "mc/alias.hpp"
7#include <Kokkos_Core.hpp>
8#include <common/has_serialize.hpp>
9#include <cstdint>
10#include <mc/prng/prng.hpp>
11#include <mc/traits.hpp>
12#include <stdexcept>
13#include <type_traits>
14
15namespace
16{
17
18 using TeamPolicy = Kokkos::TeamPolicy<ComputeSpace>;
19 using TeamMember = TeamPolicy::member_type;
20
21 template <ModelType M> struct FillGapFunctor
22 {
23
24 FillGapFunctor(MC::ParticleStatus _status,
25 M::SelfParticle _model,
26 MC::ParticlePositions _position,
27 MC::ParticleAges _ages,
28 std::size_t _to_remove,
29 std::size_t _last_used_index)
30 : status(std::move(_status)), model(std::move(_model)), position(std::move(_position)),
31 ages(std::move(_ages)), offset("offset"), to_remove(_to_remove),
32 last_used_index(_last_used_index)
33 {
34
35 Kokkos::deep_copy(offset, 0);
36 }
37
38 KOKKOS_INLINE_FUNCTION void operator()(const int i, std::size_t& update, const bool final) const
39 {
40 const bool is_dead = (status(i) != MC::Status::Idle);
41
42 std::size_t scan_index = update;
43 update += is_dead ? 1 : 0;
44
45 if (final && is_dead && scan_index < to_remove)
46 {
47
48 const auto i_to_remove = i;
49 auto idx_to_move = last_used_index - Kokkos::atomic_fetch_add(&offset(), 1);
50 while (status(idx_to_move) != MC::Status::Idle ||
51 idx_to_move == static_cast<std::size_t>(i_to_remove))
52 {
53 idx_to_move = last_used_index - Kokkos::atomic_fetch_add(&offset(), 1);
54 }
55
56 status(i_to_remove) = MC::Status::Idle;
57 Kokkos::atomic_exchange(&position(i_to_remove), position(idx_to_move));
58 // TODO Use hierachical parallism here, thread range is likely to work
59 for (std::size_t i_properties = 0; i_properties < M::n_var; ++i_properties)
60 {
61 model(i_to_remove, i_properties) = model(idx_to_move, i_properties);
62 }
63 ages(i_to_remove, 0) = ages(idx_to_move, 0);
64 ages(i_to_remove, 1) = ages(idx_to_move, 1);
65 }
66 }
67
68 MC::ParticleStatus status;
69 M::SelfParticle model;
70 MC::ParticlePositions position;
72 Kokkos::View<std::size_t, ComputeSpace> offset;
73 std::size_t to_remove;
74 std::size_t last_used_index;
75 };
76
77 template <ModelType M> struct InsertFunctor
78 {
79 InsertFunctor(std::size_t _original_size,
80 M::SelfParticle _model,
81 MC::ParticlePositions _position,
82 M::SelfParticle _buffer_model,
83 MC::ParticlePositions _buffer_position)
84 : original_size(_original_size), model(std::move(_model)), position(std::move(_position)),
85 buffer_model(std::move(_buffer_model)), buffer_position(std::move(_buffer_position))
86 {
87 }
88
89 // KOKKOS_INLINE_FUNCTION void operator()(const int i) const
90 // {
91 // position(original_size + i) = buffer_position(i);
92 // for (std::size_t j = 0; j < M::n_var; ++j)
93 // {
94 // model(original_size + i, j) = buffer_model(i, j);
95 // }
96 // }
97
98 // TODO try this functor
99 // TODO find a way to organise data to not copy non needed data (like contribs). Split model in
100 // two arrays?
101
102 KOKKOS_INLINE_FUNCTION
103 void operator()(const TeamMember& team) const
104 {
105 auto range = M::n_var;
106 const int i = team.league_rank();
107
108 position(original_size + i) = buffer_position(i);
109 Kokkos::parallel_for(Kokkos::TeamVectorRange(team, range),
110 [&](const int& j) { model(original_size + i, j) = buffer_model(i, j); });
111 }
112
113 std::size_t original_size;
114 M::SelfParticle model;
115 MC::ParticlePositions position;
116 M::SelfParticle buffer_model;
117 MC::ParticlePositions buffer_position;
118 };
119
120}; // namespace
121
122namespace MC
123{
134 template <ModelType Model> class ParticlesContainer
135 {
136 public:
137 static constexpr double buffer_ratio =
138 1;
140
143 using UsedModel = Model;
144
145 explicit ParticlesContainer(std::size_t n_particle);
146 ParticlesContainer(); //=default;
147
148 Model::SelfParticle model;
160
165
171 [[nodiscard]] KOKKOS_INLINE_FUNCTION std::size_t n_particles() const
172 {
173 return n_used_elements;
174 };
175
176 // TODO Add logic off internal_dead_counter header
177 std::size_t counter = 0; //FIXME !!
178
179 KOKKOS_INLINE_FUNCTION void get_contributions(std::size_t idx,
180 const ContributionView& contributions) const;
181
182 [[nodiscard]] KOKKOS_INLINE_FUNCTION bool
183 handle_division(const MC::KPRNG::pool_type& random_pool, std::size_t idx1) const;
184
185 void clean_dead(std::size_t to_remove);
186
188
189 [[maybe_unused]] [[nodiscard]] auto get_buffer_index() const
190 {
191 return buffer_index();
192 }
193
194 template <class Archive> void save(Archive& ar) const
195 {
197 serialize_view(ar, weights);
198 serialize_view(ar, position);
199 serialize_view(ar, status);
200 serialize_view(ar, model);
201 serialize_view(ar, ages);
202 }
203
204 template <class Archive> void load(Archive& ar)
205 {
207 deserialize_view(ar, weights);
208 deserialize_view(ar, position);
209 deserialize_view(ar, status);
210 deserialize_view(ar, model);
211 deserialize_view(ar, ages);
212
214 }
215
216 [[nodiscard]] KOKKOS_INLINE_FUNCTION double get_weight(std::size_t idx) const;
217
218 [[nodiscard]] double get_allocation_factor() const
219 {
220 return allocation_factor;
221 }
222
223 [[nodiscard]] std::size_t capacity() const
224 {
226 }
227
228 private:
229 Model::SelfParticle buffer_model;
231 Kokkos::View<uint64_t, Kokkos::SharedSpace> buffer_index;
235
236 static constexpr double default_allocation_factor = 2;
237 void __allocate__(std::size_t new_size);
239 void __shrink__(std::size_t new_size, bool force);
240 };
241
242 template <ModelType Model>
243 KOKKOS_INLINE_FUNCTION bool
245 std::size_t idx1) const
246 {
247 if (Kokkos::atomic_load(&buffer_index()) < buffer_model.extent(0))
248 {
249 const auto idx2 = Kokkos::atomic_fetch_add(&buffer_index(), 1);
250 Model::division(random_pool, idx1, idx2, model, buffer_model);
251 buffer_position(idx2) = position(idx1);
252
253 return true;
254 }
255 Kokkos::printf("%ld %ld\r\n", Kokkos::atomic_load(&buffer_index()), buffer_model.extent(0));
256 return false;
257 }
258
259 template <ModelType Model> void ParticlesContainer<Model>::merge_buffer()
260 {
261 PROFILE_SECTION("ParticlesContainer::merge_buffer")
262 const auto original_size = n_used_elements;
263 const auto n_add_item = buffer_index();
264 if (n_add_item == 0)
265 {
266 return;
267 }
268 __allocate__(original_size + n_add_item);
269 // Merge position EZ
270
271 auto get_policy_insert = [=]()
272 {
273 if constexpr (std::is_same_v<Kokkos::DefaultHostExecutionSpace,
274 Kokkos::DefaultExecutionSpace>)
275 {
276 return TeamPolicy(n_add_item / Kokkos::num_threads(), Kokkos::AUTO, Model::n_var);
277 }
278 else
279 {
280 return TeamPolicy(n_add_item, Kokkos::AUTO, Model::n_var);
281 }
282 };
283
284 Kokkos::parallel_for(
285 "InsertMerge",
286 TeamPolicy(n_add_item, Kokkos::AUTO, Model::n_var),
287 InsertFunctor<Model>(original_size, model, position, buffer_model, buffer_position));
288
289 buffer_index() = 0;
290 n_used_elements += n_add_item;
291 __allocate_buffer__();
292
293 }
294
295 //TODO Merge reduce duplicate __allocate__ and __shrink__
296 template <ModelType Model>
297 void ParticlesContainer<Model>::__allocate__(const std::size_t new_size)
298 {
299 PROFILE_SECTION("ParticlesContainer::__allocate__")
300 if (new_size > 0)
301 {
302 if (new_size >= n_allocated_elements)
303 {
304 const auto new_allocated_size =
305 static_cast<std::size_t>(std::ceil(static_cast<double>(new_size) * allocation_factor));
306 n_allocated_elements = new_allocated_size;
307 Kokkos::resize(position, n_allocated_elements);
308 Kokkos::resize(model, n_allocated_elements, Model::n_var); // use 2nd dim resize if dynamic
309 Kokkos::resize(status, n_allocated_elements);
310 Kokkos::resize(ages, n_allocated_elements);
311 if constexpr (ConstWeightModelType<Model>)
312 {
313 Kokkos::resize(weights, 1);
314 }
315 else
316 {
317 Kokkos::resize(weights, n_allocated_elements);
318 }
319 }
320 }
321 }
322
323 template <ModelType Model>
324 void ParticlesContainer<Model>::__shrink__(std::size_t new_size, bool force)
325 {
326 PROFILE_SECTION("ParticlesContainer::__shrink__")
327 if (new_size > 0 && (new_size > n_used_elements || force))
328 {
329 const auto new_allocated_size =
330 static_cast<std::size_t>(std::ceil(static_cast<double>(new_size) * allocation_factor));
331 n_allocated_elements = new_allocated_size;
332 Kokkos::resize(position, n_allocated_elements);
333 Kokkos::resize(model, n_allocated_elements, Model::n_var); // use 2nd dim resize if dynamic
334 Kokkos::resize(status, n_allocated_elements);
335 Kokkos::resize(ages, n_allocated_elements);
336 if constexpr (ConstWeightModelType<Model>)
337 {
338 Kokkos::resize(weights, 1);
339 }
340 else
341 {
342 Kokkos::resize(weights, n_allocated_elements);
343 }
344 }
345 }
346
347 template <ModelType Model> void ParticlesContainer<Model>::__allocate_buffer__()
348 {
349 PROFILE_SECTION("ParticlesContainer::__allocate_buffer__")
350 auto buffer_size = buffer_position.extent(0);
351 if (static_cast<double>(buffer_size) / static_cast<double>(n_allocated_elements) < buffer_ratio)
352 {
353 buffer_size = static_cast<std::size_t>(
354 std::ceil(static_cast<double>(n_allocated_elements) * buffer_ratio));
355
356 // Realloc because not needed to keep buffer as it has been copied
357 Kokkos::realloc(buffer_position, buffer_size);
358 Kokkos::realloc(buffer_model, buffer_size, Model::n_var); // use 2nd dim resize if dynamic
359 buffer_index() = 0;
360 }
361 }
362
363 template <ModelType M>
365 : model(Kokkos::view_alloc(Kokkos::WithoutInitializing, "particle_model"), 0),
366 position(Kokkos::view_alloc(Kokkos::WithoutInitializing, "particle_position"), 0),
367 status(Kokkos::view_alloc(Kokkos::WithoutInitializing, "particle_status"), 0),
368 weights(Kokkos::view_alloc(Kokkos::WithoutInitializing, "particle_weigth"), 0),
369 ages(Kokkos::view_alloc(Kokkos::WithoutInitializing, "particle_age"), 0),
370 buffer_model("buffer_particle_model", 0),
371 buffer_position("buffer_particle_position", 0), // Dont allocate now
372 buffer_index("buffer_index"), allocation_factor(default_allocation_factor),
373 n_allocated_elements(0), n_used_elements(n_particle)
374 {
375
376 __allocate__(n_particle);
378 }
379
380 template <ModelType M>
382 : model(Kokkos::view_alloc(Kokkos::WithoutInitializing, "particle_model"), 0),
383 position(Kokkos::view_alloc(Kokkos::WithoutInitializing, "particle_position"), 0),
384 status(Kokkos::view_alloc(Kokkos::WithoutInitializing, "particle_status"), 0),
385 weights(Kokkos::view_alloc(Kokkos::WithoutInitializing, "particle_weigth"), 0),
386 ages(Kokkos::view_alloc(Kokkos::WithoutInitializing, "particle_age"), 0),
387 buffer_model(Kokkos::view_alloc(Kokkos::WithoutInitializing, "buffer_particle_model"), 0),
388 buffer_position(Kokkos::view_alloc(Kokkos::WithoutInitializing, "buffer_particle_model")),
389 buffer_index("buffer_index"), allocation_factor(default_allocation_factor),
390 n_allocated_elements(0), n_used_elements(0)
391
392 {
393 }
394
395 template <ModelType M> void ParticlesContainer<M>::clean_dead(std::size_t to_remove)
396 {
397 PROFILE_SECTION("ParticlesContainer::remove_dead")
398 if (to_remove == 0)
399 {
400 return;
401 }
402
403 if (to_remove == n_used_elements)
404 {
405 __shrink__(0, true);
406 n_used_elements = 0;
407 }
408 else if (to_remove > n_used_elements)
409 {
410 throw std::runtime_error(
411 "clean_dead: Error in kernel cannot remove more element than existing");
412 }
413 else
414 {
415
416 const auto new_used_item = n_used_elements - to_remove;
417
418 const auto last_used_index = n_used_elements - 1;
419 Kokkos::parallel_scan(
420 "find_and_fill_gap",
421 Kokkos::RangePolicy<ComputeSpace>(0, n_used_elements),
422 FillGapFunctor<M>(status, model, position, ages, to_remove, last_used_index));
423
424 Kokkos::fence();
425 KOKKOS_ASSERT(this->position.extent(0) == n_allocated_elements);
426 KOKKOS_ASSERT(this->model.extent(0) == n_allocated_elements);
427 KOKKOS_ASSERT(this->status.extent(0) == n_allocated_elements);
428 n_used_elements = new_used_item;
429 if (static_cast<double>(n_used_elements) / static_cast<double>(n_allocated_elements) <= 0.1)
430 {
431 // std::cout << "SHRINK" << std::endl;
432 __shrink__(n_used_elements * 2, false);
433 }
434 };
435 }
436
437 template <ModelType M>
438 KOKKOS_INLINE_FUNCTION void
440 const ContributionView& contributions) const
441 {
442 static_assert(ConstWeightModelType<M>, "ModelType: Const apply_weight()");
443 const double weight = get_weight(idx);
444 M::contribution(idx, position(idx), weight, model, contributions);
445 }
446
447 template <ModelType M>
448 [[nodiscard]] KOKKOS_INLINE_FUNCTION double
450 {
451 if constexpr (ConstWeightModelType<M>)
452 {
453 return weights(0);
454 }
455 else
456 {
457 return weights(idx);
458 }
459 }
460
461} // namespace MC
462
463#endif
Kokkos::Random_XorShift1024_Pool< Kokkos::DefaultExecutionSpace > pool_type
Definition prng.hpp:17
Main owning object for Monte-Carlo particles.
Definition particles_container.hpp:135
~ParticlesContainer()=default
Default destructor.
void __allocate__(std::size_t new_size)
Definition particles_container.hpp:297
Model::SelfParticle model
Definition particles_container.hpp:148
KOKKOS_INLINE_FUNCTION void get_contributions(std::size_t idx, const ContributionView &contributions) const
Definition particles_container.hpp:439
KOKKOS_INLINE_FUNCTION std::size_t n_particles() const
Gets the number of particles in the container.
Definition particles_container.hpp:171
ParticlesContainer(std::size_t n_particle)
Definition particles_container.hpp:364
double allocation_factor
Definition particles_container.hpp:232
ParticleWeigths weights
Definition particles_container.hpp:151
double get_allocation_factor() const
Definition particles_container.hpp:218
Kokkos::View< uint64_t, Kokkos::SharedSpace > buffer_index
Definition particles_container.hpp:231
void __allocate_buffer__()
Definition particles_container.hpp:347
ParticlesContainer()
Definition particles_container.hpp:381
KOKKOS_INLINE_FUNCTION bool handle_division(const MC::KPRNG::pool_type &random_pool, std::size_t idx1) const
Definition particles_container.hpp:244
ParticlesContainer & operator=(ParticlesContainer &&)=default
void __shrink__(std::size_t new_size, bool force)
Definition particles_container.hpp:324
uint64_t n_used_elements
Definition particles_container.hpp:234
auto get_buffer_index() const
Definition particles_container.hpp:189
void merge_buffer()
Definition particles_container.hpp:259
Model::SelfParticle buffer_model
Definition particles_container.hpp:229
KOKKOS_INLINE_FUNCTION double get_weight(std::size_t idx) const
Definition particles_container.hpp:449
std::size_t capacity() const
Definition particles_container.hpp:223
ParticlePositions buffer_position
Definition particles_container.hpp:230
MC::ParticlePositions position
Definition particles_container.hpp:149
void save(Archive &ar) const
Definition particles_container.hpp:194
static constexpr double default_allocation_factor
Definition particles_container.hpp:236
static constexpr double buffer_ratio
Definition particles_container.hpp:137
ParticleAges ages
Definition particles_container.hpp:152
std::size_t counter
Definition particles_container.hpp:177
Model UsedModel
Alias for the model used by the container.
Definition particles_container.hpp:143
void clean_dead(std::size_t to_remove)
Definition particles_container.hpp:395
ParticlesContainer & operator=(const ParticlesContainer &)=default
ParticlesContainer(const ParticlesContainer &)=default
Default copy and move constructors and assignment operators.
void load(Archive &ar)
Definition particles_container.hpp:204
ParticlesContainer(ParticlesContainer &&)=default
std::size_t n_allocated_elements
Definition particles_container.hpp:233
MC::ParticleStatus status
Definition particles_container.hpp:150
Concept to check if a model type has uniform_weight
Definition traits.hpp:173
Namespace that contains classes and structures related to Monte Carlo (MC) simulations.
Definition alias.hpp:9
Kokkos::View< double *, ComputeSpace > ParticleWeigths
Definition alias.hpp:21
Kokkos::View< Status *, ComputeSpace > ParticleStatus
Definition alias.hpp:20
Kokkos::View< double *[2], Kokkos::LayoutLeft, ComputeSpace > ParticleAges
Definition alias.hpp:22
Kokkos::View< uint64_t *, ComputeSpace > ParticlePositions
Definition alias.hpp:19
Kokkos::Experimental::ScatterView< double **, Kokkos::LayoutRight > ContributionView
Definition alias.hpp:31