From b670c7f05d2ad4a2dba71efc471a46ffaa23491f Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Wed, 11 Apr 2018 11:12:14 -0500
Subject: [PATCH] Conditionally update coordinates using linear interpolation
 when there are not many beads

---
 cadnano_segments.py | 20 ++++++++++----------
 run.py              | 10 ++++++----
 segmentmodel.py     | 46 +++++++++++++++++++++++++++++++++++++++------
 3 files changed, 56 insertions(+), 20 deletions(-)

diff --git a/cadnano_segments.py b/cadnano_segments.py
index 32c2104..626ceea 100644
--- a/cadnano_segments.py
+++ b/cadnano_segments.py
@@ -526,15 +526,15 @@ def run_simulation_protocol( output_name, job_id, json_data,
                              sequence=None,
                              remove_long_bonds=False,
                              gpu = 0,
-                             directory=None
+                             directory=None,
+                             coarse_steps = 1e6,
+                             fine_steps = 1e6
                          ):
 
-    coarseSteps = 1e5+1
-    fineSteps = 1e5+1
-
     ret = None
     d_orig = os.getcwd()
     try:
+        directory = "/dev/shm/origami-job-robot_v1.9_bent--0f4axf4k"
         if directory is None:
             import tempfile
             directory = tempfile.mkdtemp(prefix='origami-%s-' % job_id, dir='/dev/shm/')
