Add widget for MSD analysis#35
Conversation
Documentation build overview
14 files changed ·
|
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files
🚀 New features to boost your workflow:
|
There was a problem hiding this comment.
@PardhavMaradani, the PR looks like a good start for other lag-time based analyses. There are some issues with the current implementation that I have pointed out. Once resolved, it should be good to go.
Let me know if you need any help or clarification with any points raised or want to discuss possible solutions.
- Add a new SlidingWindowMSD class that is O(n) for each window
|
Hi @amruthesht , updated the code based on your example in #34 . I had to make a few minor changes to get this working correctly (verified that the plots were identical to the MDA's MSD ones):
Could you please review and let me know if there are any changes needed. Thanks |
|
Hi @amruthesht , Few additional thoughts to the previous comment: I tried the following changes shown in the 'Show code' section below: Show codeclass SlidingWindowMSD:
"""Sliding Window MSD
Calculate MSD for a sliding window of frames
"""
def __init__(self, u: mda.Universe, select: str = "all", msd_type: str = "xyz"):
self.u = u
self.select = select
self.msd_type = msd_type
self._parse_msd_type()
self.ag = u.select_atoms(self.select)
self.window_size = self.u.trajectory.buffer_size
# msd_values, w_sum, g_sum, g_count
self.msd_data = defaultdict(
lambda: [deque(maxlen=self.window_size), 0.0, 0.0, 0]
)
self.msd_data[0] = [deque([0]), 0.0, 0.0, 1]
def _parse_msd_type(self):
"""Sets up the desired dimensionality of the MSD."""
keys = {
"x": [0],
"y": [1],
"z": [2],
"xy": [0, 1],
"xz": [0, 2],
"yz": [1, 2],
"xyz": [0, 1, 2],
}
self._dim = keys[self.msd_type.lower()]
# pylint: disable=too-many-locals
def run(self, parallel: bool = False) -> tuple:
"""Run MSD for the current window"""
time_current = self.u.trajectory.ts.time
positions_current = self.ag.positions[:, self._dim]
for i in range(len(self.u.trajectory) - 1):
ts = self.u.trajectory[i] # set the buffered trajectory frame
delta_t = round(time_current - ts.time, 2)
disp = positions_current - self.ag.positions[:, self._dim]
msd = np.mean(np.sum(disp**2, axis=1))
msd_values, w_sum, g_sum, g_count = self.msd_data[delta_t]
if len(msd_values) == self.window_size:
w_sum -= msd_values[0]
msd_values.append(msd)
w_sum += msd
self.msd_data[delta_t] = [msd_values, w_sum, g_sum + msd, g_count + 1]
delta_t_values = sorted(self.msd_data.keys())
# v[0] = msd_values, v[1] = w_sum, v[2] = g_sum, v[3] = g_count
window_avgs, global_avgs = zip(
*[
(v[1] / len(v[0]), v[2] / v[3])
for dt in delta_t_values
for v in [self.msd_data[dt]]
]
)
return (
delta_t_values,
window_avgs,
global_avgs,
self.msd_data if parallel else None,
)With the above, we can plot both the global averages and the current window averages as shown in an example below: 
The mean calculation is also optimized like how you suggested in the comments in issue #34 code. Could you please share your thoughts on this. Thanks Update: I was able to get VACF for a sliding window (using this reference from MDAnalysis/transport-analysis) using the same approach as above - keep a running sum, delete the oldest when it gets removed from buffer, add the sum for the new entity, return the normalized avg. If this approach looks ok, I can create a PR for that after this one. Thanks |
|
@PardhavMaradani, thank you for responding promptly and taking into account the example code. I would let you choose the data structure that best suits how you envision the code, but the current implementation looks good to me. As for the averages, a windowed average doesn't make much sense here. The statistics for such a windowed average would be a little confusing. You might want to just stick with the updating global average as you have it. This way, the user can see a live-updating plot that updates the average based on new data as it is read by MDA. As for your point here
I was unsure how the final indexing in Overall, I would say the implementation looks good to me, great job! As an aside, since you have made this code contribution here, I would suggest opening up an issue on MDA, asking to support a FixedReferenceMSD apart from the |
- Remove windowed averages - Use pure numpy arrays instead of dicts and loops
|
Hi @amruthesht, thanks for reviewing the changes. I removed the windowed averages. I've added support for Here is an example with different combinations for large selections: 
Here is one for smaller selections: 
If the intention was not to display them, please let me know, I can remove it.
Sure. I'll make a note of this and open this along with other issues that come out of this project. Thanks Edit: I updated the plots. The previous ones seemed a bit off. I tried with both ( |
|
@PardhavMaradani, sorry I missed your comment and changes, but this looks good to me. As for the individual atom MSDs, you can keep a toggleable option to have them or not. The idea is to access those values in a data structure as the results of the MSD routine, which the user can use for other analyses if needed. As for the analysis, you are right that the user would be expected to provide or transform trajectories to remove any jumps. But that is something the user should be expected to do. But since you are providing MSD as a pre-defined widget, it makes sense to hardcode it as a part of the widget. The 'hook' in the MSD observed at long lag times is mostly an artifact of lower stats, which should ideally improve as the simulation runs longer and more data arrives (more averaging). |
2 participants
Changes made in this Pull Request:
PR Checklist
Hi @orbeckst, @jeremyleung521, @amruthesht, @HeydenLabASU,
This PR adds a new 'MSD Analysis' widget based on issue #34 raised by @amruthesht that uses
MDAnalysis.analysis.msd.EinsteinMSD.The widget runs every-frame and uses the timestep values from the sliding window buffer to create the plot. In the screenshots shown below, the 'Batch size' / 'Buffer size' is configured as 100 in 'Settings > Universe configuration' and hence the lag time shows up to 2 ps (each step is 0.02 ps in this example simulation).
Here is an example plot with all default values (
select='all', msd_type='xyz', fft=False, non_linear=False):Here is the same plot with a log scale (configurable from the GUI):
Here is the plot with
fft=True:Here is the plot with
non_linear=True:Here are the configurable inputs for this widget:
I added parallel support for this as it seemed to take a bit of time.
Here is a plot with a different selection phrase:
Thanks