diff --git a/mrdna/segmentmodel.py b/mrdna/segmentmodel.py index 16872d4a9448af9bbf32f97a4d38f042b74d6580..f642d6e118a61c85f59e4b2033cae393c6cf205c 100644 --- a/mrdna/segmentmodel.py +++ b/mrdna/segmentmodel.py @@ -286,7 +286,11 @@ class SegmentParticle(PointParticle): ret_list = [] hist_list = [] ## Find intrahelical locations in seg that we might pass through - for c,A,B in seg.get_connections_and_locations("intrahelical"): + conn_locs = seg.get_connections_and_locations("intrahelical") + if isinstance(target_seg, SingleStrandedSegment): + tmp = seg.get_connections_and_locations("sscrossover") + conn_locs = conn_locs + list(filter(lambda x: x[2].container == target_seg, tmp)) + for c,A,B in conn_locs: if B.container in visited_segs: continue dx = seg.contour_to_nt_pos( A.address, round_nt=False ) - seg.contour_to_nt_pos( contour_in_seg, round_nt=False) dx = np.abs(dx) @@ -742,7 +746,10 @@ class Segment(ConnectableElement, Group): existing_beads0 = {l.particle for l in self.locations if l.particle is not None} existing_beads = sorted( list(existing_beads0), key=lambda b: b.get_contour_position(self) ) - + + if self.num_nt == 1 and all([l.particle is not None for l in self.locations]): + return + for b in existing_beads: assert(b.parent is not None) @@ -1458,12 +1465,13 @@ class SegmentModel(ArbdModel): # pdb.set_trace() cabs = s.get_connections_and_locations("intrahelical") + if isinstance(s,SingleStrandedSegment): + cabs = cabs + [[c,A,B] for c,A,B in s.get_connections_and_locations("sscrossover") if A.address == 0 or A.address == 1] if np.any( [B.particle is None for c,A,B in cabs] ): print( "WARNING: none type found in connection, skipping" ) cabs = [e for e in cabs if e[2].particle is not None] + beads = set(s.beads + [A.particle for c,A,B in cabs]) - if len(beads) <= 1: - pdb.set_trace() ## Add nearby beads for c,A,B in cabs: @@ -1484,6 +1492,7 @@ class SegmentModel(ArbdModel): beads = list(filter(lambda x: x.type_.name[0] == "D", beads)) contours = [] + remove = [] for b in beads: try: # if s.name == "61-1" and b.parent.name == "60-1": @@ -1492,9 +1501,9 @@ class SegmentModel(ArbdModel): contours.append(b.get_contour_position(s)) except: ## ignore failed attempts - beads.remove(b) - - contours = [b.get_contour_position(s) for b in beads] + remove.append(b) + beads = list(filter(lambda x: x not in remove, 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-8 ): @@ -1639,15 +1648,6 @@ 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 ("49-3","49-4") and s2.name in ("49-3","49-4"): - # 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 - # if isinstance(s1, DoubleStrandedSegment) and isinstance(s2, DoubleStrandedSegment) and A.on_fwd_strand == False: continue - ## if A.on_fwd_strand == False: continue # TODO verify this avoids double-counting - ## TODO determine whether any logic is needed to prevent double-counting assert( A.particle is None ) assert( B.particle is None ) @@ -1655,10 +1655,8 @@ class SegmentModel(ArbdModel): ## TODO: offload the work here to s1 # TODOTODO a1,a2 = [l.address for l in (A,B)] - # a1,a2 = [a - s.nt_pos_to_contour(0.5) if a == 0 else a + s.nt_pos_to_contour(0.5) for a,s in zip((a1,a2),(s1,s2))] for a in (a1,a2): assert( np.isclose(a,0) or np.isclose(a,1) ) - # a1,a2 = [a - s.nt_pos_to_contour(0) if a == 0 else a + s.nt_pos_to_contour(0) for a,s in zip((a1,a2),(s1,s2))] ## TODO improve this for combinations of ssDNA and dsDNA (maybe a1/a2 should be calculated differently) @@ -2364,7 +2362,6 @@ class SegmentModel(ArbdModel): def add_strand_if_needed(seg,is_fwd): history = [] if not strands_cover_segment(seg, is_fwd): - pdb.set_trace() pos = nt = find_nt_not_in_strand(seg, is_fwd) s = Strand(is_circular = True) history = _recursively_build_strand(s, seg, pos, is_fwd)