diff --git a/segmentmodel.py b/segmentmodel.py index 85faa80da3f0700e52b2b5fbcfb697bb369c6b54..3b57a50a30b788a7db6b4283ea5fa82562c6c555 100644 --- a/segmentmodel.py +++ b/segmentmodel.py @@ -20,16 +20,16 @@ TODO: + fix handling of crossovers for atomic representation + map to atomic representation + add nicks - - transform ssDNA nucleotides + + transform ssDNA nucleotides - shrink ssDNA - - shrink dsDNA backbone + + shrink dsDNA backbone + make orientation continuous - - sequence + + sequence - handle circular dna + ensure crossover bead potentials aren't applied twice + remove performance bottlenecks - test for large systems - - assign sequence + + assign sequence - ENM - rework Location class - remove recursive calls @@ -236,193 +236,16 @@ class SegmentParticle(PointParticle): self.intrahelical_neighbors.append(b) b.intrahelical_neighbors.append(self) - - # def get_nt_position(self,seg): - # if seg == self.parent: - # return seg.contour_to_nt_pos(self.contour_position) - # else: - # cl = [e for e in self.parent.get_connections_and_locations() if e[2].container is seg] - - # dc = [(self.contour_position - A.address)**2 for c,A,B in cl] - - # if len(dc) == 0: - # import pdb - # pdb.set_trace() - - # i = np.argmin(dc) - # c,A,B = cl[i] - # ## TODO: generalize, removing np.abs and conditional - # delta_nt = np.abs( A.container.contour_to_nt_pos(self.contour_position - A.address) ) - # B_nt_pos = seg.contour_to_nt_pos(B.address) - # if B.address < 0.5: - # return B_nt_pos-delta_nt - # else: - # return B_nt_pos+delta_nt - - def get_nt_position_old(self, seg, nearest_position_on_seg=None): - if seg == self.parent: - return seg.contour_to_nt_pos(self.contour_position) - else: - - def get_nt_pos(contour1, seg1, seg2): - cl = [e for e in seg1.get_connections_and_locations("intrahelical") if e[2].container is seg2] - dc = [(contour1 - A.address)**2 for c,A,B in cl] - if len(dc) == 0: return None,None - - i = np.argmin(dc) - c,A,B = cl[i] - - ## TODO: generalize, removing np.abs and conditional - # delta_nt = np.abs( seg1.contour_to_nt_pos(contour1 - A.address) ) - delta_nt = seg1.contour_to_nt_pos(contour1) - seg1.contour_to_nt_pos(A.address) - B_nt_pos = seg2.contour_to_nt_pos(B.address) - dirA = 2*(A.address > 0.5)-1 - dirB = 2*(B.address < 0.5)-1 - direction = dirA*dirB - if c.type_ in ("intrahelical",): - pass - elif c.type_ in ("crossover","terminal-crossover"): - direction *= -1 - else: - raise Exception("Unhandled connection type {}".format(c.type_)) - distance = np.abs(delta_nt)+np.abs(B_nt_pos) - return B_nt_pos + direction*delta_nt, distance - - pos,dist = get_nt_pos(self.contour_position, self.parent, seg) - if pos is None: - ## Particle is not directly connected - visited_segs = set(seg) - positions = [] - for l in self.locations: - if l.container == self.parent: continue - if l.connection.type_ != "intrahelical": continue - pos0,dist = get_nt_pos(self.contour_position, self.parent, l.container) - assert(pos0 is not None) - pos0 = l.container.nt_pos_to_contour(pos0) - pos,dist = get_nt_pos( pos0, l.container, seg ) - if pos is not None: - positions.append( pos ) - assert( len(positions) > 0 ) - if len(positions) > 1: - import pdb - pdb.set_trace() - pos = positions[0] - return pos - def get_nt_position(self, seg): """ Returns the "address" of the nucleotide relative to seg in nucleotides, taking the shortest (intrahelical) contour length route to seg """ - if seg == self.parent: return seg.contour_to_nt_pos(self.contour_position) else: pos = self.get_contour_position(seg) return seg.contour_to_nt_pos(pos) - - # def get_nt_pos(contour1, seg1, seg2): - # cl = [e for e in seg1.get_connections_and_locations("intrahelical") if e[2].container is seg2] - # dc = [(contour1 - A.address)**2 for c,A,B in cl] - # if len(dc) == 0: return None,None - - # i = np.argmin(dc) - # c,A,B = cl[i] - - # ## TODO: generalize, removing np.abs and conditional - # # delta_nt = np.abs( seg1.contour_to_nt_pos(contour1 - A.address) ) - # delta_nt = seg1.contour_to_nt_pos(contour1) - seg1.contour_to_nt_pos(A.address) - # B_nt_pos = seg2.contour_to_nt_pos(B.address) - # dirA = 2*(A.address > 0.5)-1 - # dirB = 2*(B.address < 0.5)-1 - # direction = dirA*dirB - # distance = np.abs(delta_nt)+np.abs(B_nt_pos) - # return B_nt_pos + direction*delta_nt, distance - # pos,dist = get_nt_pos(self.contour_position, self.parent, seg) - - cutoff = 30*3 - target_seg = seg - - ## depth-first search - ## TODO cache distances to nearby locations? - def descend_search_tree(seg, contour_in_seg, distance=0, visited_segs=[]): - nonlocal cutoff - - if seg == target_seg: - # pdb.set_trace() - ## Found a segment in our target - sign = (contour_in_seg == 1) - (contour_in_seg == 0) - assert( sign in (1,-1) ) - if distance < cutoff: # TODO: check if this does anything - cutoff = distance - return [[distance, contour_in_seg+sign*seg.nt_pos_to_contour(distance)]] - if distance > cutoff: - return None - - ret_list = [] - ## Find intrahelical locations in seg that we might pass through - for c,A,B in seg.get_connections_and_locations("intrahelical"): - if B.container in visited_segs: continue - dx = seg.contour_to_nt_pos( np.abs(A.address-contour_in_seg), - round_nt=False) - results = descend_search_tree( B.container, B.address, - distance+dx, visited_segs + [seg] ) - if results is not None: - ret_list.extend( results ) - return ret_list - - results = descend_search_tree(self.parent, self.contour_position) - if results is None or len(results) == 0: - raise Exception("Could not find location in segment") # TODO better error - return sorted(results,key=lambda x:x[0])[0][1] - - # def descend_search_tree(seg, target_seg, distance, from_contour): - # if distance > cutoff: - # return None - # ret_list = [] - - # ## Find intrahelical locations in seg that we might pass through - # for c,A,B in seg.get_connections_and_locations("intrahelical"): - # if np.isclose(A.address, from_contour): continue - # if B.container == target_seg: - # ## Found a segment in our target - # ... - # else: - # ret_list.append( [B.container, seg.num_nts] ] - - - # curr_seg = seg - # dist = 0 - # for c,A,B in curr_seg.get_connections_and_locations("intrahelical"): - # while - - # cl = [e for e in seg1.get_connections_and_locations("intrahelical") if e[2].container is seg2] - - - - # if dist is not None and dist < closest_distance: - # closest_distance - # if pos is None: - # ## Particle is not directly connected - # visited_segs = set(seg) - # positions = [] - # for l in self.locations: - # if l.container == self.parent: continue - # if l.connection.type_ != "intrahelical": continue - # pos0,dist = get_nt_pos(self.contour_position, self.parent, l.container) - # assert(pos0 is not None) - # pos0 = l.container.nt_pos_to_contour(pos0) - # pos,dist = get_nt_pos( pos0, l.container, seg ) - # if pos is not None: - # positions.append( pos ) - # assert( len(positions) > 0 ) - # if len(positions) > 1: - # import pdb - # pdb.set_trace() - # pos = positions[0] - # return pos - - def get_contour_position(self,seg): if seg == self.parent: return self.contour_position