Switch refactor for210

[Iximiel]

Description

I refactored the switching function in order to make it easier to edit/add new functions by transforming the choices in a strategies.

By doing so I removed the if else chains in SwitchingFunction::calculate() and in similar methods.

I changed basic/rt20 to a series of tests for all the various switching functions.

These modification brings some performance changes that impact for example COORDINATION in the following way:

The time data is the mean value of 200 runs on at 100 frame trajectory of 108 atoms, and is taken from the plumed timing report. The plotted data is the value of the row 4 Calculating (forward loop)). The red data points are the time of the new implementation over the time of the original, and they are referred to the right axis The vertical lines for the points are the standard deviation

The benches were run with the following commands:

key	COORDINATION command
cosinus	COORDINATION @GROUPS@ SWITCH={COSINUS R_0=2.6 }
cubic	COORDINATION @GROUPS@ SWITCH={CUBIC D_0=0 D_MAX=2.6 }
exp	COORDINATION @GROUPS@ SWITCH={EXP R_0=0.8 D_0=0.5 D_MAX=2.6 }
fastrational-MMeq2NN	COORDINATION @GROUPS@ SWITCH={RATIONAL R_0=1.3 NN=8 D_MAX=2.6 }
fastrational	COORDINATION @GROUPS@ SWITCH={RATIONAL R_0=1.3 NN=6 MM=10 D_MAX=2.6 }
fastrationalTemplated-MMeq2NN	COORDINATION @GROUPS@ SWITCH={RATIONAL R_0=1.3 NN=6 D_MAX=2.6 }
gaussian	COORDINATION @GROUPS@ SWITCH={GAUSSIAN R_0=1.0 D_0=0.3 D_MAX=2.6 }
lepton	COORDINATION @GROUPS@ SWITCH={MATHEVAL R_0=1.3 D_MAX=2.6 FUNC=1/(1+x^6) }
q	COORDINATION @GROUPS@ SWITCH={Q R_0=1.0 D_0=0.3 BETA=1.0 LAMBDA=1.0 REF=1.3 D_MAX=2.6 }
rational-MMeq2NN	COORDINATION @GROUPS@ SWITCH={RATIONAL R_0=1.3 NN=5 D_MAX=2.6 }
rational	COORDINATION @GROUPS@ SWITCH={RATIONAL R_0=1.3 NN=5 MM=11 D_MAX=2.6 }
smap	COORDINATION @GROUPS@ SWITCH={SMAP R_0=1.3 A=3 B=2 D_MAX=2.6 }
tanh	COORDINATION @GROUPS@ SWITCH={TANH R_0=1.3 D_MAX=2.6 }

For the _ns variants the NOSTRETCH keyword have been used.

Some clarification on the names for the rational function benchmarks:

• fastRational is when both NN and MM are even and so there is no need to use the square root of the distance
• MMeq2NN stand for MM=2*NN and gives room to an optimization of the rational
• fastrationalTemplated-MMeq2NN for some exponent combinations there is a compile-time optimized implementation (NN=6,MM=12 here)

There are still some corners that need to be smoothed:

• updating the documentation to make the user aware of the pre-compiled rational
• adding more options for the pre-compiled rationals, currently there is only 6/12 and 12/24
• maybe add a SwitchingFunction::set(std::unique_ptr<baseSwitch>&&) for declaring a new switch within any CV

Target release

I would like my code to appear in release 2.10

Type of contribution

• [X] changes to code or doc authored by PLUMED developers, or additions of code in the core or within the default modules
• [ ] changes to a module not authored by you
• [ ] new module contribution or edit of a module authored by you

Copyright

• [ ] I agree to transfer the copyright of the code I have written to the PLUMED developers or to the author of the code I am modifying.

• [ ] the module I added or modified contains a COPYRIGHT file with the correct license information. Code should be released under an open source license. I also used the command cd src && ./header.sh mymodulename in order to make sure the headers of the module are correct.

Tests

• [x] I added a new regtest or modified an existing regtest to validate my changes.
• [ ] I verified that all regtests are passed successfully on GitHub Actions.

[GiovanniBussi]

@Iximiel this is very useful I think 👏

A few comments:

• I noticed you added many regtest, are they passing also using the old implementation? I mean: you were aiming at identical results right?
• I noticed some comments in the Q variable. I think its derivatives were fixed recently. I guess you used the latest (fixed) version.
• I really do not like the vector of Lepton functions (that I implemented indeed :-)). Do you think there could be a nicer solution? Otherwise we can leave it as is
• stretch and nostretch have identical performances, as expected. If you wanted to optimize this as well I think you should have used a template variable. However, I don't think it's worth: using the nostretch version is discouraged and I don't think we need to optimize it.

updating the documentation to make the user aware of the pre-compiled rational

I think adding a sentence with "most commonly used rational functions are better optimized and might run faster" would be sufficient. People would not tune these coefficients for computational efficiency, but rather to have proper tail behaviors.

adding more options for the pre-compiled rationals, currently there is only 6/12 and 12/24

I've seen using other small numbers, such as n=4. I would just do something like N=2,4,6,8,10,12 (with M=2N), it's just a few more copy and paste if I understood properly the code

maybe add a SwitchingFunction::set(std::unique_ptr<baseSwitch>&&) for declaring a new switch within any CV

I am not sure I understand the syntax. Shall we perhaps have another register (as we have for CVs)? In the end I am not sure it's worth. We are talking about 1D functions, people can already play a lot with lepton, and then, if something becomes performance critical, one can directly modify the code. I think it's much less interesting wrt adding new CVs or biasing methods, for which we used registers to facilitate extensibility.

