zemanj
commented
5 years ago Trajectory Reading Performance Impact
I tested the runtimes of trajectory reading (from a RAM disk) for various trajectory formats.
tl;dr: As expected, DCD becomes significantly slower when holding positions in C-contiguous arrays. Reading any other format becomes either faster or the change is negligible.
Test code:
import os
import shutil
from collections import defaultdict
from warnings import filterwarnings
from timeit import timeit
import MDAnalysis as mda
from MDAnalysisTests.datafiles import *
filterwarnings('ignore')
###### HELPER FUNCTIONS ######
def copy_trajectories(destination, trajectories):
for i in range(len(trajectories)):
trjfile, name = trajectories[i]
if isinstance(trjfile, tuple):
trjfile = list(trjfile)
for j in range(len(trjfile)):
dest = os.path.join(destination, os.path.basename(trjfile[j]))
shutil.copyfile(trjfile[j], dest)
trjfile[j] = dest
trjfile = tuple(trjfile)
trajectories[i] = (trjfile, name)
else:
dest = os.path.join(destination, os.path.basename(trjfile))
shutil.copyfile(trjfile, dest)
trajectories[i] = (dest, name)
def universe(trjfile, name):
if isinstance(trjfile, tuple):
return mda.Universe(*trjfile, format=formats[name])
return mda.Universe(trjfile, format=formats[name])
def readall(trjfile, name):
trj = universe(trjfile, name).trajectory
for _ in trj:
pass
def get_atoms_frames(trjfile, name):
u = universe(trjfile, name)
return u.atoms.n_atoms, u.trajectory.n_frames
def time_trajectory_read(trjfile, name):
n_atoms, n_frames = get_atoms_frames(trjfile, name)
normfactor = 1.0 / (n_runs * n_frames * n_atoms)
# Read once:
readall(trjfile, name)
# Get runtime per frame:
t = timeit(lambda: readall(trjfile, name), number=n_runs) * normfactor
# Print result:
print("{:14s} {:.3e} {:>5} {:>6} {:>5}"
"".format(name, t, n_runs, n_frames, n_atoms))
###### SETUP ######
n_runs = 100
formats = defaultdict(lambda: None)
formats['LAMMPSDUMP'] = 'LAMMPSDUMP'
formats['DLP'] = 'HISTORY'
ramdisk = '/media/ramdisk'
trajectories = [(CRD, 'CRD'),
(DCD, 'DCD'),
(DLP_HISTORY, 'DLP'),
(DMS, 'DMS'),
(GMS_ASYMSURF, 'GMS'),
(GRO, 'GRO'),
(GSD, 'GSD'),
(LAMMPSdata, 'LAMMPS'),
(LAMMPSDUMP, 'LAMMPSDUMP'),
(MMTF, 'MMTF'),
(MMTF_gz, 'MMTF_gz'),
(mol2_molecules, 'MOL2'),
(NCDF, 'NCDF'),
(PDB_multiframe, 'PDB_multiframe'),
(PDBQT_input, 'PDBQT'),
(PDB_xvf, 'PDB_xvf'),
(PQR, 'PQR'),
(TRR, 'TRR'),
((TRZ_psf, TRZ), 'TRZ'),
(TXYZ, 'TXYZ'),
(XTC, 'XTC'),
(XYZ, 'XYZ')]
copy_trajectories(ramdisk, trajectories)
###### RUNTIME MEASUREMENT ######
print("Format s/frame/atom runs frames atoms")
print("--------------------------------------------------")
for trjfile, name in trajectories:
time_trajectory_read(trjfile, name)Comparison of Test Results: new (C-contiguous) vs. old (F-contiguous)
new, C-cont. old, F-cont. (old-new)/old
Format s/frame/atom s/frame/atom relative gain
------------------------------------------------------
CRD 1.060E-05 1.084E-05 2.21%
DCD 2.554E-08 2.171E-08 -17.64%
DLP 1.284E-05 1.311E-05 2.06%
DMS 1.725E-05 1.821E-05 5.27%
GMS 6.183E-05 6.282E-05 1.58%
GRO 7.737E-06 7.962E-06 2.83%
GSD 5.728E-07 5.732E-07 0.07%
LAMMPS 1.462E-05 1.477E-05 1.02%
LAMMPSDUMP 3.420E-05 3.415E-05 -0.15%
MMTF 7.418E-06 7.407E-06 -0.15%
MMTF_gz 4.304E-06 4.347E-06 0.99%
MOL2 3.149E-06 3.175E-06 0.82%
NCDF 2.225E-08 2.262E-08 1.64%
PDB_multiframe 2.276E-06 2.315E-06 1.68%
PDBQT 5.230E-06 5.253E-06 0.44%
PDB_xvf 1.015E-05 1.024E-05 0.88%
PQR 4.238E-06 4.245E-06 0.16%
TRR 1.537E-07 1.555E-07 1.16%
TRZ 7.711E-07 8.087E-07 4.65%
TXYZ 1.495E-04 1.485E-04 -0.67%
XTC 4.778E-08 4.868E-08 1.85%
XYZ 1.512E-06 1.518E-06 0.40%
------------------------------------------------------
Average 1.584E-05 1.593E-05 0.58%
orbeckst
kain88-de
tylerjereddy
richardjgowers
Is your feature request related to a problem? Please describe. Currently, all functions that operate (especially in-place) on atom positions obtain a copy of atom positions via
AtomGroup.positions, regardless of whether positions are modified in- or out-of-place, or they are merely input data that is not altered. In some (if not many) cases, a tremendous speed-up could be achieved if the respective functions could obtain a view instead of a copy of an AtomGroup's positions. At the moment, this is not possible for two reasons:AtomGroup.positionsgets its data fromUniverse.trajectory.ts.positionsvia fancy indexing, and fancy indexing means copying.Universe.trajectory.ts.positionsrefers tocoordinates.base.Timestep._pos, which is a Fortran-contiguous array. Therefore, potentially obtained views of that array are unsuitable for use with C-level functions.Describe the solution you'd like Change
coordinates.base.Timestep.orderfrom'F'to'C'and provide a way of obtaining a (if possible, C-contiguous) view of an AtomGroup's positions.Describe alternatives you've considered Invent computers with two-dimensional main memory...
Additional context Two candidates that come to my mind in terms of speed-up are
*Group.center*()and*Group.wrap()when the index array*Group.atoms.ixis a (contiguous) range.A possible drawback of changing the positions/velocities/forces memory layout is that reading and (probably to a lesser extent) writing of DCD trajectories will be a bit slower since DCD stores positions etc. as Fortran-contiguous arrays.
If you have any additional concerns, please let me know!