@@ -583,18 +583,18 @@ def run_simulation_protocol( output_name, job_id, json_data,
             output_prefix = "%s-0" % output_name
             full_output_prefix = "%s/%s" % (output_directory,output_prefix)
             ## TODO Remove long bonds
-            model.simulate( outputPrefix = output_prefix, numSteps=0.05*coarseSteps, **simargs )
+            model.simulate( outputPrefix = output_prefix, numSteps=0.05*coarse_steps, **simargs )
             coordinates = readArbdCoords('%s.0.restart' % full_output_prefix)
             output_prefix = "%s-1" % output_name
             model._update_segment_positions(coordinates)
             # model._clear_beads()
             # # model._generate_bead_model( 7, 4, False )
             # model._generate_bead_model( 7, 4, False )
-            model.simulate( outputPrefix = output_prefix, numSteps=0.95*coarseSteps, **simargs )
+            model.simulate( outputPrefix = output_prefix, numSteps=0.95*coarse_steps, **simargs )
         else:
             output_prefix = "%s-1" % output_name
             full_output_prefix = "%s/%s" % (output_directory,output_prefix)
-            model.simulate( outputPrefix = output_prefix, numSteps=coarseSteps, **simargs )
+            # model.simulate( outputPrefix = output_prefix, numSteps=coarse_steps, **simargs )
         coordinates = readArbdCoords('%s.0.restart' % full_output_prefix)
 
 
@@ -605,7 +605,7 @@ def run_simulation_protocol( output_name, job_id, json_data,
         model._update_segment_positions(coordinates)
         model._clear_beads()
         model._generate_bead_model( 1, 1, local_twist=True, escapable_twist=True )
-        model.simulate( outputPrefix = output_prefix, numSteps=fineSteps, **simargs )
+        # model.simulate( outputPrefix = output_prefix, numSteps=fine_steps, **simargs )
         coordinates = readAvgArbdCoords('%s.psf' % output_prefix,'%s.pdb' % output_prefix, '%s.0.dcd' % full_output_prefix, rmsdThreshold=1)
 
 
@@ -615,7 +615,7 @@ def run_simulation_protocol( output_name, job_id, json_data,
         model._update_segment_positions(coordinates)
         model._clear_beads()
         model._generate_bead_model( 1, 1, local_twist=True, escapable_twist=False )
-        model.simulate( outputPrefix = output_prefix, numSteps=fineSteps, **simargs )
+        model.simulate( outputPrefix = output_prefix, numSteps=fine_steps, **simargs )
         coordinates = readAvgArbdCoords('%s.psf' % output_name,'%s.pdb' % output_prefix, '%s.0.dcd' % full_output_prefix )
 
 
@@ -625,7 +625,7 @@ def run_simulation_protocol( output_name, job_id, json_data,
         model._update_segment_positions(coordinates)
         model._clear_beads()
         model._generate_atomic_model(scale=0.25)
-        model.atomic_simulate( outputPrefix = full_output_prefix )
+        model.atomic_simulate( outputPrefix = output_prefix )
 
         ret = directory
         
diff --git a/run.py b/run.py
index 03b1bb6..499f49c 100644
--- a/run.py
+++ b/run.py
@@ -18,7 +18,9 @@ if __name__ == '__main__':
         except:
             print("WARNING: skipping unreadable json file {}".format(f))
             continue
-        try:
-            run_simulation_protocol( out, "job-"+out+"-", data, gpu=0 )
-        except:
-            print("ERROR: failed to simulate {}".format(f))
+            
+        run_simulation_protocol( out, "job-"+out+"-", data, gpu=0 )
+        # try:
+        #     run_simulation_protocol( out, "job-"+out+"-", data, gpu=0 )
+        # except:
+        #     print("ERROR: failed to simulate {}".format(f))
diff --git a/segmentmodel.py b/segmentmodel.py
index 5875ecb..8ba4461 100644
--- a/segmentmodel.py
+++ b/segmentmodel.py
@@ -668,6 +668,9 @@ class Segment(ConnectableElement, Group):
         # conn_locs = self.get_contour_sorted_connections_and_locations()
         # locs = [A for c,A,B in conn_locs]
         # existing_beads = [l.particle for l in locs if l.particle is not None]
+        # if self.name == "51-2":
+        #     pdb.set_trace()
+
         existing_beads = {l.particle for l in self.locations if l.particle is not None}
         existing_beads = sorted( list(existing_beads), key=lambda b: b.get_contour_position(self) )
         
@@ -722,6 +725,8 @@ class Segment(ConnectableElement, Group):
             last = eb1
             for j in range(num_beads):
                 s = ds*(j+1) + s0
+                # if self.name in ("51-2","51-3"):
+                #     print(" adding bead at {}".format(s))
                 b = self._generate_one_bead(s,nts)
 
                 last.make_intrahelical_neighbor(b)
@@ -732,6 +737,8 @@ class Segment(ConnectableElement, Group):
 
         if eb2.parent == self:
             tmp_children.append(eb2)
+        # if self.name in ("51-2","51-3"):
+        #     pdb.set_trace()
         self._rebuild_children(tmp_children)
 
     def _regenerate_beads(self, max_nts_per_bead=4, ):
@@ -1292,7 +1299,11 @@ class SegmentModel(ArbdModel):
             beads = list(beads)
                         
             ## Skip beads that are None (some locations were not assigned a particle to avoid double-counting) 
-            beads = [b for b in beads if b is not None]
+            # beads = [b for b in beads if b is not None]
+            beads = list(filter(lambda x: x is not None, beads))
+            if isinstance(self, DoubleStrandedSegment):
+                beads = list(filter(lambda x: x.type_.name[0] == "D", beads))
+
             contours = []
             for b in beads:
                 try:
@@ -1316,10 +1327,31 @@ class SegmentModel(ArbdModel):
             """ Get positions """
             positions = bead_coordinates[ids,:].T
 
-            try:
-                tck, u = interpolate.splprep( positions, u=contours, s=0, k=3 )
-            except:
-                tck, u = interpolate.splprep( positions, u=contours, s=0, k=1 )
+            ret = interpolate.splprep( positions, u=contours, s=0, k=1, full_output=1 )
+            tck = ret[0][0]
+
+            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]
+            #     if ret[2] > 0:
+            #         pdb.set_trace()
+
+            # else:
+            #     try:
+            #         ret = interpolate.splprep( positions, u=contours, s=0, k=3, full_output=1 )
+            #         tck = ret[0][0]
+            #         if ret[2] > 0:
+            #             pdb.set_trace()
+            #     except:
+            #         ret = interpolate.splprep( positions, u=contours, s=0, k=1, full_output=1 )
+            #         tck = ret[0][0]
+            #         if ret[2] > 0:
+            #             pdb.set_trace()
 
             s.position_spline_params = tck
 
@@ -1404,6 +1436,8 @@ class SegmentModel(ArbdModel):
         ## Loop through all connections, generating beads at appropriate locations
         for c,A,B in self.get_connections("intrahelical"):
             s1,s2 = [l.container for l in (A,B)]
+            # if s1.name in ("51-2","51-3") and s2.name in ("51-2","51-3"):
+            #     pdb.set_trace()
             print("Working on {}".format(c))
             ## TODO be more elegant!
             # if isinstance(s1, DoubleStrandedSegment) and isinstance(s2, DoubleStrandedSegment) and A.on_fwd_strand == False: continue
@@ -1498,7 +1532,7 @@ class SegmentModel(ArbdModel):
             b1,b2 = [l.particle for l in (A,B)]
             if b1 is b2:
                 ## already handled by Segment._generate_beads
-                continue
+                pass
             else:
                 for b in (b1,b2): assert( b is not None )
                 b1.make_intrahelical_neighbor(b2)
-- 
GitLab