ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
core/observables/PidObservable.hpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2016-2022 The ESPResSo project
3 *
4 * This file is part of ESPResSo.
5 *
6 * ESPResSo is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * ESPResSo is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20#pragma once
21
23
24#include "Observable.hpp"
25#include "Particle.hpp"
26#include "ParticleTraits.hpp"
27
28#include <utils/Vector.hpp>
29#include <utils/flatten.hpp>
30
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>
35
36#include <algorithm>
37#include <cassert>
38#include <cstddef>
39#include <iterator>
40#include <type_traits>
41#include <utility>
42#include <vector>
43
44namespace Observables {
45
47 std::vector<std::reference_wrapper<Particle const>>;
48
49/** Particle-based observable.
50 *
51 * Base class for observables extracting raw data from particle subsets and
52 * returning either the data or a statistic derived from it.
53 */
54class PidObservable : virtual public Observable {
55 /** Identifiers of particles measured by this observable */
56 std::vector<int> m_ids;
57
58 virtual std::vector<double>
59 evaluate(boost::mpi::communicator const &comm,
60 ParticleReferenceRange const &local_particles,
62
63public:
64 explicit PidObservable(std::vector<int> ids) : m_ids(std::move(ids)) {}
65 std::vector<double>
66 operator()(boost::mpi::communicator const &comm) const final;
67 std::vector<int> const &ids() const { return m_ids; }
68};
69
70namespace detail {
71/**
72 * Recursive implementation for finding the shape of a given `std::vector` of
73 * types. A vector of extents is constructed starting at
74 * the template specialization for `std::vector<T>`.
75 */
76template <class T> struct shape_impl;
77
78template <> struct shape_impl<double> {
79 static std::vector<std::size_t> eval(std::size_t /* n_part */) { return {1}; }
80};
81template <class _, std::size_t N> struct shape_impl<Utils::Vector<_, N>> {
82 static std::vector<std::size_t> eval(std::size_t /* n_part */) { return {N}; }
83};
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));
88
89 return ret;
90 }
91};
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);
95 }
96};
97
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{};
102
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);
112 }
113 }
114 // get vector of indices that sorts the data vectors
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);
122 auto const i =
123 static_cast<std::size_t>(std::distance(pid_begin, pid_pos));
124 argsort.emplace_back(iota[i]);
125 }
126 }
127 return argsort;
128}
129
130/** Get the positions of all the particles in the system in the order
131 * specified by \a sorted_pids. Only the head node returns a vector
132 * of positions, all other nodes return an empty vector.
133 * This requires several MPI communications to construct.
134 */
135inline auto
136get_all_particle_positions(boost::mpi::communicator const &comm,
137 ParticleReferenceRange const &local_particles,
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());
146
147 for (auto const &particle : local_particles) {
148 if (use_folded_positions) {
149 local_positions.emplace_back(traits.position_folded(particle));
150 } else {
151 local_positions.emplace_back(traits.position(particle));
152 }
153 local_pids.emplace_back(traits.id(particle));
154 }
155
156 auto const argsort = detail::get_argsort(comm, local_pids, sorted_pids);
157
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);
161
162 if (comm.rank() != 0) {
163 return std::vector<pos_type>();
164 }
165
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);
171 }
172 }
173 assert(global_positions_flattened.size() == sorted_pids.size());
174
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]);
179 }
180
181 return positions_sorted;
182}
183} // namespace detail
184
185/**
186 * This class implements an interface to the `particle_observables` library that
187 * implements necessary algorithms needed for observables that are based on
188 * single particle properties.
189 * @tparam ObsType An observables composed of an algorithm from
190 * @ref src/particle_observables/include/particle_observables/algorithms.hpp
191 * and two particle properties.
192 *
193 * Example usage:
194 * @code{.cpp}
195 * using namespace ParticleObservables;
196 * using CenterOfMass = ParticleObservable<WeightedAverage<Position, Mass>>;
197 * @endcode
198 */
199template <class ObsType> class ParticleObservable : public PidObservable {
200public:
202 std::vector<std::size_t> shape() const override {
203 using std::declval;
204
205 return detail::shape_impl<decltype(declval<ObsType>()(
206 declval<ParticleReferenceRange const &>()))>::eval(ids().size());
207 }
208
209 template <typename T> struct is_map : std::false_type {};
210 template <typename T>
211 struct is_map<ParticleObservables::Map<T>> : std::true_type {};
212
213 std::vector<double>
214 evaluate(boost::mpi::communicator const &comm,
215 ParticleReferenceRange const &local_particles,
216 ParticleObservables::traits<Particle> const &) const override {
217 if constexpr (is_map<ObsType>::value) {
218 std::vector<double> local_traits{};
219 local_traits.reserve(local_particles.size());
220 Utils::flatten(ObsType{}(local_particles),
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));
226
227 std::vector<std::vector<double>> global_traits{};
228 boost::mpi::gather(comm, local_traits, global_traits, 0);
229
230 auto const argsort = detail::get_argsort(comm, local_pids, ids());
231
232 if (comm.rank() != 0) {
233 return {};
234 }
235
236 // get total size of the global traits vector
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(); });
240
241 auto const n_dims = size / ids().size();
242
243 std::vector<double> global_traits_flattened{};
244 global_traits_flattened.reserve(size);
245
246 for (auto const &vec : global_traits) {
247 for (auto const val : vec) {
248 global_traits_flattened.emplace_back(val);
249 }
250 }
251
252 std::vector<double> output{};
253 output.reserve(size);
254
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]);
258 }
259 }
260 return output;
261 } else {
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);
266
267 return result.first;
268 }
269 }
270};
271
272} // namespace Observables
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
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.
Definition flatten.hpp:64
std::vector< T, allocator< T > > vector
Definition vector.hpp:52