From 3a7e63fc671ffee571cbb42d424015e8aee77aa0 Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Sat, 15 Sep 2018 17:42:17 -0500
Subject: [PATCH] Fix spline interpolation; added routines for transforming
entire model
---
mrdna/coords.py | 2 +-
mrdna/segmentmodel.py | 41 ++++++++++++++++++++++++++++-------------
2 files changed, 29 insertions(+), 14 deletions(-)
diff --git a/mrdna/coords.py b/mrdna/coords.py
index b74906d..41628e8 100644
--- a/mrdna/coords.py
+++ b/mrdna/coords.py
@@ -150,7 +150,7 @@ def quaternion_from_matrix( R ):
def rotationAboutAxis(axis,angle, normalizeAxis=True):
if normalizeAxis: axis = axis / np.linalg.norm(axis)
- angle *= 0.5 * np.pi/180
+ angle = angle * 0.5 * np.pi/180
cos = np.cos( angle )
sin = np.sin( angle )
q = [cos] + [sin*x for x in axis]
diff --git a/mrdna/segmentmodel.py b/mrdna/segmentmodel.py
index 0682714..61ddfca 100644
--- a/mrdna/segmentmodel.py
+++ b/mrdna/segmentmodel.py
@@ -420,7 +420,7 @@ class Segment(ConnectableElement, Group):
ids = list(filter(lambda i: position_filter(r[i,:]), ids))
return ids
- def translate(self, translation_vector, about=None, position_filter=None, contour_filter=None):
+ def translate(self, translation_vector, position_filter=None, contour_filter=None):
dr = np.array(translation_vector)
tck, u = self.position_spline_params
r = self.contour_to_position(u)
@@ -429,7 +429,7 @@ class Segment(ConnectableElement, Group):
## Translate
r[ids,:] = r[ids,:] + dr[np.newaxis,:]
- self.set_position_splines(r,u)
+ self.set_splines(u,r)
def rotate(self, rotation_matrix, about=None, position_filter=None, contour_filter=None):
tck, u = self.position_spline_params
@@ -445,16 +445,16 @@ class Segment(ConnectableElement, Group):
## TODO: do this more efficiently
r[ids,:] = np.array([rotation_matrix.dot(r[i,:]-dr) + dr for i in range(r.shape[0])])
- self.set_position_splines(r,u)
+ self.set_splines(u,r)
if self.quaternion_spline_params is not None:
## TODO: performance: don't shift between quaternion and matrix representations so much
tck, u = self.quaternion_spline_params
- orientations = [self.contour_to_orientation(v)) for v in u]
+ orientations = [self.contour_to_orientation(v) for v in u]
for i in ids:
orientations[i,:] = rotation_matrix.dot(orientations[i])
quats = [quaternion_from_matrix(o) for o in orientations]
- self.set_orientation_splines(quats,u)
+ self.set_orientation_splines(u, quats)
def _set_splines_from_ends(self):
self.quaternion_spline_params = None
@@ -1386,6 +1386,7 @@ class SegmentModel(ArbdModel):
self._bonded_potential = dict() # cache for bonded potentials
self._generate_strands()
+ self.potentials = []
self._generate_bead_model( max_basepairs_per_bead, max_nucleotides_per_bead, local_twist, escapable_twist)
self.useNonbondedScheme( nbDnaScheme )
@@ -1540,6 +1541,26 @@ class SegmentModel(ArbdModel):
print("WARNING: called deprecated command '_update_segment_positions; use 'update_splines' instead")
return self.update_splines(bead_coordinates)
+ ## Operations on spline coordinates
+ def translate(self, translation_vector, position_filter=None):
+ for s in self.segments:
+ s.translate(translation_vector, position_filter=position_filter)
+
+ def rotate(self, rotation_matrix, about=None, position_filter=None):
+ for s in self.segments:
+ s.rotate(rotation_matrix, about=about, position_filter=position_filter)
+
+ def get_center(self, include_ssdna=False):
+ if include_ssdna:
+ segments = self.segments
+ else:
+ segments = list(filter(lambda s: isinstance(s,DoubleStrandedSegment),
+ self.segments))
+ centers = [s.get_center() for s in segments]
+ weights = [s.num_nt*2 if isinstance(s,DoubleStrandedSegment) else s.num_nt for s in segments]
+ # centers,weights = [np.array(a) for a in (centers,weights)]
+ return np.average( centers, axis=0, weights=weights)
+
def update_splines(self, bead_coordinates):
""" Set new function for each segments functions
contour_to_position and contour_to_orientation """
@@ -1626,14 +1647,8 @@ class SegmentModel(ArbdModel):
# print(" {}[{:03f}]: {:0.3f}".format(b.parent.name, b.contour_position, c) )
- ret = interpolate.splprep( positions, u=contours, s=0, k=1, full_output=1 )
- tck = ret[0][0]
+ tck, u = interpolate.splprep( positions, u=contours, s=0, k=1 )
- ret = interpolate.splprep( positions, u=contours, s=0, k=1, full_output=1 )
- tck2 = ret[0][0]
-
- assert( tck[0][0] == tck2[0][0] )
-
# if len(beads) < 8:
# ret = interpolate.splprep( positions, u=contours, s=0, k=1, full_output=1 )
# tck = ret[0][0]
@@ -1652,7 +1667,7 @@ class SegmentModel(ArbdModel):
# if ret[2] > 0:
# pdb.set_trace()
- s.position_spline_params = (tck,ret[1][0])
+ s.position_spline_params = (tck,u)
""" Get twist """
cb = [e for e in zip(contours,beads) if 'orientation_bead' in e[1].__dict__]
--
GitLab