Return correct shapes from Correlator

[jngrad]

Fixes #3577

Description of changes:

• reshape Correlator output arrays based on the dimensions of the input observables
• return only the correlated values from the Correlator observable
• get time lags and sample sizes from dedicated class methods

[review-notebook-app]

Check out this pull request on 

Review Jupyter notebook visual diffs & provide feedback on notebooks.

---

Powered by ReviewNB

[jngrad]

@KaiSzuttor while working on this PR, I realized the output of the Correlator class could be ambiguous. From the user perspective, correlation happens on observables with shapes, but internally we correlate flattened arrays, which can lead to different results. For example, if we correlate the positions of particles with square_distance_componentwise resp. scalar_product, we get shapes [X, N, 3] resp. [X, 1] (instead of [X, N, 1]). Below is a MWE on a modified version of this PR, where the two exceptions scalar_product requires vectors, but observable ... is a matrix in src/core/accumulators/Correlator.cpp were removed:

import espressomd
from espressomd.accumulators import Correlator
from espressomd.observables import ParticlePositions
import numpy as np

system = espressomd.System(box_l=[1.0, 1.0, 1.0])
system.time_step = 0.05
system.cell_system.skin = 0.4
system.part.add(pos=[(0, 0, 0), (0.5, .5, 0)])
system.thermostat.set_langevin(kT=1.37, gamma=2.4, seed=42)

pos_obs = ParticlePositions(ids=(0,1))

c_pos1 = Correlator(obs1=pos_obs, tau_lin=16, tau_max=20., delta_N=10,
                   corr_operation="square_distance_componentwise",
                   compress1="discard1")
c_pos2 = Correlator(obs1=pos_obs, tau_lin=16, tau_max=20., delta_N=10,
                   corr_operation="scalar_product", compress1="discard1")
system.auto_update_accumulators.add(c_pos1)
system.auto_update_accumulators.add(c_pos2)

system.integrator.run(10000)

c_pos1.finalize()
c_pos2.finalize()
print(c_pos1.result().shape)  # output: (33, 2, 3)
print(c_pos2.result().shape)  # output: (33, 1)

For the tensor_product resp. scalar_product I've added extra checks to prevent passing two matrices as input, as we don't actually do Kronecker products resp. per-axis scalar products. Do you think it's the correct approach? Do we have a use for implementing the Kronecker product and per-axis scalar product (which is a simpler form of fcs_acf)?

[KaiSzuttor]

so you mean the shape is a function not only of the observable type but also of the operation? Isn't that expected?

[jngrad]

sure, but the scalar_product currently ignores the observable shape:

espresso/src/core/accumulators/Correlator.cpp

Lines 234 to 235 in f22be4c

	} else if (corr_operation_name == "scalar_product") {
	m_dim_corr = 1;

this is surprising given that we already have the code logic to take the observable shape into account

[KaiSzuttor]

So, the correlator has a bug which needs a fix before we can go on with returning the correct shape, right?

[jngrad]

I just tried to implement a scalar_product_last_axis that would allow us to take e.g. a ParticlePositions(ids=(0,1)) and return the scalar product of the individual particles, e.g. a correlation shape of [N_tau, 2, 1], but it can't be done with our current infrastructure. The correlation functions take two flat arrays as input, so we cannot know whether we have two particles with 3 coordinates or 3 particles with 2 coordinates. The fcs_acf function "trusts" that the last dimension has size 3. It also checks the array dimension is a multiple of 3, but it's to avoid out-of-bounds errors, not to guarantee the user took the correct observable.
We also cannot provide information about the flat array dimensions as additional function arguments, because we would have to provide the same arguments to all other functions (they need to have the same signature because we store a function pointer).

This gives us 2 options: throw an error when doing a scalar product of matrices (that's currently implemented in acfa0b4) which is an API breaking change, or not throw an error and trust the user can figure out why the correlation shape is 1. Same thing for tensor product.

[KaiSzuttor]

so our correlation function behaves like numpy correlate and assumes 1d data... I think there is no generic solution to correlating multidimensional data because the algorithms are not implemented in a generic fashion.

[jngrad]

What should we do then? Right now I throw an error if the user accidentally passes a matrix instead of a vector to the scalar product, because it was not originally designed to handle this case. There was no mechanism to prevent it because it didn't really make sense back then to calculate the scalar product of a ParticlePositions(ids=(0,1)), as the user guide made it clear these positions were flattened. But now, the user guide shows these positions are stored in 2D matrices (but only at the script interface level!), which can give the impression to users that they should be able to "scalar product" these matrices.

I cannot think of an application where one would want to do a scalar product on flattened matrices, but maybe there is one, in which case this would no longer be possible because we do not offer the possibility to reshape an observable before passing it to the correlator.

We already have an m_correlation_args designed to help with this kind of situation, unfortunately it's constrained to a Vector3d:

• using a boost::variant instead of Vector3d (aa0b4f8) won't work because the script interface would need a visitor pattern
• storing the dimensionality of the data in e.g. the first index m_correlation_args would work, although it is now visible to the user and mutable (which is not good design)
• replacing the function pointer by a functor means the correlation functions have a state that has to be serialized
• adding extra parameters to all correlation functions doesn't scale; the core of the problem lies in the tight coupling between the Correlator class and its operators; right now they are implemented as free functions to be completely decoupled, which prevents us from generalizing the infrastructure

[KaiSzuttor]

i'm not so sure if this is a cat 1 issue to be resolved next... let's wait for @RudolfWeeber

[jngrad]

i'm not so sure if this is a cat 1 issue to be resolved next...

Agreed, refactoring the Correlator framework is out of scope for this PR. I've removed the API-breaking matrix assertions in 56764cb, we can always get back to that idea in the future. Let's move forward with this PR.

[KaiSzuttor]

https://numpy.org/doc/stable/reference/generated/numpy.c_.html

[KaiSzuttor]

https://numpy.org/doc/stable/reference/generated/numpy.c_.html

[KaiSzuttor]

LGTM