ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
Correlator.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-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#include "Correlator.hpp"
20
21#include <utils/Vector.hpp>
22#include <utils/math/sqr.hpp>
24
25#include <boost/archive/binary_iarchive.hpp>
26#include <boost/archive/binary_oarchive.hpp>
27#include <boost/iostreams/device/array.hpp>
28#include <boost/iostreams/stream.hpp>
29#include <boost/serialization/string.hpp>
30#include <boost/serialization/vector.hpp>
31
32#include <algorithm>
33#include <array>
34#include <cassert>
35#include <cmath>
36#include <cstddef>
37#include <functional>
38#include <numeric>
39#include <sstream>
40#include <stdexcept>
41#include <string>
42#include <vector>
43
44namespace {
45int min(int i, unsigned int j) { return std::min(i, static_cast<int>(j)); }
46} // namespace
47
48namespace Accumulators {
49/** Compress computing arithmetic mean: A_compressed=(A1+A2)/2 */
50std::vector<double> compress_linear(std::vector<double> const &A1,
51 std::vector<double> const &A2) {
52 assert(A1.size() == A2.size());
53 std::vector<double> A_compressed(A1.size());
54
55 std::transform(A1.begin(), A1.end(), A2.begin(), A_compressed.begin(),
56 [](double a, double b) -> double { return 0.5 * (a + b); });
57
58 return A_compressed;
59}
60
61/** Compress discarding the 1st argument and return the 2nd */
62std::vector<double> compress_discard1(std::vector<double> const &A1,
63 std::vector<double> const &A2) {
64 assert(A1.size() == A2.size());
65 std::vector<double> A_compressed(A2);
66 return A_compressed;
67}
68
69/** Compress discarding the 2nd argument and return the 1st */
70std::vector<double> compress_discard2(std::vector<double> const &A1,
71 std::vector<double> const &A2) {
72 assert(A1.size() == A2.size());
73 std::vector<double> A_compressed(A1);
74 return A_compressed;
75}
76
77std::vector<double> scalar_product(std::vector<double> const &A,
78 std::vector<double> const &B,
79 Utils::Vector3d const &) {
80 if (A.size() != B.size()) {
81 throw std::runtime_error(
82 "Error in scalar product: The vector sizes do not match");
83 }
84
85 auto const result = std::inner_product(A.begin(), A.end(), B.begin(), 0.0);
86 return {result};
87}
88
89std::vector<double> componentwise_product(std::vector<double> const &A,
90 std::vector<double> const &B,
91 Utils::Vector3d const &) {
92 std::vector<double> C(A.size());
93 if (A.size() != B.size()) {
94 throw std::runtime_error(
95 "Error in componentwise product: The vector sizes do not match");
96 }
97
98 std::transform(A.begin(), A.end(), B.begin(), C.begin(), std::multiplies<>());
99
100 return C;
101}
102
103std::vector<double> tensor_product(std::vector<double> const &A,
104 std::vector<double> const &B,
105 Utils::Vector3d const &) {
106 std::vector<double> C(A.size() * B.size());
107 auto C_it = C.begin();
108
109 for (double a : A) {
110 for (double b : B) {
111 *(C_it++) = a * b;
112 }
113 }
114
115 return C;
116}
117
118std::vector<double> square_distance_componentwise(std::vector<double> const &A,
119 std::vector<double> const &B,
120 Utils::Vector3d const &) {
121 if (A.size() != B.size()) {
122 throw std::runtime_error(
123 "Error in square distance componentwise: The vector sizes do not "
124 "match.");
125 }
126
127 std::vector<double> C(A.size());
128
129 std::transform(
130 A.begin(), A.end(), B.begin(), C.begin(),
131 [](double a, double b) -> double { return Utils::sqr(a - b); });
132
133 return C;
134}
135
136// note: the argument name wsquare denotes that its value is w^2 while the user
137// sets w
138std::vector<double> fcs_acf(std::vector<double> const &A,
139 std::vector<double> const &B,
140 Utils::Vector3d const &wsquare) {
141 if (A.size() != B.size()) {
142 throw std::runtime_error(
143 "Error in fcs_acf: The vector sizes do not match.");
144 }
145
146 auto const C_size = A.size() / 3;
147 assert(3 * C_size == A.size());
148
149 std::vector<double> C(C_size, 0);
150
151 for (std::size_t i = 0; i < C_size; i++) {
152 for (int j = 0; j < 3; j++) {
153 auto const &a = A[3 * i + j];
154 auto const &b = B[3 * i + j];
155
156 C[i] -= Utils::sqr(a - b) / wsquare[j];
157 }
158 }
159
160 std::transform(C.begin(), C.end(), C.begin(),
161 [](double c) -> double { return std::exp(c); });
162
163 return C;
164}
165
166void Correlator::initialize_operations() {
167 // Class members are assigned via the initializer list
168
169 if (m_tau_lin == 1) { // use the default
170 m_tau_lin = static_cast<int>(std::ceil(m_tau_max / m_dt));
171 m_tau_lin += m_tau_lin % 2;
172 }
173
174 if (m_tau_lin < 2) {
175 throw std::runtime_error("tau_lin must be >= 2");
176 }
177
178 if (m_tau_lin % 2) {
179 throw std::runtime_error("tau_lin must be divisible by 2");
180 }
181
182 if (m_tau_max <= m_dt) {
183 throw std::runtime_error("tau_max must be >= delta_t (delta_N too large)");
184 }
185 // set hierarchy depth which can accommodate at least m_tau_max
186 if ((m_tau_max / m_dt) < m_tau_lin) {
187 m_hierarchy_depth = 1;
188 } else {
189 auto const operand = (m_tau_max / m_dt) / double(m_tau_lin - 1);
190 assert(operand > 0.);
191 m_hierarchy_depth = static_cast<int>(std::ceil(1. + std::log2(operand)));
192 }
193
194 assert(A_obs);
195 assert(B_obs);
196 dim_A = A_obs->n_values();
197 dim_B = B_obs->n_values();
198
199 if (dim_A == 0u) {
200 throw std::runtime_error("dimension of first observable has to be >= 1");
201 }
202 if (dim_B == 0u) {
203 throw std::runtime_error("dimension of second observable has to be >= 1");
204 }
205
206 // choose the correlation operation
207 if (corr_operation_name == "componentwise_product") {
208 m_dim_corr = dim_A;
209 m_shape = A_obs->shape();
210 corr_operation = &componentwise_product;
211 m_correlation_args = Utils::Vector3d{0, 0, 0};
212 } else if (corr_operation_name == "tensor_product") {
213 m_dim_corr = dim_A * dim_B;
214 m_shape.clear();
215 m_shape.emplace_back(dim_A);
216 m_shape.emplace_back(dim_B);
217 corr_operation = &tensor_product;
218 m_correlation_args = Utils::Vector3d{0, 0, 0};
219 } else if (corr_operation_name == "square_distance_componentwise") {
220 m_dim_corr = dim_A;
221 m_shape = A_obs->shape();
222 corr_operation = &square_distance_componentwise;
223 m_correlation_args = Utils::Vector3d{0, 0, 0};
224 } else if (corr_operation_name == "fcs_acf") {
225 // note: user provides w=(wx,wy,wz) but we want to use
226 // wsquare=(wx^2,wy^2,wz^2)
227 if (not(m_correlation_args_input > Utils::Vector3d::broadcast(0.))) {
228 throw std::runtime_error("missing parameter for fcs_acf: w_x w_y w_z");
229 }
230 m_correlation_args = Utils::hadamard_product(m_correlation_args_input,
231 m_correlation_args_input);
232 if (dim_A % 3u)
233 throw std::runtime_error("dimA must be divisible by 3 for fcs_acf");
234 m_dim_corr = dim_A / 3u;
235 m_shape = A_obs->shape();
236 if (m_shape.back() != 3u)
237 throw std::runtime_error(
238 "the last dimension of dimA must be 3 for fcs_acf");
239 m_shape.pop_back();
240 corr_operation = &fcs_acf;
241 } else if (corr_operation_name == "scalar_product") {
242 m_dim_corr = 1u;
243 m_shape.clear();
244 m_shape.emplace_back(1u);
245 corr_operation = &scalar_product;
246 m_correlation_args = Utils::Vector3d{0, 0, 0};
247 } else {
248 throw std::invalid_argument("correlation operation '" +
249 corr_operation_name + "' not implemented");
250 }
251
252 // Choose the compression function
253 if (compressA_name == "discard2") {
254 compressA = &compress_discard2;
255 } else if (compressA_name == "discard1") {
256 compressA = &compress_discard1;
257 } else if (compressA_name == "linear") {
258 compressA = &compress_linear;
259 } else {
260 throw std::invalid_argument("unknown compression method '" +
261 compressA_name + "' for first observable");
262 }
263
264 if (compressB_name == "discard2") {
265 compressB = &compress_discard2;
266 } else if (compressB_name == "discard1") {
267 compressB = &compress_discard1;
268 } else if (compressB_name == "linear") {
269 compressB = &compress_linear;
270 } else {
271 throw std::invalid_argument("unknown compression method '" +
272 compressB_name + "' for second observable");
273 }
274}
275
276void Correlator::initialize_buffers() {
277 using index_type = decltype(result)::index;
278
279 A.resize(std::array<int, 2>{{m_hierarchy_depth, m_tau_lin + 1}});
280 std::fill_n(A.data(), A.num_elements(), std::vector<double>(dim_A, 0));
281 B.resize(std::array<int, 2>{{m_hierarchy_depth, m_tau_lin + 1}});
282 std::fill_n(B.data(), B.num_elements(), std::vector<double>(dim_B, 0));
283
284 n_data = 0;
285 A_accumulated_average = std::vector<double>(dim_A, 0);
286 B_accumulated_average = std::vector<double>(dim_B, 0);
287
288 auto const n_result = n_values();
289 n_sweeps = std::vector<std::size_t>(n_result, 0);
290 n_vals = std::vector<long>(m_hierarchy_depth, 0);
291
292 result.resize(std::array<std::size_t, 2>{{n_result, m_dim_corr}});
293 for (index_type i = 0; i < static_cast<index_type>(n_result); i++) {
294 for (index_type j = 0; j < static_cast<index_type>(m_dim_corr); j++) {
295 result[i][j] = 0.;
296 }
297 }
298
299 newest = std::vector<long>(m_hierarchy_depth, m_tau_lin);
300
301 tau.resize(n_result);
302 for (int i = 0; i < m_tau_lin + 1; i++) {
303 tau[i] = i;
304 }
305
306 for (int j = 1; j < m_hierarchy_depth; j++) {
307 for (int k = 0; k < m_tau_lin / 2; k++) {
308 tau[m_tau_lin + 1 + (j - 1) * m_tau_lin / 2 + k] =
309 (k + (m_tau_lin / 2) + 1) * (1 << j);
310 }
311 }
312}
313
314void Correlator::update(boost::mpi::communicator const &comm) {
315 if (finalized) {
316 throw std::runtime_error(
317 "No data can be added after finalize() was called.");
318 }
319
320 if (comm.rank() != 0) {
321 // worker nodes just need to update the observables and exit
322 A_obs->operator()(comm);
323 if (A_obs != B_obs) {
324 B_obs->operator()(comm);
325 }
326
327 return;
328 }
329
330 // We must now go through the hierarchy and make sure there is space for the
331 // new datapoint. For every hierarchy level we have to decide if it is
332 // necessary to move something
333 int highest_level_to_compress = -1;
334
335 t++;
336
337 // Let's find out how far we have to go back in the hierarchy to make space
338 // for the new value
339 {
340 auto const max_depth = m_hierarchy_depth - 1;
341 int i = 0;
342 while (true) {
343 if (i >= max_depth or n_vals[i] <= m_tau_lin) {
344 break;
345 }
346 auto const modulo = 1 << (i + 1);
347 auto const remainder = (t - (m_tau_lin + 1) * (modulo - 1) - 1) % modulo;
348 if (remainder != 0) {
349 break;
350 }
351 highest_level_to_compress += 1;
352 i++;
353 }
354 }
355
356 // Now we know we must make space on the levels 0..highest_level_to_compress
357 // Now let's compress the data level by level.
358
359 for (int i = highest_level_to_compress; i >= 0; i--) {
360 // We increase the index indicating the newest on level i+1 by one (plus
361 // folding)
362 newest[i + 1] = (newest[i + 1] + 1) % (m_tau_lin + 1);
363 n_vals[i + 1] += 1;
364 A[i + 1][newest[i + 1]] =
365 (*compressA)(A[i][(newest[i] + 1) % (m_tau_lin + 1)],
366 A[i][(newest[i] + 2) % (m_tau_lin + 1)]);
367 B[i + 1][newest[i + 1]] =
368 (*compressB)(B[i][(newest[i] + 1) % (m_tau_lin + 1)],
369 B[i][(newest[i] + 2) % (m_tau_lin + 1)]);
370 }
371
372 newest[0] = (newest[0] + 1) % (m_tau_lin + 1);
373 n_vals[0]++;
374
375 A[0][newest[0]] = A_obs->operator()(comm);
376 if (A_obs != B_obs) {
377 B[0][newest[0]] = B_obs->operator()(comm);
378 } else {
379 B[0][newest[0]] = A[0][newest[0]];
380 }
381
382 // Now we update the cumulated averages and variances of A and B
383 n_data++;
384 for (std::size_t k = 0; k < dim_A; k++) {
385 A_accumulated_average[k] += A[0][newest[0]][k];
386 }
387
388 for (std::size_t k = 0; k < dim_B; k++) {
389 B_accumulated_average[k] += B[0][newest[0]][k];
390 }
391
392 using index_type = decltype(result)::index;
393 // Now update the lowest level correlation estimates
394 for (long j = 0; j < min(m_tau_lin + 1, n_vals[0]); j++) {
395 auto const index_new = newest[0];
396 auto const index_old = (newest[0] - j + m_tau_lin + 1) % (m_tau_lin + 1);
397 auto const temp =
398 (corr_operation)(A[0][index_old], B[0][index_new], m_correlation_args);
399 assert(temp.size() == m_dim_corr);
400
401 n_sweeps[j]++;
402 for (index_type k = 0; k < static_cast<index_type>(m_dim_corr); k++) {
403 result[j][k] += temp[k];
404 }
405 }
406 // Now for the higher ones
407 for (int i = 1; i < highest_level_to_compress + 2; i++) {
408 for (long j = (m_tau_lin + 1) / 2 + 1; j < min(m_tau_lin + 1, n_vals[i]);
409 j++) {
410 auto const index_new = newest[i];
411 auto const index_old = (newest[i] - j + m_tau_lin + 1) % (m_tau_lin + 1);
412 auto const index_res =
413 m_tau_lin + (i - 1) * m_tau_lin / 2 + (j - m_tau_lin / 2 + 1) - 1;
414 auto const temp = (corr_operation)(A[i][index_old], B[i][index_new],
415 m_correlation_args);
416 assert(temp.size() == m_dim_corr);
417
418 n_sweeps[index_res]++;
419 for (index_type k = 0; k < static_cast<index_type>(m_dim_corr); k++) {
420 result[index_res][k] += temp[k];
421 }
422 }
423 }
424}
425
426int Correlator::finalize(boost::mpi::communicator const &comm) {
427 using index_type = decltype(result)::index;
428 if (finalized) {
429 throw std::runtime_error("Correlator::finalize() can only be called once.");
430 }
431 // We must now go through the hierarchy and make sure there is space for the
432 // new datapoint. For every hierarchy level we have to decide if it is
433 // necessary to move something
434
435 // mark the correlation as finalized
436 finalized = true;
437
438 // worker nodes don't need to do anything
439 if (comm.rank() != 0) {
440 return 0;
441 }
442
443 for (int ll = 0; ll < m_hierarchy_depth - 1; ll++) {
444 long vals_ll; // number of values remaining in the lowest level
445 if (n_vals[ll] > m_tau_lin + 1)
446 vals_ll = m_tau_lin + n_vals[ll] % 2;
447 else
448 vals_ll = n_vals[ll];
449
450 while (vals_ll) {
451 // Check, if we will want to push the value from the lowest level
452 auto highest_level_to_compress = (vals_ll % 2) ? ll : -1;
453
454 // Let's find out how far we have to go back in the hierarchy to make
455 // space for the new value
456 {
457 auto const max_depth = m_hierarchy_depth - 1;
458 int i = ll + 1; // lowest level for which to check for compression
459 while (highest_level_to_compress > -1) {
460 if (i >= max_depth or n_vals[i] % 2 == 0 or n_vals[i] <= m_tau_lin) {
461 break;
462 }
463 highest_level_to_compress += 1;
464 i++;
465 }
466 }
467 vals_ll -= 1;
468
469 // Now we know we must make space on the levels
470 // 0..highest_level_to_compress
471 // Now let's compress the data level by level.
472
473 for (int i = highest_level_to_compress; i >= ll; i--) {
474 // We increase the index indicating the newest on level i+1 by one (plus
475 // folding)
476 newest[i + 1] = (newest[i + 1] + 1) % (m_tau_lin + 1);
477 n_vals[i + 1] += 1;
478
479 (*compressA)(A[i][(newest[i] + 1) % (m_tau_lin + 1)],
480 A[i][(newest[i] + 2) % (m_tau_lin + 1)]);
481 (*compressB)(B[i][(newest[i] + 1) % (m_tau_lin + 1)],
482 B[i][(newest[i] + 2) % (m_tau_lin + 1)]);
483 }
484 newest[ll] = (newest[ll] + 1) % (m_tau_lin + 1);
485
486 // We only need to update correlation estimates for the higher levels
487 for (int i = ll + 1; i < highest_level_to_compress + 2; i++) {
488 for (long j = (m_tau_lin + 1) / 2 + 1;
489 j < min(m_tau_lin + 1, n_vals[i]); j++) {
490 auto const index_new = newest[i];
491 auto const index_old =
492 (newest[i] - j + m_tau_lin + 1) % (m_tau_lin + 1);
493 auto const index_res =
494 m_tau_lin + (i - 1) * m_tau_lin / 2 + (j - m_tau_lin / 2 + 1) - 1;
495
496 auto const temp = (corr_operation)(A[i][index_old], B[i][index_new],
497 m_correlation_args);
498 assert(temp.size() == m_dim_corr);
499
500 n_sweeps[index_res]++;
501 for (index_type k = 0; k < static_cast<index_type>(m_dim_corr); k++) {
502 result[index_res][k] += temp[k];
503 }
504 }
505 }
506 }
507 }
508 return 0;
509}
510
511std::vector<double> Correlator::get_correlation() {
512 using index_type = decltype(result)::index;
513 auto const n_result = n_values();
514 std::vector<double> res(n_result * m_dim_corr);
515
516 for (std::size_t i = 0; i < n_result; i++) {
517 auto const index = static_cast<index_type>(m_dim_corr * i);
518 for (index_type k = 0; k < static_cast<index_type>(m_dim_corr); k++) {
519 if (n_sweeps[i]) {
520 res[index + k] = result[static_cast<index_type>(i)][k] /
521 static_cast<double>(n_sweeps[i]);
522 }
523 }
524 }
525 return res;
526}
527
528std::vector<double> Correlator::get_lag_times() const {
529 std::vector<double> res(n_values());
530 std::ranges::transform(tau, res.begin(),
531 [dt = m_dt](auto const &a) { return a * dt; });
532 return res;
533}
534
536 std::stringstream ss;
537 boost::archive::binary_oarchive oa(ss);
538
539 oa << t;
540 oa << m_dt;
541 oa << m_shape;
542 oa << m_correlation_args_input;
543 oa << A;
544 oa << B;
545 oa << result;
546 oa << n_sweeps;
547 oa << n_vals;
548 oa << newest;
549 oa << A_accumulated_average;
550 oa << B_accumulated_average;
551 oa << n_data;
552
553 return ss.str();
554}
555
556void Correlator::set_internal_state(std::string const &state) {
557 namespace iostreams = boost::iostreams;
558 iostreams::array_source src(state.data(), state.size());
559 iostreams::stream<iostreams::array_source> ss(src);
560 boost::archive::binary_iarchive ia(ss);
561
562 ia >> t;
563 ia >> m_dt;
564 ia >> m_shape;
565 ia >> m_correlation_args_input;
566 ia >> A;
567 ia >> B;
568 ia >> result;
569 ia >> n_sweeps;
570 ia >> n_vals;
571 ia >> newest;
572 ia >> A_accumulated_average;
573 ia >> B_accumulated_average;
574 ia >> n_data;
575 initialize_operations();
576 m_system = nullptr;
577}
578
579} // namespace Accumulators
Vector implementation and trait types for boost qvm interoperability.
void const * m_system
for bookkeeping purposes
std::string get_internal_state() const final
void set_internal_state(std::string const &) final
std::vector< double > get_lag_times() const
int finalize(boost::mpi::communicator const &comm)
At the end of data collection, go through the whole hierarchy and correlate data left there.
std::vector< double > get_correlation()
Return correlation result.
void update(boost::mpi::communicator const &comm) override
The function to process a new datapoint of A and B.
static DEVICE_QUALIFIER constexpr Vector< T, N > broadcast(typename Base::value_type const &value)
Create a vector that has all entries set to the same value.
Definition Vector.hpp:110
std::vector< double > componentwise_product(std::vector< double > const &A, std::vector< double > const &B, Utils::Vector3d const &)
std::vector< double > tensor_product(std::vector< double > const &A, std::vector< double > const &B, Utils::Vector3d const &)
std::vector< double > compress_linear(std::vector< double > const &A1, std::vector< double > const &A2)
Compress computing arithmetic mean: A_compressed=(A1+A2)/2.
std::vector< double > scalar_product(std::vector< double > const &A, std::vector< double > const &B, Utils::Vector3d const &)
std::vector< double > compress_discard1(std::vector< double > const &A1, std::vector< double > const &A2)
Compress discarding the 1st argument and return the 2nd.
std::vector< double > compress_discard2(std::vector< double > const &A1, std::vector< double > const &A2)
Compress discarding the 2nd argument and return the 1st.
std::vector< double > fcs_acf(std::vector< double > const &A, std::vector< double > const &B, Utils::Vector3d const &wsquare)
std::vector< double > square_distance_componentwise(std::vector< double > const &A, std::vector< double > const &B, Utils::Vector3d const &)
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
Definition sqr.hpp:28
auto hadamard_product(Vector< T, N > const &a, Vector< U, N > const &b)
Definition Vector.hpp:379
int min(int i, unsigned int j)