From 43af36b0010a64954c4343fa314750e6aa5c1632 Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Tue, 4 Sep 2018 14:06:57 -0500
Subject: [PATCH] Updated algorithm that fits splines into beads to extend
farther into intrahelically connected segments, and to use linear
interpolation for quaternions, resulting in better all-atom models
---
dnarbd/segmentmodel.py | 25 ++++++++++++++++---------
1 file changed, 16 insertions(+), 9 deletions(-)
diff --git a/dnarbd/segmentmodel.py b/dnarbd/segmentmodel.py
index cf82cd3..cf834d6 100644
--- a/dnarbd/segmentmodel.py
+++ b/dnarbd/segmentmodel.py
@@ -1410,8 +1410,9 @@ class SegmentModel(ArbdModel):
filter_fn = lambda x: x is not None and x not in beads
bs = list( filter( filter_fn, B.particle.intrahelical_neighbors ) )
beads.update(bs)
- bs = list( filter( filter_fn, [n for b in bs for n in b.intrahelical_neighbors] ) )
- beads.update(bs)
+ for i in range(3):
+ bs = list( filter( filter_fn, [n for b in bs for n in b.intrahelical_neighbors] ) )
+ beads.update(bs)
beads = list(beads)
@@ -1433,12 +1434,21 @@ class SegmentModel(ArbdModel):
beads.remove(b)
contours = [b.get_contour_position(s) for b in beads]
+ contours = np.array(contours, dtype=np.float16) # deliberately use low precision
contours,ids = np.unique(contours, return_index=True)
- if np.any( (contours[:-1] - contours[1:])**2 < 1e-6 ):
+ if np.any( (contours[:-1] - contours[1:])**2 < 1e-8 ):
pdb.set_trace()
## TODO: keep closest beads beyond +-1.5 if there are fewer than 2 beads
- tmp = [ci for ci in zip(contours,ids) if np.abs(ci[0]-0.5) < 1]
+ tmp = []
+ dist = 1
+ while len(tmp) < 5 and dist < 3:
+ tmp = [ci for ci in zip(contours,ids) if np.abs(ci[0]-0.5) < dist]
+ dist += 0.1
+
+ if len(tmp) <= 1:
+ raise Exception("Failed to fit spline into segment {}".format(s))
+
contours = [ci[0] for ci in tmp]
ids = [ci[1] for ci in tmp]
@@ -1514,11 +1524,8 @@ class SegmentModel(ArbdModel):
# pdb.set_trace()
quats = np.array(quats)
- ## TODOTODO test smoothing
- try:
- tck, u = interpolate.splprep( quats.T, u=contours, s=3, k=3 )
- except:
- tck, u = interpolate.splprep( quats.T, u=contours, s=0, k=1 )
+ # tck, u = interpolate.splprep( quats.T, u=contours, s=3, k=3 ) ;# cubic spline not as good
+ tck, u = interpolate.splprep( quats.T, u=contours, s=0, k=1 )
s.quaternion_spline_params = tck
def _generate_bead_model(self,
--
GitLab