[Iximiel]

• I noticed you added many regtest, are they passing also using the old implementation? I mean: you were aiming at identical results right?

yes the regtest where make reseted with the old implementation, I left the date in a comment the config of the tests

• I noticed some comments in the Q variable. I think its derivatives were fixed recently. I guess you used the latest (fixed) version.

I ran the q test with the 2.10a to be sure about that

maybe add a SwitchingFunction::set(std::unique_ptr&&) for declaring a new switch within any CV

I am not sure I understand the syntax. Shall we perhaps have another register (as we have for CVs)? In the end I am not sure it's worth. We are talking about 1D functions, people can already play a lot with lepton, and then, if something becomes performance critical, one can directly modify the code. I think it's much less interesting wrt adding new CVs or biasing methods, for which we used registers to facilitate extensibility.

I don't want to propose a register for the switching function, the idea is that set would look like:

SwitchingFunction::set(std::unique_ptr<baseSwitch>&& ptr) {
function = ptr;
}

so that if one needs an ad hoc switch and does not want it depending on lepton in an eventual LOADed cv.cpp can write something that inherits from baseSwitch and put in that cv

• I really do not like the vector of Lepton functions (that I implemented indeed :-)). Do you think there could be a nicer solution? Otherwise we can leave it as is

I'll see what I can do about it

[Iximiel]

• I really do not like the vector of Lepton functions (that I implemented indeed :-)). Do you think there could be a nicer solution? Otherwise we can leave it as is

I refactored the thing by "compacting" the logic in a class As now is private in the switch implementation, if it is useful it may be moved to the lepton interface

[codecov-commenter]

Codecov Report

Attention: Patch coverage is 84.96732% with 69 lines in your changes are missing coverage. Please review.

Project coverage is 83.29%. Comparing base (4ceefe1) to head (ca0a0f7). Report is 1 commits behind head on master.

:exclamation: Current head ca0a0f7 differs from pull request most recent head 932fab0. Consider uploading reports for the commit 932fab0 to get more accurate results

Files	Patch %	Lines
src/cltools/SwitchingPlotter.cpp	25.75%	49 Missing :warning:
src/tools/SwitchingFunction.cpp	96.92%	11 Missing :warning:
src/adjmat/ContactMatrixShortcut.cpp	45.45%	6 Missing :warning:
src/tools/Grid.cpp	72.72%	3 Missing :warning:

:exclamation: Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1038      +/-   ##
==========================================
- Coverage   83.33%   83.29%   -0.04%     
==========================================
  Files         626      628       +2     
  Lines       59788    59987     +199     
==========================================
+ Hits        49822    49968     +146     
- Misses       9966    10019      +53     

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[GiovanniBussi]

@Iximiel if you are able to solve the problem with the numerical instabilities of the derivatives and it turns out to be a simple recipe such as:

• change the boundaries
• use a Taylor expansion in the middle Maybe we can backport it to v2.8?

Another independent thought. Perhaps something even better than a Taylor expansion (for our purposes) could be:

• choose the boundaries
• replace the function in the middle with a linear interpolation between the boundaries.

If the function can be computed in a numerically stable way at the boundaries, then it's sufficient to call it twice when the SwitchingFunction is created to find the expression of the straight line fitting those two points, and save the two coefficients (A*x+B). The value at r=R_0 will not be anymore exactly NN/MM, but this is not a problem. The derivative will be A over the entire range.

Even better: one could compute both the values and the derivatives at the boundary, and interpolate with a cubic spline. This will make also the first derivative continuous. If you have 4 conditions (two values and two derivatives) you can fit a polynomial function such as Ax^3+Bx^3+C*x+D, for which the derivative is analytical.

These approaches do not require to do the algebraic derivatives. An additional advantage is that they replace the function with another function that is very close to the original one (if we pick the boundaries well) but continuous. Continuity is actually a very important property for us (and it's the reason for the stretching thing)

[Iximiel]

@GiovanniBussi

@Iximiel if you are able to solve the problem with the numerical instabilities of the derivatives and it turns out to be a simple recipe such as:

• change the boundaries

• use a Taylor expansion in the middle Maybe we can backport it to v2.8?

I think it is possible, I touched not too many files and this could "just" be rebased on that. I just have to rework the C++17 shortcuts, like all the if constexpr that I used, and then I just have to prepare a small PR to revert the code to c++17 after the merge into master In the discussion of the error i the CI I posted some data using Taylor in the middle of the gap

Another independent thought. Perhaps something even better than a Taylor expansion (for our purposes) could be:

• choose the boundaries

• replace the function in the middle with a linear interpolation between the boundaries.

If the function can be computed in a numerically stable way at the boundaries, then it's sufficient to call it twice when the SwitchingFunction is created to find the expression of the straight line fitting those two points, and save the two coefficients (A*x+B). The value at r=R_0 will not be anymore exactly NN/MM, but this is not a problem. The derivative will be A over the entire range.

Even better: one could compute both the values and the derivatives at the boundary, and interpolate with a cubic spline. This will make also the first derivative continuous. If you have 4 conditions (two values and two derivatives) you can fit a polynomial function such as A_x^3+B_x^3+C*x+D, for which the derivative is analytical.

It should work, We just need to measure where the numerical instabilities stops, to put reasonable boundaries there

[GiovanniBussi]

Great! For backporting I just meant the derivative calculation, because it is a bug, not the full implementation... thanks!

[GiovanniBussi]

@Iximiel is this ready for merging?

[Iximiel]

@GiovanniBussi I updated the test to be correct to with the change of the rational, if it passes also the intel CI I think it is ready to push