fixes #268: ignore loopouts when relaxing Helix rolls

scadnano/scadnano.py

@@ -1763,35 +1763,44 @@ def backbone_angle_at_offset(self, offset: int, forward: bool, geometry: Geometr
         angle %= 360
         return angle
 
-    def crossovers(self) -> List[Tuple[int, int, bool]]:
+    def crossover_addresses(self) -> List[Tuple[int, int, bool]]:
         """
         :return:
-            list of triples (`offset`, `helix_idx`, `forward`) of all crossovers incident to this
+            list of triples (`helix_idx`, `offset`, `forward`) of all crossovers incident to this
             :any:`Helix`, where `offset` is the offset of the crossover and `helix_idx` is the
             :data:`Helix.idx` of the other :any:`Helix` incident to the crossover.
         """
-        crossovers: List[Tuple[int, int, bool]] = []
+        addresses: List[Tuple[int, int, bool]] = []
         for domain in self.domains:
             strand = domain.strand()
             domains_on_strand = strand.bound_domains()
             num_domains = len(domains_on_strand)
             domain_idx = domains_on_strand.index(domain)
+            domain_idx_in_substrands = strand.domains.index(domain)
 
             # if not first domain, then there is a crossover to the previous domain
             if domain_idx > 0:
-                offset = domain.offset_5p()
-                other_domain = domains_on_strand[domain_idx - 1]
-                other_helix_idx = other_domain.helix
-                crossovers.append((offset, other_helix_idx, domain.forward))
+                # ... unless there's a loopout between them
+                previous_substrand = strand.domains[domain_idx_in_substrands - 1]
+                if isinstance(previous_substrand, Domain):
+                    offset = domain.offset_5p()
+                    other_domain = domains_on_strand[domain_idx - 1]
+                    assert previous_substrand == other_domain
+                    other_helix_idx = other_domain.helix
+                    addresses.append((other_helix_idx, offset, domain.forward))
 
             # if not last domain, then there is a crossover to the next domain
             if domain_idx < num_domains - 1:
-                offset = domain.offset_3p()
-                other_domain = domains_on_strand[domain_idx + 1]
-                other_helix_idx = other_domain.helix
-                crossovers.append((offset, other_helix_idx, domain.forward))
+                # ... unless there's a loopout between them
+                next_substrand = strand.domains[domain_idx_in_substrands + 1]
+                if isinstance(next_substrand, Domain):
+                    offset = domain.offset_3p()
+                    other_domain = domains_on_strand[domain_idx + 1]
+                    assert next_substrand == other_domain
+                    other_helix_idx = other_domain.helix
+                    addresses.append((other_helix_idx, offset, domain.forward))
 
-        return crossovers
+        return addresses
 
     def relax_roll(self, helices: Dict[int, Helix], grid: Grid, geometry: Geometry) -> None:
         """
@@ -1806,7 +1815,7 @@ def compute_relaxed_roll(self, helices: Dict[int, Helix], grid: Grid, geometry:
         rather than changing the field :data:`Helix.roll`.
         """
         angles = []
-        for offset, helix_idx, forward in self.crossovers():
+        for helix_idx, offset, forward in self.crossover_addresses():
             other_helix = helices[helix_idx]
             angle_of_other_helix = angle_from_helix_to_helix(self, other_helix, grid, geometry)
             crossover_angle = self.backbone_angle_at_offset(offset, forward, geometry)

tests/scadnano_tests.py

@@ -8215,6 +8215,78 @@ def test_3_helix_2_crossovers(self) -> None:
         self.assertAlmostEqual(exp_h1_roll, design3h.helices[1].roll)
         self.assertAlmostEqual(exp_h2_roll, design3h.helices[2].roll)
 
+    def test_3_helix_2_crossovers_1_loopout_crossovers_method(self) -> None:
+        '''
+          0         1          2
+          012345678901234567890
+        0 [---+[------+[------+
+              |       |        \
+        1 [---+       |        |
+                      |        /
+        2      <------+<------+
+        '''
+        helices = [sc.Helix(max_offset=60) for _ in range(3)]
+        helices[2].grid_position = (1, 0)
+        design3h = sc.Design(helices=helices, grid=sc.square)
+        design3h.draw_strand(0, 0).move(5).cross(1).move(-5)
+        design3h.draw_strand(0, 5).move(8).cross(2).move(-8)
+        design3h.draw_strand(0, 13).move(8).loopout(2, 3).move(-8)
+
+        crossover_addresses_h0 = helices[0].crossover_addresses()
+        crossover_addresses_h1 = helices[1].crossover_addresses()
+        crossover_addresses_h2 = helices[2].crossover_addresses()
+
+        self.assertEqual(2, len(crossover_addresses_h0))
+        self.assertEqual(1, len(crossover_addresses_h1))
+        self.assertEqual(1, len(crossover_addresses_h2))
+
+        self.assertEqual((1, 4, True), crossover_addresses_h0[0])
+        self.assertEqual((2, 12, True), crossover_addresses_h0[1])
+
+        self.assertEqual((0, 4, False), crossover_addresses_h1[0])
+
+        self.assertEqual((0, 12, False), crossover_addresses_h2[0])
+
+    def test_3_helix_2_crossovers_1_loopout(self) -> None:
+        '''
+          0         1          2
+          012345678901234567890
+        0 [---+[------+[------+
+              |       |        \
+        1 [---+       |        |
+                      |        /
+        2      <------+<------+
+        '''
+        helices = [sc.Helix(max_offset=60) for _ in range(3)]
+        helices[2].grid_position = (1, 0)
+        design3h = sc.Design(helices=helices, grid=sc.square)
+        design3h.draw_strand(0, 0).move(5).cross(1).move(-5)
+        design3h.draw_strand(0, 5).move(8).cross(2).move(-8)
+        design3h.draw_strand(0, 13).move(8).loopout(2, 3).move(-8)
+        f1 = 4 / 10.5
+        f2 = 12 / 10.5
+        a1 = f1 * 360 % 360
+        a2 = f2 * 360 % 360
+
+        # rules for angles:
+        # - add 150 if on reverse strand to account of minor groove
+        # - subtract angle of helix crossover is connecting to
+
+        ave_h0 = (a1 - 180 + a2 - 90) / 2  # helix 1 at 180 degrees, helix 2 at 90 degrees
+        exp_h0_roll = (-ave_h0) % 360
+
+        ave_h1 = a1 + 150  # helix 0 at 0 degrees relative to helix 1
+        exp_h1_roll = (-ave_h1) % 360
+
+        ave_h2 = a2 + 150 - (-90)  # helix 0 at -90 degrees relative to helix 2
+        exp_h2_roll = (-ave_h2) % 360
+
+        design3h.relax_helix_rolls()
+
+        self.assertAlmostEqual(exp_h0_roll, design3h.helices[0].roll)
+        self.assertAlmostEqual(exp_h1_roll, design3h.helices[1].roll)
+        self.assertAlmostEqual(exp_h2_roll, design3h.helices[2].roll)
+
     def test_3_helix_6_crossover(self) -> None:
         '''
           0         1         2         3         4         5         6
@@ -8275,7 +8347,7 @@ def test_2_helix_no_crossover(self) -> None:
 
         1 <---]<--------]
         '''
-        initial_roll = 30
+        initial_roll = 30.0
         helices = [sc.Helix(max_offset=60, roll=initial_roll) for _ in range(2)]
         design2h = sc.Design(helices=helices, grid=sc.square)
         design2h.draw_strand(0, 0).move(5)
@@ -8313,11 +8385,11 @@ def test_2_helix_3_crossover(self) -> None:
     def test_helix_crossovers(self) -> None:
         ############################################
         # 3-helix design with 3 strands
-        xs0 = self.design3helix3strand.helices[0].crossovers()
+        xs0 = self.design3helix3strand.helices[0].crossover_addresses()
         self.assertEqual(len(xs0), 3)
-        o0, h0, f0 = xs0[0]
-        o1, h1, f1 = xs0[1]
-        o2, h2, f2 = xs0[2]
+        h0, o0, f0 = xs0[0]
+        h1, o1, f1 = xs0[1]
+        h2, o2, f2 = xs0[2]
         self.assertEqual(o0, 4)
         self.assertEqual(o1, 14)
         self.assertEqual(o2, 26)
@@ -8328,31 +8400,31 @@ def test_helix_crossovers(self) -> None:
         self.assertEqual(f1, True)
         self.assertEqual(f2, True)
 
-        xs1 = self.design3helix3strand.helices[1].crossovers()
+        xs1 = self.design3helix3strand.helices[1].crossover_addresses()
         self.assertEqual(len(xs1), 2)
-        o0, h0, f0 = xs1[0]
-        o1, h1, f1 = xs1[1]
+        h0, o0, f0 = xs1[0]
+        h1, o1, f1 = xs1[1]
         self.assertEqual(o0, 4)
         self.assertEqual(o1, 26)
         self.assertEqual(h0, 0)
         self.assertEqual(h1, 0)
         self.assertEqual(f0, False)
         self.assertEqual(f1, False)
 
-        xs2 = self.design3helix3strand.helices[2].crossovers()
+        xs2 = self.design3helix3strand.helices[2].crossover_addresses()
         self.assertEqual(len(xs2), 1)
-        o0, h0, f0 = xs2[0]
+        h0, o0, f0 = xs2[0]
         self.assertEqual(o0, 14)
         self.assertEqual(h0, 0)
         self.assertEqual(h0, False)
 
         ############################################
         # 2-helix design
-        xs0 = self.design2h.helices[0].crossovers()
+        xs0 = self.design2h.helices[0].crossover_addresses()
         self.assertEqual(len(xs0), 3)
-        o0, h0, f0 = xs0[0]
-        o1, h1, f1 = xs0[1]
-        o2, h2, f2 = xs0[2]
+        h0, o0, f0 = xs0[0]
+        h1, o1, f1 = xs0[1]
+        h2, o2, f2 = xs0[2]
         self.assertEqual(o0, 4)
         self.assertEqual(o1, 14)
         self.assertEqual(o2, 26)
@@ -8363,11 +8435,11 @@ def test_helix_crossovers(self) -> None:
         self.assertEqual(h1, True)
         self.assertEqual(h2, True)
 
-        xs1 = self.design2h.helices[1].crossovers()
+        xs1 = self.design2h.helices[1].crossover_addresses()
         self.assertEqual(len(xs1), 3)
-        o0, h0, f0 = xs1[0]
-        o1, h1, f1 = xs1[1]
-        o2, h2, f2 = xs1[2]
+        h0, o0, f0 = xs1[0]
+        h1, o1, f1 = xs1[1]
+        h2, o2, f2 = xs1[2]
         self.assertEqual(o0, 4)
         self.assertEqual(o1, 14)
         self.assertEqual(o2, 26)