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>
45int min(
int i,
unsigned int j) {
return std::min(i,
static_cast<int>(j)); }
51 std::vector<double>
const &A2) {
52 assert(A1.size() == A2.size());
53 std::vector<double> A_compressed(A1.size());
55 std::transform(A1.begin(), A1.end(), A2.begin(), A_compressed.begin(),
56 [](
double a,
double b) ->
double { return 0.5 * (a + b); });
63 std::vector<double>
const &A2) {
64 assert(A1.size() == A2.size());
65 std::vector<double> A_compressed(A2);
71 std::vector<double>
const &A2) {
72 assert(A1.size() == A2.size());
73 std::vector<double> A_compressed(A1);
78 std::vector<double>
const &B,
80 if (A.size() != B.size()) {
81 throw std::runtime_error(
82 "Error in scalar product: The vector sizes do not match");
85 auto const result = std::inner_product(A.begin(), A.end(), B.begin(), 0.0);
90 std::vector<double>
const &B,
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");
98 std::transform(A.begin(), A.end(), B.begin(), C.begin(), std::multiplies<>());
104 std::vector<double>
const &B,
106 std::vector<double> C(A.size() * B.size());
107 auto C_it = C.begin();
119 std::vector<double>
const &B,
121 if (A.size() != B.size()) {
122 throw std::runtime_error(
123 "Error in square distance componentwise: The vector sizes do not "
127 std::vector<double> C(A.size());
130 A.begin(), A.end(), B.begin(), C.begin(),
131 [](
double a,
double b) ->
double { return Utils::sqr(a - b); });
138std::vector<double>
fcs_acf(std::vector<double>
const &A,
139 std::vector<double>
const &B,
141 if (A.size() != B.size()) {
142 throw std::runtime_error(
143 "Error in fcs_acf: The vector sizes do not match.");
146 auto const C_size = A.size() / 3;
147 assert(3 * C_size == A.size());
149 std::vector<double> C(C_size, 0);
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];
160 std::transform(C.begin(), C.end(), C.begin(),
161 [](
double c) ->
double { return std::exp(c); });
166void Correlator::initialize_operations() {
169 if (m_tau_lin == 1) {
170 m_tau_lin =
static_cast<int>(std::ceil(m_tau_max / m_dt));
171 m_tau_lin += m_tau_lin % 2;
175 throw std::runtime_error(
"tau_lin must be >= 2");
179 throw std::runtime_error(
"tau_lin must be divisible by 2");
182 if (m_tau_max <= m_dt) {
183 throw std::runtime_error(
"tau_max must be >= delta_t (delta_N too large)");
186 if ((m_tau_max / m_dt) < m_tau_lin) {
187 m_hierarchy_depth = 1;
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)));
196 dim_A = A_obs->n_values();
197 dim_B = B_obs->n_values();
200 throw std::runtime_error(
"dimension of first observable has to be >= 1");
203 throw std::runtime_error(
"dimension of second observable has to be >= 1");
207 if (corr_operation_name ==
"componentwise_product") {
209 m_shape = A_obs->shape();
212 }
else if (corr_operation_name ==
"tensor_product") {
213 m_dim_corr = dim_A * dim_B;
215 m_shape.emplace_back(dim_A);
216 m_shape.emplace_back(dim_B);
219 }
else if (corr_operation_name ==
"square_distance_componentwise") {
221 m_shape = A_obs->shape();
224 }
else if (corr_operation_name ==
"fcs_acf") {
228 throw std::runtime_error(
"missing parameter for fcs_acf: w_x w_y w_z");
231 m_correlation_args_input);
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");
241 }
else if (corr_operation_name ==
"scalar_product") {
244 m_shape.emplace_back(1u);
248 throw std::invalid_argument(
"correlation operation '" +
249 corr_operation_name +
"' not implemented");
253 if (compressA_name ==
"discard2") {
255 }
else if (compressA_name ==
"discard1") {
257 }
else if (compressA_name ==
"linear") {
260 throw std::invalid_argument(
"unknown compression method '" +
261 compressA_name +
"' for first observable");
264 if (compressB_name ==
"discard2") {
266 }
else if (compressB_name ==
"discard1") {
268 }
else if (compressB_name ==
"linear") {
271 throw std::invalid_argument(
"unknown compression method '" +
272 compressB_name +
"' for second observable");
276void Correlator::initialize_buffers() {
277 using index_type =
decltype(result)::index;
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));
285 A_accumulated_average = std::vector<double>(dim_A, 0);
286 B_accumulated_average = std::vector<double>(dim_B, 0);
289 n_sweeps = std::vector<std::size_t>(n_result, 0);
290 n_vals = std::vector<long>(m_hierarchy_depth, 0);
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++) {
299 newest = std::vector<long>(m_hierarchy_depth, m_tau_lin);
301 tau.resize(n_result);
302 for (
int i = 0; i < m_tau_lin + 1; i++) {
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);
316 throw std::runtime_error(
317 "No data can be added after finalize() was called.");
320 if (comm.rank() != 0) {
322 A_obs->operator()(comm);
323 if (A_obs != B_obs) {
324 B_obs->operator()(comm);
333 int highest_level_to_compress = -1;
340 auto const max_depth = m_hierarchy_depth - 1;
343 if (i >= max_depth or n_vals[i] <= m_tau_lin) {
346 auto const modulo = 1 << (i + 1);
347 auto const remainder = (t - (m_tau_lin + 1) * (modulo - 1) - 1) % modulo;
348 if (remainder != 0) {
351 highest_level_to_compress += 1;
359 for (
int i = highest_level_to_compress; i >= 0; i--) {
362 newest[i + 1] = (newest[i + 1] + 1) % (m_tau_lin + 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)]);
372 newest[0] = (newest[0] + 1) % (m_tau_lin + 1);
375 A[0][newest[0]] = A_obs->operator()(comm);
376 if (A_obs != B_obs) {
377 B[0][newest[0]] = B_obs->operator()(comm);
379 B[0][newest[0]] = A[0][newest[0]];
384 for (std::size_t k = 0; k < dim_A; k++) {
385 A_accumulated_average[k] += A[0][newest[0]][k];
388 for (std::size_t k = 0; k < dim_B; k++) {
389 B_accumulated_average[k] += B[0][newest[0]][k];
392 using index_type =
decltype(result)::index;
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);
398 (corr_operation)(A[0][index_old], B[0][index_new], m_correlation_args);
399 assert(temp.size() == m_dim_corr);
402 for (index_type k = 0; k < static_cast<index_type>(m_dim_corr); k++) {
403 result[j][k] += temp[k];
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]);
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],
416 assert(temp.size() == m_dim_corr);
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];
427 using index_type =
decltype(result)::index;
429 throw std::runtime_error(
"Correlator::finalize() can only be called once.");
439 if (comm.rank() != 0) {
443 for (
int ll = 0; ll < m_hierarchy_depth - 1; ll++) {
445 if (n_vals[ll] > m_tau_lin + 1)
446 vals_ll = m_tau_lin + n_vals[ll] % 2;
448 vals_ll = n_vals[ll];
452 auto highest_level_to_compress = (vals_ll % 2) ? ll : -1;
457 auto const max_depth = m_hierarchy_depth - 1;
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) {
463 highest_level_to_compress += 1;
473 for (
int i = highest_level_to_compress; i >= ll; i--) {
476 newest[i + 1] = (newest[i + 1] + 1) % (m_tau_lin + 1);
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)]);
484 newest[ll] = (newest[ll] + 1) % (m_tau_lin + 1);
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;
496 auto const temp = (corr_operation)(A[i][index_old], B[i][index_new],
498 assert(temp.size() == m_dim_corr);
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];
512 using index_type =
decltype(result)::index;
514 std::vector<double> res(n_result * m_dim_corr);
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++) {
520 res[index + k] = result[
static_cast<index_type
>(i)][k] /
521 static_cast<double>(n_sweeps[i]);
529 std::vector<double> res(
n_values());
530 std::ranges::transform(tau, res.begin(),
531 [
dt = m_dt](
auto const &a) { return a * dt; });
536 std::stringstream ss;
537 boost::archive::binary_oarchive oa(ss);
542 oa << m_correlation_args_input;
549 oa << A_accumulated_average;
550 oa << B_accumulated_average;
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);
565 ia >> m_correlation_args_input;
572 ia >> A_accumulated_average;
573 ia >> B_accumulated_average;
575 initialize_operations();
Vector implementation and trait types for boost qvm interoperability.
void const * m_system
for bookkeeping purposes
std::size_t n_values() const
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.
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.
auto hadamard_product(Vector< T, N > const &a, Vector< U, N > const &b)
int min(int i, unsigned int j)