add R_c (convergence fraction) convergence analysis

[xiki-tempula]

Fix #104

I have reimplemented the R_c implemented by @VOD555.
I have also implemented the A_c based on my understanding of the paper.

I have done a speed check, for a 32 window simulation where each window has 10,000 data points.

from alchemlyb.convergence import A_c
import time

data = pd.Series(data=range(10000))
data.attrs['temperature'] = 310
data.attrs['energy_unit'] = 'kcal/mol'
a=time.time()
value = A_c([data, ] * 32)
b=time.time()

Takes 0.01 second, so I think the performance is satisfying.

Includes commits from

• Fix dhdl2series and u_nk2series would not reattach the unit Fix dhdl2series and u_nk2series would not reattach the unit #248
• Change the plot_convergence so that it could take output from R_c Make the error optional for plot_convergence #249
• Remove the 'dhdl' method in decorrelating the u_nk Remove the 'dhdl' method in decorrelating the u_nk #250

[codecov]

Codecov Report

Merging #239 (2628740) into master (3809fdf) will increase coverage by 0.02%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master     #239      +/-   ##
==========================================
+ Coverage   98.65%   98.68%   +0.02%     
==========================================
  Files          26       26              
  Lines        1711     1750      +39     
  Branches      352      359       +7     
==========================================
+ Hits         1688     1727      +39     
  Misses          3        3              
  Partials       20       20              

Impacted Files	Coverage Δ
src/alchemlyb/preprocessing/subsampling.py	100.00% <ø> (ø)
src/alchemlyb/convergence/__init__.py	100.00% <100.00%> (ø)
src/alchemlyb/convergence/convergence.py	100.00% <100.00%> (ø)

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

[xiki-tempula]

@orbeckst Do you mind review this PR, please?
Note that the coverage is not as high as we want due to

for i in range(out_length):
    if all(conv[i:]):
        return i / length

This loop checks of the last i points in conv are true. conv is true when the backward and the forward are all with a certain distance to the end point. Given that the last forward point must be the same as the last backward point as they both cover the whole trajectory, conv[-1] is always true and thus, the else part of this loop is never covered, which lowers the coverage.

[orbeckst]

I might be able to look at it near the end of the week. Sorry, a bit busy.

@VOD555 , could you please review this PR? It's your method after all and you know it best. Although you cannot block or approve in alchemlyb, just state in your review what — if anything — needs to be done. I'll then use your review as a basis. Thank you!

[orbeckst]

Commenting on #239 (comment) : How about adding

else: # pragma: no cover
   # This loop checks of the last i points in conv are true. conv is true when the backward 
   # and the forward are all with a certain distance to the end point. Given that the last forward 
   # point must be the same as the last backward point as they both cover the whole trajectory, 
   # conv[-1] is always true and thus, the else part of this loop is never executed.
   raise AssertionError("Reached forbidden code path in convergence.R_c")

to make clear that the else part is never executed. Probably shorten the comment a bit... I just copied and pasted.

[VOD555]

It's better to also return g_forward and g_backward so that one can draw the convergence plot with these two series.

[xiki-tempula]

@VOD555 Thanks for the review. I have returned the forward and backward as dataframe now.

[orbeckst]

@xiki-tempula sorry, I'll look at this PR next week; I pretty much burned my "time to work on alchemlyb" on looking at all the other PRs.

@VOD555 , could you please have another look and review? Thank you!!

[orbeckst]

@xiki-tempula could you please resolve the conflicts with all the other recently merged PRs?

[orbeckst]

Overall this looks very good. I have some questions and would also like @VOD555 's input.

For consistency I'd keep input in $k_B T$ but if you have good reasons for absolute energies please say so.

[orbeckst]

Is this faster than abs?

[orbeckst]

use math mode for $A_c$ and equations like $A_c = 1$

[orbeckst]

math mode $R_c$

[orbeckst]

This function computes $R_c$ ...

[orbeckst]

Do you really mean R_c or should it be A_c.

I'd also put a SeeAlso to R_c and state that the R_c are calculated for the whole set using R_c().

[orbeckst]

in kT ?

[xiki-tempula]

@orbeckst Thanks very much for the review. Sorry for the delay. I have addressed all the issues.

[orbeckst]

Just minor doc things, otherwise looks good.

@VOD555 do you want to take a last look?

[xiki-tempula]

@orbeckst Thanks for the review. I have checked that the reference could be resolved.

[orbeckst]

Thanks for addressing all my nitpicks, much appreciated.

[orbeckst]

I made one change to the docs (wasn't in your edits) and I am now waiting for CI.

[orbeckst]

Thank you @xiki-tempula , a very welcome contribution!

[orbeckst]

I added @VOD555 as a "co-authored" by to the squash-merge because you said that you used his code as an initial basis and because he participated in the review. I put an extra note in the commit message that this version has your performance improvement included.