31#include <boost/mpi/collectives/gather.hpp>
32#include <boost/mpi/collectives/reduce.hpp>
33#include <boost/serialization/utility.hpp>
34#include <boost/serialization/vector.hpp>
47 std::vector<std::reference_wrapper<Particle const>>;
56 std::vector<int> m_ids;
58 virtual std::vector<double>
66 operator()(boost::mpi::communicator
const &comm)
const final;
67 std::vector<int>
const &
ids()
const {
return m_ids; }
76template <
class T>
struct shape_impl;
78template <>
struct shape_impl<double> {
79 static std::vector<std::size_t> eval(std::size_t ) {
return {1}; }
81template <
class _, std::
size_t N>
struct shape_impl<
Utils::Vector<_, N>> {
82 static std::vector<std::size_t> eval(std::size_t ) {
return {N}; }
84template <
class T>
struct shape_impl<std::
vector<T>> {
85 static std::vector<std::size_t> eval(std::size_t n_part) {
86 std::vector<std::size_t> ret{n_part};
87 std::ranges::copy(shape_impl<T>::eval(n_part), std::back_inserter(ret));
92template <
class T,
class U>
struct shape_impl<std::pair<T, U>> {
93 static std::vector<std::size_t> eval(std::size_t n_part) {
94 return shape_impl<T>::eval(n_part);
98inline auto get_argsort(boost::mpi::communicator
const &comm,
99 std::vector<int>
const &local_pids,
100 std::vector<int>
const &sorted_pids) {
101 std::vector<unsigned int> argsort{};
103 std::vector<std::vector<int>> global_pids;
104 boost::mpi::gather(comm, local_pids, global_pids, 0);
105 if (comm.rank() == 0) {
106 auto const n_part = sorted_pids.size();
107 std::vector<int> unsorted_pids;
108 unsorted_pids.reserve(n_part);
109 for (
auto const &vec : global_pids) {
110 for (
auto const pid : vec) {
111 unsorted_pids.emplace_back(pid);
115 std::vector<unsigned int> iota(n_part);
116 std::iota(iota.begin(), iota.end(), 0u);
117 argsort.reserve(n_part);
118 auto const pid_begin = std::begin(unsorted_pids);
119 auto const pid_end = std::end(unsorted_pids);
120 for (
auto const pid : sorted_pids) {
121 auto const pid_pos = std::find(pid_begin, pid_end, pid);
123 static_cast<std::size_t
>(std::distance(pid_begin, pid_pos));
124 argsort.emplace_back(iota[i]);
136get_all_particle_positions(boost::mpi::communicator
const &comm,
138 std::vector<int>
const &sorted_pids,
140 bool use_folded_positions =
false) {
141 using pos_type =
decltype(
traits.position(std::declval<Particle>()));
142 std::vector<pos_type> local_positions{};
143 std::vector<int> local_pids{};
144 local_positions.reserve(local_particles.size());
145 local_pids.reserve(local_particles.size());
147 for (
auto const &particle : local_particles) {
148 if (use_folded_positions) {
149 local_positions.emplace_back(
traits.position_folded(particle));
151 local_positions.emplace_back(
traits.position(particle));
153 local_pids.emplace_back(
traits.id(particle));
156 auto const argsort = detail::get_argsort(comm, local_pids, sorted_pids);
158 std::vector<std::vector<pos_type>> global_positions{};
159 global_positions.reserve(
static_cast<std::size_t
>(comm.size()));
160 boost::mpi::gather(comm, local_positions, global_positions, 0);
162 if (comm.rank() != 0) {
163 return std::vector<pos_type>();
166 std::vector<pos_type> global_positions_flattened{};
167 global_positions_flattened.reserve(sorted_pids.size());
168 for (
auto const &vec : global_positions) {
169 for (
auto const &pos : vec) {
170 global_positions_flattened.emplace_back(pos);
173 assert(global_positions_flattened.size() == sorted_pids.size());
175 std::vector<pos_type> positions_sorted{};
176 positions_sorted.reserve(sorted_pids.size());
177 for (
auto const i : argsort) {
178 positions_sorted.emplace_back(global_positions_flattened[i]);
181 return positions_sorted;
202 std::vector<std::size_t>
shape()
const override {
205 return detail::shape_impl<decltype(declval<ObsType>()(
206 declval<ParticleReferenceRange const &>()))>::eval(
ids().size());
209 template <
typename T>
struct is_map : std::false_type {};
210 template <
typename T>
218 std::vector<double> local_traits{};
219 local_traits.reserve(local_particles.size());
221 std::back_inserter(local_traits));
222 std::vector<int> local_pids{};
223 local_pids.reserve(local_particles.size());
225 std::back_inserter(local_pids));
227 std::vector<std::vector<double>> global_traits{};
228 boost::mpi::gather(comm, local_traits, global_traits, 0);
230 auto const argsort = detail::get_argsort(comm, local_pids,
ids());
232 if (comm.rank() != 0) {
237 auto const size = std::accumulate(
238 global_traits.begin(), global_traits.end(), 0u,
239 [](
auto const acc,
auto const &vec) { return acc + vec.size(); });
241 auto const n_dims = size /
ids().size();
243 std::vector<double> global_traits_flattened{};
244 global_traits_flattened.reserve(size);
246 for (
auto const &vec : global_traits) {
247 for (
auto const val : vec) {
248 global_traits_flattened.emplace_back(val);
252 std::vector<double> output{};
253 output.reserve(size);
255 for (
auto const i : argsort) {
256 for (std::size_t j = 0; j < n_dims; ++j) {
257 output.emplace_back(global_traits_flattened[i * n_dims + j]);
262 auto observable = ObsType{};
263 auto const local_result = observable(local_particles);
264 std::remove_const_t<
decltype(local_result)> result{};
265 boost::mpi::reduce(comm, local_result, result, observable, 0);
Vector implementation and trait types for boost qvm interoperability.
Base class for observables.
This class implements an interface to the particle_observables library that implements necessary algo...
std::vector< double > evaluate(boost::mpi::communicator const &comm, ParticleReferenceRange const &local_particles, ParticleObservables::traits< Particle > const &) const override
std::vector< std::size_t > shape() const override
Particle-based observable.
PidObservable(std::vector< int > ids)
std::vector< int > const & ids() const
std::vector< double > operator()(boost::mpi::communicator const &comm) const final
virtual std::vector< double > evaluate(boost::mpi::communicator const &comm, ParticleReferenceRange const &local_particles, const ParticleObservables::traits< Particle > &traits) const =0
std::vector< std::reference_wrapper< Particle const > > ParticleReferenceRange
void flatten(Range const &v, OutputIterator out)
Flatten a range of ranges.
std::vector< T, allocator< T > > vector