Use global counter in thermostats #3884

jngrad · 2020-09-03T20:00:53Z

Description of changes:

Remove RNG correlation stemming from seed offsets (fixes Thermostat noise for different seeds is correlated #3585)
- seeds are now used as keys
- a global, monotonically increasing counter is used in the thermostats
- the only way to reset this counter is to create a new System
Remove RNG correlation stemming from resetting sim_time or time_step during simulations with SD (fixes Noise behavior of stokesian dynamics and other thermostats differ #3840)
- the SD thermostat now uses the same RNG interface as other thermostats
Accelerate RNG unit tests (fixes Unit tests thermostats_test and random_test are very slow #3573)
- they now take 2 sec to run in coverage and sanitizer builds in CI
Split thermostats from integrators
- better separation of concerns
- globals from integrate.hpp are passed to thermostats as function parameters

In the previous implementation, thermostat noise for different seeds was correlated, the seed just shifted the sequence. For example with seed offset X, the cross-correlation used to have a peak at lag X. Co-authored-by: Jean-Noël Grad <[email protected]>

Co-authored-by: Jean-Noël Grad <[email protected]>

The SD external library uses Philox internally for both the CPU and GPU implementations. SD now uses the global counter instead of the simulation time divided by the time step, which are both mutable and can lead to correlated sequences when changed by the user.

Remove overhead from the std::function wrapper and heap allocation.

Reduce sample sizes and increase tolerances.

KaiSzuttor · 2020-09-04T11:48:30Z

wouldn't it be a good idea to have a global that monotonically increases with time. It may be useful also outside of the thermostat context at some point?

jngrad · 2020-09-04T12:28:17Z

wouldn't it be a good idea to have a global that monotonically increases with time. It may be useful also outside of the thermostat context at some point?

Yes, this is how the global counter is implemented. But we cannot rewind it by resetting the system time. It could be made a property of the integrator, or a property of the system to be able to re-use it elsewhere.

This counter is incremented by the integrator and is completely dissociated from the thermostat infrastructure.

Pass the current value of the integrator_counter to the force kernels as a parameter instead of relying on external linkage. There are two exceptions: thermalized bonds and constraints, which were not designed to be extensible in this way, and still rely on external linkage.

Header file thermostat.hpp no longer includes integrate.hpp.

Fixes clang-analyzer-deadcode.DeadStores warning.

jngrad · 2020-09-08T18:59:47Z

Using the same global counter in all thermostats solves the correlation issue but creates new problems. Some thermostats have specific needs: DPD needs extra incrementation when a DPD particle interacts with a constraint, and LB is not incremented when the integrator is steepest descent. Using a global counter incremented by the integrator means thermostat.cpp must include integrator.hpp and thermostat users such as src/core/constraints/ShapeBasedConstraint.cpp must include both integrator.hpp and thermostat.hpp. Having one counter per thermostat removes these issues.

jngrad · 2020-09-08T20:16:20Z

New attempt without the global counter in #3888.

KaiSzuttor · 2020-09-09T08:31:55Z

IIRC, the idea of the counter based RNG is that by knowing the value of the counter and two particle ids we determine one random number (in fact we could get 4). If we need more than one random number per particle pair and integration, the whole logic brakes down, right?

KaiSzuttor · 2020-09-09T08:33:12Z

LB is not incremented when the integrator is steepest descent

what do you mean by that?

KaiSzuttor · 2020-09-09T08:33:33Z

DPD needs extra incrementation when a DPD particle interacts with a constraint

why?

KaiSzuttor · 2020-09-09T08:36:10Z

Using a global counter incremented by the integrator means thermostat.cpp must include integrator.hpp

I think the counter state should be passed as an argument to all methods that need to know about it

KaiSzuttor · 2020-09-09T08:44:13Z

src/core/constraints/ShapeBasedConstraint.cpp

-        force1 +=
-            dpd_pair_force(p, part_rep, ia_params, dist_vec, dist, dist * dist);
-        // Additional use of DPD here requires counter increase
-        dpd_rng_counter_increment();


this counter has never been incremented twice... the control flow is in this if branch or in the one below

KaiSzuttor · 2020-09-09T08:47:21Z

src/core/integrate.hpp

@@ -144,3 +151,9 @@ int integrate_set_npt_isotropic(double ext_pressure, double piston,
                                bool zdir_rescale, bool cubic_box);

 #endif
+
+extern Utils::Counter<uint64_t> integrator_counter;


those declarations should come before the #endif of the header guard

jngrad · 2020-09-09T10:56:13Z

IIRC, the idea of the counter based RNG is that by knowing the value of the counter and two particle ids we determine one random number (in fact we could get 4). If we need more than one random number per particle pair and integration, the whole logic brakes down, right?

True for consumers of the thermostats. However, the thermostats functions can generate more than 1 random number per tuple (counter, seed, pid1, pid2) using the salt, because they are the one to call the RNG functions. This is used e.g. in NpT with NPTISO0_HALF_STEP1 and NPTISO0_HALF_STEP2. This cannot be used in DPD because the salt has to be known at compile time, while the DPD thermostat needs to generate an arbitrary number of random numbers per unique tuple, depending on the number of shape-based constraints in the system.

On a side-note, adding a new salt means shifting all the following salts in the RNGSalt enum, which is a silent change: scripts from one espresso minor version will generate different random numbers with the next espresso minor version where salts are shifted.

LB is not incremented when the integrator is steepest descent

what do you mean by that?

The function that increments the LB counters (CPU and GPU) is not called when using steepest descent:

espresso/src/core/integrate.cpp

Lines 267 to 268 in 2fbea6e

    
           if (integ_switch != INTEG_METHOD_STEEPEST_DESCENT) { 
        
             lb_lbfluid_propagate();

With a global counter, the LB RNG sequence will be incremented anyway, just like the other thermostats. Looking again at the code, this should not pose an issue.

DPD needs extra incrementation when a DPD particle interacts with a constraint

why?

Looking at the code, I would say it prevents correlation in the noise when you have more than one shape-based constraint. The shape-based constraint has a particle representation Particle ShapeBasedConstraint::part_rep which is default-constructed, i.e. its particle id is always -1, presumably to ensure real particles are numbered continuously.

This makes DPD quite different, because its counter gets out of sync with the number of integration steps, and gets incremented by system.integrator.run(0) while other thermostats don't. This is also why we have this logic at line 191 of the checkpoint test:

espresso/testsuite/python/test_checkpoint.py

Lines 187 to 191 in 2fbea6e

    
           def test_thermostat_DPD(self): 
        
               thmst = system.thermostat.get_state()[0] 
        
               self.assertEqual(thmst['type'], 'DPD') 
        
               self.assertEqual(thmst['kT'], 1.0) 
        
               self.assertEqual(thmst['seed'], 42 + 6)

Using a global counter incremented by the integrator means thermostat.cpp must include integrator.hpp

I think the counter state should be passed as an argument to all methods that need to know about it

I tried in 486612d but this means adding an extra argument to all classes in the constraint infrastructure just for DPD. Or adding a getter function to get the value of the global counter variable (to avoid exposing it via extern), in which case you need to include integrate.hpp.

KaiSzuttor · 2020-09-09T11:12:19Z

This cannot be used in DPD because the salt has to be known at compile time, while the DPD thermostat needs to generate an arbitrary number of random numbers per unique tuple, depending on the number of shape-based constraints in the system.

here the id of the particle representation of the constraint has to be used

KaiSzuttor · 2020-09-09T11:13:20Z

Looking at the code, I would say it prevents correlation in the noise when you have more than one shape-based constraint. The shape-based constraint has a particle representation Particle ShapeBasedConstraint::part_rep which is default-constructed, i.e. its particle id is always -1, presumably to ensure real particles are numbered continuously.

Then this needs to change.

KaiSzuttor · 2020-09-09T11:55:40Z

This cannot be used in DPD because the salt has to be known at compile time, while the DPD thermostat needs to generate an arbitrary number of random numbers per unique tuple, depending on the number of shape-based constraints in the system.

This would also not be the right approach. Actually, the interface for the random numbers should be such that it cannot be used the wrong way... Then we would see that we neither have the correct RNG interface for the NPT case nor have we solved the special case of constraint interaction via DPD (which should be handled just like any other particle particle interaction). The real issue here is that the particle can be default initialized which leads to the issue that you don't know which information of the particle representation in the constraint can actually be used and which not.

KaiSzuttor · 2020-09-09T11:58:02Z

With a global counter, the LB RNG sequence will be incremented anyway, just like the other thermostats. Looking again at the code, this should not pose an issue.

isn't the correct behavior to not propagate the system's state at all in case of steepest descent since there is no notion of time here? not sure about that

KaiSzuttor · 2020-09-09T11:59:42Z

I tried in 486612d but this means adding an extra argument to all classes in the constraint infrastructure just for DPD. Or adding a getter function to get the value of the global counter variable (to avoid exposing it via extern), in which case you need to include integrate.hpp.

What's wrong which the additional argument? Following your reasoning we should have more global variables.

KaiSzuttor · 2020-09-09T12:04:46Z

On a side-note, adding a new salt means shifting all the following salts in the RNGSalt enum, which is a silent change: scripts from one espresso minor version will generate different random numbers with the next espresso minor version where salts are shifted.

Not an issue if you add the salt at the end?

KaiSzuttor · 2020-09-09T14:38:32Z

offline discussion: the issue with the DPD interaction on ShapeBasedConstraint has it's origin in the fact that the particle representation in the Constraint is not a real particle and thus has not a meaningful particle id (which is in fact a part of the counter based RNG we are using). One possible solution could be to actually create the particle and only store the pointer to the particle in the constraint. The interaction of this solution with other parts of the code, however, is not clear.

KaiSzuttor · 2020-09-09T15:24:38Z

src/core/random.hpp

 *
 * @return Vector of uniform random numbers.
 */
 template <RNGSalt salt, size_t N = 3,
          std::enable_if_t<(N > 1) and (N <= 4), int> = 0>
-auto noise_uniform(uint64_t counter, int key1, int key2 = 0) {
+auto noise_uniform(uint64_t counter, uint32_t seed, int key1, int key2 = 0) {


i think what is missing is an abstraction for retrieving the random numbers specifically for the case of particle ids

jngrad · 2020-09-14T13:35:29Z

The DPD interaction currently prevents the use of a single counter, we'll need to find a long-term solution based on e.g. uuids (see #3892). Until then, the correlation bugs can be fixed by #3888, where thermostats still have their own counter.

fweik and others added 7 commits September 3, 2020 19:37

utils: Fixed parameter type for uint64_t u32_to_u64

9f9fc26

core: Use a single global counter for RNGs

cc97747

Co-authored-by: Jean-Noël Grad <[email protected]>

docs: Document unified thermostat infrastructure

7374dce

core: Make statistical unit tests faster

cc60ab9

Remove overhead from the std::function wrapper and heap allocation.

core: Make statistical unit tests faster

ea3f600

Reduce sample sizes and increase tolerances.

jngrad added Testcase Core Improvement BugFix labels Sep 3, 2020

jngrad added this to the Espresso 4.2 milestone Sep 3, 2020

tests: Increase RNG checks tolerance

54f8c74

jngrad added 9 commits September 4, 2020 20:35

tests: Increase RNG checks tolerance

8feab33

core: Move global counter to the integrator

2ca4b0e

This counter is incremented by the integrator and is completely dissociated from the thermostat infrastructure.

core: Break integrator/thermostat cyclic dependency

4e415bb

Header file thermostat.hpp no longer includes integrate.hpp.

core: Use vector operations

4be1ab1

core: Move NpT functions to dedicated file

5f785c5

core: Split thermostats from integrators

1f60e87

core: Remove temporary variable

48b7622

Fixes clang-analyzer-deadcode.DeadStores warning.

core: Header cleanup

50d097a

KaiSzuttor reviewed Sep 9, 2020

View reviewed changes

RudolfWeeber removed this from the Espresso 4.2 milestone Oct 26, 2020

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Use global counter in thermostats #3884

Use global counter in thermostats #3884

jngrad commented Sep 3, 2020 •

edited

Loading

KaiSzuttor commented Sep 4, 2020

jngrad commented Sep 4, 2020

jngrad commented Sep 8, 2020

jngrad commented Sep 8, 2020

KaiSzuttor commented Sep 9, 2020 •

edited

Loading

KaiSzuttor commented Sep 9, 2020

KaiSzuttor commented Sep 9, 2020

KaiSzuttor commented Sep 9, 2020

KaiSzuttor Sep 9, 2020

KaiSzuttor Sep 9, 2020

jngrad commented Sep 9, 2020

KaiSzuttor commented Sep 9, 2020

KaiSzuttor commented Sep 9, 2020

KaiSzuttor commented Sep 9, 2020

KaiSzuttor commented Sep 9, 2020

KaiSzuttor commented Sep 9, 2020

KaiSzuttor commented Sep 9, 2020

KaiSzuttor commented Sep 9, 2020

KaiSzuttor Sep 9, 2020

jngrad commented Sep 14, 2020

Use global counter in thermostats #3884

Are you sure you want to change the base?

Use global counter in thermostats #3884

Conversation

jngrad commented Sep 3, 2020 • edited Loading

KaiSzuttor commented Sep 4, 2020

jngrad commented Sep 4, 2020

jngrad commented Sep 8, 2020

jngrad commented Sep 8, 2020

KaiSzuttor commented Sep 9, 2020 • edited Loading

KaiSzuttor commented Sep 9, 2020

KaiSzuttor commented Sep 9, 2020

KaiSzuttor commented Sep 9, 2020

KaiSzuttor Sep 9, 2020

Choose a reason for hiding this comment

KaiSzuttor Sep 9, 2020

Choose a reason for hiding this comment

jngrad commented Sep 9, 2020

KaiSzuttor commented Sep 9, 2020

KaiSzuttor commented Sep 9, 2020

KaiSzuttor commented Sep 9, 2020

KaiSzuttor commented Sep 9, 2020

KaiSzuttor commented Sep 9, 2020

KaiSzuttor commented Sep 9, 2020

KaiSzuttor commented Sep 9, 2020

KaiSzuttor Sep 9, 2020

Choose a reason for hiding this comment

jngrad commented Sep 14, 2020

jngrad commented Sep 3, 2020 •

edited

Loading

KaiSzuttor commented Sep 9, 2020 •

edited

Loading