Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
i0 = trimmed_assignments[0][0]
if i0 == 1:
for m in trimmed_assignments:
m *= -1
m += 1
pairs = msm.draw_samples(trimmed_assignments, 2000)
samples = map_drawn_samples(pairs, data)
mu = np.mean(samples, axis=1)
eq(mu, np.array([[0., 0., 0.0], [25., 25., 25.]]), decimal=1)
# We should make sure we can sample from Trajectory objects too...
# Create a fake topology with 1 atom to match our input dataset
top = md.Topology.from_dataframe(pd.DataFrame({"serial":[0], "name":["HN"], "element":["H"], "resSeq":[1], "resName":"RES", "chainID":[0]}), bonds=np.zeros(shape=(0, 2), dtype='int'))
trajectories = [md.Trajectory(x[:, np.newaxis], top) for x in data] # np.newaxis reshapes the data to have a 40000 frames, 1 atom, 3 xyz
trj_samples = map_drawn_samples(pairs, trajectories)
mu = np.array([t.xyz.mean(0)[0] for t in trj_samples])
eq(mu, np.array([[0., 0., 0.0], [25., 25., 25.]]), decimal=1)
def test_superpose_refinds():
# make one frame far from the origin
normal = np.random.randn(1, 100, 3)
normal_xyz = normal.copy()
flipped = np.zeros_like(normal)
flipped[:, :50, :] = normal[:, 50:, :]
flipped[:, 50:, :] = normal[:, :50, :]
flipped_xyz = flipped.copy()
normal = md.Trajectory(xyz=normal, topology=None)
flipped = md.Trajectory(xyz=flipped, topology=None)
normal.superpose(flipped, atom_indices=np.arange(0, 50), ref_atom_indices=np.arange(50, 100))
eq(normal.xyz, normal_xyz)
flipped.superpose(normal, atom_indices=np.arange(50, 100), ref_atom_indices=np.arange(0, 50))
eq(flipped.xyz, flipped_xyz)
normal.superpose(flipped)
assert not np.allclose(normal.xyz, normal_xyz)
def test_closest_contact():
box_size = np.array([3.0, 4.0, 5.0])
traj = md.Trajectory(xyz=xyz * box_size, topology=None)
_verify_closest_contact(traj)
traj.unitcell_lengths = np.array([box_size for i in range(N_FRAMES)])
traj.unitcell_angles = np.array([[90.0, 90.0, 90.0] for i in range(N_FRAMES)])
_verify_closest_contact(traj)
traj.unitcell_angles = np.array([[80.0, 90.0, 100.0] for i in range(N_FRAMES)])
_verify_closest_contact(traj)
* What the dictionary actually contains
* ``paths_dict[idxs][type_of_path]["proj"]`` : ndarray of shape (n_points, proj_dim) with the
coordinates of the projection along the path
* paths_dict[idxs][type_of_path]["proj"] : geometries along the path
idata :
list of ndarrays with the the data in :obj:`projected_trajectories`
"""
if not isinstance(MD_trajectories, list):
MD_trajectories = [MD_trajectories]
if isinstance(MD_trajectories[0], _md.Trajectory):
src = MD_trajectories
else:
src = _source(MD_trajectories, top=MD_top)
idata = _data_from_input(projected_trajectories)
assert len(MD_trajectories) == len(idata), "Mismatch between the number of " \
"MD_trajectories and projected_trajectories: %u vs %u"%(len(MD_trajectories), len(idata))
#TODO: assert total_n_frames (strided) coincies with the n_frames in data
# What's the hightest dimensionlatiy that the input data allows?
input_dim = idata[0].shape[1]
if proj_idxs is None:
proj_idxs = _np.arange(n_projs)
else:
if isinstance(proj_idxs, int):
def _maketraj(self, positions):
Newxyz = np.zeros((1, self.ref.n_atoms, 3))
for i in range(len(positions)):
Newxyz[0,i,:] = ([positions[i]._value[0], positions[i]._value[1],
positions[i]._value[2]])
return mdj.Trajectory(Newxyz, self.ref.topology)
atom_indices = _np.arange(now.n_atoms)
else:
atom_indices = selection
# For the list of candidates, extract the closest one
history = now
for ii, cands in enumerate(path_of_candidates):
closest_to_now = _np.argmin(_md.rmsd(cands, now, atom_indices=atom_indices))
path_out.append(closest_to_now)
#print("choose frame %u from %u cands"%(path_out[-1], len(cands)))
now = cands[closest_to_now]
history = history.join(now)
if history_aware:
history.superpose(history, atom_indices=atom_indices)
xyz = history.xyz.mean(0)
now = _md.Trajectory(xyz, history.top)
return path_out
states = np.unique(np.hstack(([np.unique(disctraj) for disctraj in disctrajs])))
states = np.setdiff1d(states, [-1]) # exclude invalid states
cluster = [None] * (np.max(states) + 1)
for disctraj, traj, i in zip(disctrajs, trajs, xrange(len(trajs))):
assert len(disctraj) == traj.n_frames, 'Length of disctraj[%d] doesn\'t match number of frames in traj[%d].' % (
i, i)
for s in states:
match = (disctraj == s)
if np.count_nonzero(match) > 0:
if cluster[s] is None:
cluster[s] = traj.xyz[match, :, :]
else:
cluster[s] = np.concatenate((cluster[s], traj.xyz[match, :, :]), axis=0)
for i in xrange(len(cluster)):
if not cluster[i] is None:
cluster[i] = md.Trajectory(cluster[i], trajs[0].topology)
return cluster
def compute_neighbor_list(coords, neighbor_cutoff, max_num_neighbors,
periodic_box_size):
"""Computes a neighbor list from atom coordinates."""
N = coords.shape[0]
import mdtraj
traj = mdtraj.Trajectory(coords.reshape((1, N, 3)), None)
box_size = None
if periodic_box_size is not None:
box_size = np.array(periodic_box_size)
traj.unitcell_vectors = np.array(
[[[box_size[0], 0, 0], [0, box_size[1], 0], [0, 0, box_size[2]]]],
dtype=np.float32)
neighbors = mdtraj.geometry.compute_neighborlist(traj, neighbor_cutoff)
neighbor_list = {}
for i in range(N):
if max_num_neighbors is not None and len(neighbors[i]) > max_num_neighbors:
delta = coords[i] - coords.take(neighbors[i], axis=0)
if box_size is not None:
delta -= np.round(delta / box_size) * box_size
dist = np.linalg.norm(delta, axis=1)
sorted_neighbors = list(zip(dist, neighbors[i]))
sorted_neighbors.sort()
Parameters
----------
walker : object implementing the Walker interface
Returns
-------
min_distance : float
"""
cell_lengths, cell_angles = box_vectors_to_lengths_angles(walker.state['box_vectors'])
t2 = time.time()
# make a traj out of it so we can calculate distances through
# the periodic boundary conditions
walker_traj = mdj.Trajectory(walker.state['positions'],
topology=self._mdj_top,
unitcell_lengths=cell_lengths,
unitcell_angles=cell_angles)
t3 = time.time()
# calculate the distances through periodic boundary conditions
# and get hte minimum distance
min_distance = np.min(mdj.compute_distances(walker_traj,
it.product(self.ligand_idxs,
self.receptor_idxs),
periodic=self._periodic)
)
t4 = time.time()
logging.info("Make a traj: {0}; Calc dists: {1}".format(t3-t2,t4-t3))
return min_distance