Refactor integrator framework to reduce coupling

[jngrad]

Fixes #3271

Description of changes:

• represent core integrators as individual Python classes to reduce coupling in the script interface
• use an IntegratorHandle Python class to store the currently active integrator
• make the steepest descent algorithm a fully functional integrator and remove a side effect in steepest_descent() that reverted the integ_switch value after integration
• move the integrator logic from espressomd.minimize_energy.MinimizeEnergy to espressomd.integrate.SteepestDescent
• keep espressomd.minimize_energy.MinimizeEnergy as a wrapper for backward compatibility
• API change:
  • no API change for integrators
  • API change for espressomd.minimize_energy.MinimizeEnergy: now a stateless free function that takes an espresso system as argument, renamed to steepest_descent

[review-notebook-app]

Check out this pull request on 

You'll be able to see Jupyter notebook diff and discuss changes. Powered by ReviewNB.

[codecov]

Codecov Report

❗ No coverage uploaded for pull request base (python@a5e10d8). Click here to learn what that means.
The diff coverage is n/a.

@@           Coverage Diff            @@
##             python   #3390   +/-   ##
========================================
  Coverage          ?     86%           
========================================
  Files             ?     538           
  Lines             ?   25299           
  Branches          ?       0           
========================================
  Hits              ?   21818           
  Misses            ?    3481           
  Partials          ?       0

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update a5e10d8...830c76b. Read the comment docs.

[jngrad]

Feel free to comment on this WIP. Possible talking points:

• better abstraction of the integrator framework and naming convention
• although MinimizeEnergy won't hinder further refactoring of the integrator framework anymore, it could be removed to reduce coupling in the script interface and make the API change non-silent

[fweik]

I'll have a read today and will give you some feedback.

[fweik]

f_max == 0 should be allowed.

[jngrad]

it is allowed

[fweik]

Generally I don't think the minimize_energy thing should be a member of the system class. It's not part of its state, so it is confusing to put it there. For the class itself, I think the init and disable logic is unneeded and error-prone. I think it would be a better design to just have it set the the integrator at the beginning of the minimize method, and restore it at the end. This would make the object stateless, so it could be a free function. As such it is useful to keep for the future, because minimize and switch back is a common operation.

[fweik]

If the goal is to keep the current syntax for the integrators (set_...) the IntegratorHandle construction you proposed is the best we can get. So if this is the goal I'm fine with this. I think in general the integrators should not be part of the system, because this blurs the line between the physics and the MD method, (integration is something that is done to the system, not an inherent property) but I understand that there is no interest in such fundamental change at the moment.

I think the logical next step from here would be to put the core integrators into classes, akin to the
python classes derived from class Integrator.

[RudolfWeeber]

Makes things much more clear on the Python side.
That's probably as far as we should take it, until we know what the implementation in the core will look like.

I think, the legacy Minimizeenergy interface can be removed. The 4.2 release will break the interface non-silently in a few places anyway.
Is there anything in the PR which should go into a bugfix release? If so, the old interface could be removed after cherry-picking.

On a different note, integrate_vv() is no longer responsible for Velocity-Verlet alone, so the _vv should probably be removed.

Is there a reason to manually allocate the memory for the steepest descent struct?

[fweik]

I think, the legacy Minimizeenergy interface can be removed.

As I was saying, this should be kept as a free function. That is a common think, and I can easily picture how forgetting to switch back the integrator after energy minimization can be a common error mode.

@RudolfWeeber I don't think that core refactoring is in scope for this PR, let's tackle this in the next step.

[jngrad]

I think it would be a better design to just have it set the the integrator at the beginning of the minimize method, and restore it at the end. This would make the object stateless, so it could be a free function. As such it is useful to keep for the future, because minimize and switch back is a common operation.

