Various changes to segmentmodel, mostly involving deletion of DNA

mrdna/segmentmodel.py

@@ -533,6 +533,13 @@ class Segment(ConnectableElement, Group):
         last_nt = max((first_nt,last_nt))
         fist_nt = tmp
 
+        ## Sanity check
+        for l in self.locations:
+            nt = self.contour_to_nt_pos(l.address)
+            if not (np.isclose(nt,int(np.round(nt))) or l.address in (0,1)):
+                raise Exception
+            l.old_address = (l.address, self.num_nt)
+      
         if first_nt < 0 or first_nt > self.num_nt-2:
             raise ValueError("Attempted to remove DNA from {} starting at an invalid location {}".format(self, first_nt))
         if last_nt < 1 or last_nt > self.num_nt-1:
@@ -544,7 +551,7 @@ class Segment(ConnectableElement, Group):
             nt = self.contour_to_nt_pos(contours)
             first = first_nt if first_nt > 0 else -np.inf
             last = last_nt if last_nt < self.num_nt-1 else np.inf
-            bad_ids = (nt >= first) * (nt <= last)
+            bad_ids = (nt >= first-0.5) * (nt <= last+0.5)
             nt[nt>last] = nt[nt>last]-removed_nt
             nt[bad_ids] = np.nan
             return self.nt_pos_to_contour(nt)* self.num_nt/num_nt
@@ -561,11 +568,24 @@ class Segment(ConnectableElement, Group):
         if np.any(np.isnan(new_c)) and not remove_locations:
             raise Exception("Attempted to remove DNA containing locations from {} between {} and {}".format(self,first_nt,last_nt))
 
+        ## get locations at ends 
+        end_locs = [self.get_location_at( nt, is_fwd )
+                    for nt in (first_nt, last_nt) for is_fwd in (True,False)]
+
+        # print("Removing DNA:",first_nt,last_nt, self.num_nt, removed_nt)
+        # print(end_locs)
+        # print([l.connection for l in end_locs])
         for l,v in zip(list(self.locations),new_c):
             if np.isnan(v):
                 l.delete()
             else:
-                l.address = v
+                if l.address not in (0,1):
+                    # print("  Updating loc:", l.address, v, self.contour_to_nt_pos(l.address), v*num_nt-0.5 )
+                    l.address = v
+                else:
+                    idx = int((l.address == 1)*2 + l.on_fwd_strand)
+                    if end_locs[idx] is not None and l is not end_locs[idx]:
+                        l.delete()
 
         assert( len(self.locations) == np.logical_not(np.isnan(new_c)).sum() )
         tmp_nt = np.array([l.address*num_nt+0.5 for l in self.locations])
@@ -584,7 +604,9 @@ class Segment(ConnectableElement, Group):
             new_p = get_coord_function(u[ids])
             set_spline_function( new_u, new_p )
 
-        _rescale_splines(self.position_spline_params[1],
+        # u = self.position_spline_params[1]
+        u = np.linspace(0,1,self.num_nt*2+2)
+        _rescale_splines(u,
                          self.contour_to_position,
                          self.set_splines)
 
@@ -593,12 +615,26 @@ class Segment(ConnectableElement, Group):
                 os = [self.contour_to_orientation(c) for c in contours]
                 qs = [quaternion_from_matrix( o ) for o in os]
                 return np.array(qs)
-            _rescale_splines(self.quaternion_spline_params[1],
+            # u = self.quaternion_spline_params[1],
+            _rescale_splines(u,
                              _contour_to_quaternion,
                              self.set_orientation_splines)
 
+        if '_cadnano_bp_to_zidx' in self.__dict__:
+            assert( len(self._cadnano_bp_to_zidx) == self.num_nt )
+            self._cadnano_bp_to_zidx = self._cadnano_bp_to_zidx[:first_nt] + self._cadnano_bp_to_zidx[last_nt+1:]
+            assert( len(self._cadnano_bp_to_zidx) == num_nt )
+
         self.num_nt = num_nt
 
+        ## Sanity check
+        for l in self.locations:
+            nt = self.contour_to_nt_pos(l.address)
+            if not (np.isclose(nt,int(np.round(nt))) or l.address in (0,1)):
+                raise Exception
+
+        
+
     def delete(self):
         for c,loc,other in self.get_connections_and_locations():
             c.delete()
@@ -675,7 +711,7 @@ class Segment(ConnectableElement, Group):
     def contour_to_nt_pos(self, contour_pos, round_nt=False):
         nt = contour_pos*(self.num_nt) - 0.5
         if round_nt:
-            assert( np.isclose(np.around(nt),nt) )
+            assert( np.isclose(np.around(nt),nt) and contour_pos not in (0,1)  )
             nt = np.around(nt)
         return nt
 
@@ -1972,6 +2008,7 @@ class SegmentModel(ArbdModel):
         assert( isinstance(segment, Segment) )
         if copy:
             segment = deepcopy(segment)
+        segment.parent = self
         self.segments.append(segment)
         self._clear_beads()