Feature Request
Many scientifically relevant analyses require repeated access to data over a range of time points or simulation frames. The mathematical formulation behind most such analyses are time correlation-based computations (eg. autocorrelation). Typical examples of scientifically relevant quantities are MSD (mean-squared displacement), VACF (velocity autocorrelation) etc.
However, as described above, these analyses require repeated access to data over a range of time values. A MDA data structure like Universe would thus need access to data points at previous times—a limitation inherent to streaming.
Describe the solution you'd like
The solution would ideally require the Universe to store simulation frames from previous time points up to a certain user-defined number temporarily. This could be implemented via the BufferedTrajectory. However, unlike current batching strategy this would have to be implemented as FIFO Buffer with LIFO-esque consumption.
An example of the kind of analysis described above with rudimentary sample solution code, inspired by MDAnalysis' present mean-squared displacement routine.
import MDAnalysis as mda
u = mda.Universe("topol.tpr", "imd://localhost:8889")
if not isinstance(u.trajectory, BufferedTrajectory):
u.trajectory = BufferedTrajectory(u.trajectory, n_lags+1)
# n_lags is the number of lag times for the calculation of MSD
# could be done by making _maxlen initailization parameter in the class
msd_dict = collections.defaultdict(list)
msds_by_particle_dict = collections.defaultdict(list)
msd_dict[0] = [0]
msds_by_particle_dict[0.0] = [np.zeros(u.atoms.n_atoms)]
for ts in u.trajectory: # might wanna implemenet as seperate get_next_frame
# deque should already behave FIFO-esque with removing old frames when maxlen is reached
for i in range(1, len(u.trajectory)): # BufferedTrajectory needs a `len` method
delta_t = u.trajectory[0].time - u.trajectory[-i].time
# actual calculation - would be different for VACF
disp = u.trajectory[0].position - u.trajectory[-i].position
squared_disp = np.sum(disp**2, axis=1)
msd = np.mean(squared_disp)
msd_dict[delta_t].append(msd)
msds_by_particle_dict[delta_t].append(squared_disp)
delta_t_values = sorted(msd_dict.keys())
avg_msds = [np.mean(msd_dict[dt]) for dt in delta_t_values]
msds_by_particle_array = np.zeros((len(delta_t_values), u.atoms.n_atoms))
for idx, dt in enumerate(delta_t_values):
arr = np.vstack(msds_by_particle_dict[dt])
msds_by_particle_array[idx, :] = np.mean(arr, axis=0)
# this averaging can be further optimized by reusing means calculated at previous ts
# and calculating new means using msd values from newly arrived ts at relevant delta_t_values
# return avg_msds, delta_t_values, msds_by_particle_array
EDIT: Updated sample code to provided updating average results when a new frame is read into trajectory
Additional context
It might be useful to add buffered trajectory support with calculations modified in MDAnalysis to suit buffered and non-buffered universes in the future.
Feature Request
Many scientifically relevant analyses require repeated access to data over a range of time points or simulation frames. The mathematical formulation behind most such analyses are time correlation-based computations (eg. autocorrelation). Typical examples of scientifically relevant quantities are MSD (mean-squared displacement), VACF (velocity autocorrelation) etc.
However, as described above, these analyses require repeated access to data over a range of time values. A MDA data structure like
Universewould thus need access to data points at previous times—a limitation inherent to streaming.Describe the solution you'd like
The solution would ideally require the
Universeto store simulation frames from previous time points up to a certain user-defined number temporarily. This could be implemented via theBufferedTrajectory. However, unlike current batching strategy this would have to be implemented as FIFO Buffer with LIFO-esque consumption.An example of the kind of analysis described above with rudimentary sample solution code, inspired by
MDAnalysis' present mean-squared displacement routine.EDIT: Updated sample code to provided updating average results when a new frame is read into trajectory
Additional context
It might be useful to add buffered trajectory support with calculations modified in MDAnalysis to suit buffered and non-buffered universes in the future.