Although I initially planned to remove completely the MinimizeEnergy class and replace it by system.integrator.set_steepest_descent everywhere, I realized we might have a need for some Minimization feature that would take as argument a lambda function to decide when to stop the minimization. This logic is currently coded manually in a few places to integrate until inter-particle distances become large enough (ex: nacl.py#L92-L99). This could be achieved by having an extra method minimize(until=lambda: system.analysis.min_dist() < max_sigma) in the SteepestDescent(Integrator) class. The refactored MinimizeEnergy class in this WIP already has a high level of coupling with IntegratorHandle, and that's a strong sign that MinimizeEnergy should be part of the integrator abstraction.

This strategy doesn't seem satisfactory from an abstraction perspective: minimization algorithms like steepest descent, simulated annealing or genetic algorithms don't have a well-defined notion of time (i.e. at each "integration" step, the time step is different for every particle in the system) and can't really be considered as integrators. It's unfortunate the steepest descent algorithm was initially developed as a pseudo-integrator, but until we find a need to decouple the integration loop (e.g. to support other minimization algorithms), we would probably be better off by removing the MinimizeEnergy feature altogether.

Designing the minimize_energy feature as a free function means it will have to change the current integrator and restore it after minimization, which is more-or-less what the refactored MinimizeEnergy is now doing, but in a stateless form. This would be one of the few features in espresso where a function takes a System object as argument. Should we discuss this in an espresso meeting or in a dev meeting?

If the goal is to keep the current syntax for the integrators (set_...) the IntegratorHandle construction you proposed is the best we can get. So if this is the goal I'm fine with this.

Yes, that's the main goal. In the future, we could think of a new syntax like system.integrator = espressomd.Integrator.VelocityVerlet(), or introducing a Propagator class that encapsulates an Integrator instance and provides a list of compatible thermostats.

Is there anything in the PR which should go into a bugfix release?

Not really. The bug mentioned in the issue was already shipped in 4.1.1, and the side-effects mentioned in the commit messages are only relevant to this refactoring.

so the _vv should probably be removed.

Good idea.

Is there a reason to manually allocate the memory for the steepest descent struct?

Probably not, I'll have a look.

I don't think that core refactoring is in scope for this PR, let's tackle this in the next step.

The changes proposed by Rudolf seems small enough to me.

[fweik]

Designing the minimize_energy feature as a free function means it will have to change the current integrator and restore it after minimization, which is more-or-less what the refactored MinimizeEnergy is now doing, but in a stateless form. This would be one of the few features in espresso where a function takes a System object as argument.

There are a few, e.g. the polymer function. Arguable we should have more of those. I find the logically more coherent, and e.g. it makes this high level code much more easily testable e.g. with a mock system class. Personally, I also don't think that what we currently do and do not have in the python API is a a very good indicator for what is and is not a good idea.

[jngrad]

@fweik is 95549a9 what you have in mind?

[fweik]

@jngrad yes, this looks good to me. Any other opinions on this? With this the high level logic of minimize_energy can be tested in python...

[RudolfWeeber]

If we are keeping steepest descent as a free function (which I'd personally not do as per Pep8 reommendation to have peferrabley only one way to do a thing), the function should

• the name should be/contain steepest_descent rather than minimize energy, as that is, what's done
• the outcome (i.e., whether the target max force was reached) should be communicated, e.g., via return value.

[RudolfWeeber]

I agree with the goal to split system state from stuff that modifies the system. Probably, the core needs to be done first, though.

[fweik]

It seems to me that having two ways to do it for now is just the price of not changing the whole API at once. I don't have strong opinions on the name, and returning available information to the user is always a good idea.

[jngrad]

the outcome (i.e., whether the target max force was reached) should be communicated

This requires introducing a new global variable to store the state of the steepest descent, because steepest descent does not increment the simulation time and there is no simulation step counter.

[fweik]

@jngrad why can't you just return the number of steps that were executed from the function?

[jngrad]

@fweik this would require integrate() returning an int, steepest_descent() doing an MPI all_reduce, and writing a new mpi_call_all template that returns a value. I can have a look.

Also, we should make sure the return value of integrate() is only used for the steepest descent free function, and not in general for bookkeeping the total number of simulation steps, because steepest descent doesn't increment the simulation time. The canonical way of counting steps is:

espresso/src/core/io/writer/h5md_core.cpp

Line 385 in 467581d

step[0][0][0] = (int)std::round(sim_time / time_step);

[fweik]

You could always return the number of steps that are integrated. This also does not need communication, because this information is available on the head node, no?

[fweik]

LGTM except for one small issue.