Repository: UC-Davis-molecular-computing/scadnano-python-package

Issue:
deal with inter-helix-group crossovers when relaxing helix rolls
Currently, when relaxing helix rolls, all inter-helix crossovers are accounted for, but actually the calculation of the relative angle from one helix to another (based on the position of their 0 offset) only makes sense if the helices are in the same helix group, i.e., are parallel. If they are in different helix groups, we should do one of two things (possibly configurable with an option):

1. (easy) ignore those crossovers
2. (hard) do a more detailed calculation of the angle of the other helix, based not only on the helix's origin position, but also on the offsets of the two points on the helices the crossover connects.

Discussion:

Patch:
diff --git a/scadnano/scadnano.py b/scadnano/scadnano.py
index 14b96ac..6cc4a22 100644
--- a/scadnano/scadnano.py
+++ b/scadnano/scadnano.py
@@ -1764,15 +1764,32 @@ class Helix(_JSONSerializable):
         angle %= 360
         return angle
 
-    def crossover_addresses(self, allow_intrahelix=True) -> List[Tuple[int, int, bool]]:
+    def crossover_addresses(self,
+                            helices: Dict[int, Helix],
+                            allow_intrahelix: bool = True,
+                            allow_intergroup: bool = True) -> List[Tuple[int, int, bool]]:
         """
+        :param helices:
+            The dict of helices in which this :any:`Helix` is contained, that contains other helices
+            to which it might be connected by crossovers.
         :param allow_intrahelix:
             if ``False``, then do not return crossovers to the same :any:`Helix` as this :any:`Helix`
+        :param allow_intergroup:
+            if ``False``, then do not return crossovers to a :any:`Helix` in a different helix group
+            as this :any:`Helix`
         :return:
             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.
         """
+
+        def allow_crossover_to(other_helix: Helix) -> bool:
+            if not allow_intrahelix and other_helix.idx == self.idx:
+                return False
+            if not allow_intergroup and other_helix.group != self.group:
+                return False
+            return True
+
         addresses: List[Tuple[int, int, bool]] = []
         for domain in self.domains:
             strand = domain.strand()
@@ -1790,7 +1807,8 @@ class Helix(_JSONSerializable):
                     other_domain = domains_on_strand[domain_idx - 1]
                     assert previous_substrand == other_domain
                     other_helix_idx = other_domain.helix
-                    if allow_intrahelix or other_helix_idx != self.idx:
+                    other_helix = helices[other_helix_idx]
+                    if allow_crossover_to(other_helix):
                         addresses.append((other_helix_idx, offset, domain.forward))
 
             # if not last domain, then there is a crossover to the next domain
@@ -1802,7 +1820,8 @@ class Helix(_JSONSerializable):
                     other_domain = domains_on_strand[domain_idx + 1]
                     assert next_substrand == other_domain
                     other_helix_idx = other_domain.helix
-                    if allow_intrahelix or other_helix_idx != self.idx:
+                    other_helix = helices[other_helix_idx]
+                    if allow_crossover_to(other_helix):
                         addresses.append((other_helix_idx, offset, domain.forward))
 
         return addresses
@@ -1820,7 +1839,7 @@ class Helix(_JSONSerializable):
         without actually altering the field :data:`Helix.roll`.
         """
         angles = []
-        addresses = self.crossover_addresses(allow_intrahelix=False)
+        addresses = self.crossover_addresses(helices, allow_intrahelix=False)
         for helix_idx, offset, forward in addresses:
             other_helix = helices[helix_idx]
             angle_of_other_helix = angle_from_helix_to_helix(self, other_helix, grid, geometry)


Test patch:
diff --git a/tests/scadnano_tests.py b/tests/scadnano_tests.py
index 7ed6cbd..b50c487 100644
--- a/tests/scadnano_tests.py
+++ b/tests/scadnano_tests.py
@@ -8149,7 +8149,7 @@ class TestHelixRollRelax(unittest.TestCase):
           0123456789012345678901234567890123456789
         0 [---+[--------+[----------+
               |         |           |
-        1 [---+<--------+<----------+
+        1 <---+<--------+<----------+
         '''
         self.design2h = sc.Design(helices=[sc.Helix(max_offset=50) for _ in range(2)])
         # helix 0 forward
@@ -8161,7 +8161,7 @@ class TestHelixRollRelax(unittest.TestCase):
           0123456789012345678901234567890123456789
         0 [---+[--------+[----------+
               |         |           |
-        1 [---+         |<----------+
+        1 <---+         |<----------+
                         |
         2      <--------+
         '''
@@ -8177,7 +8177,7 @@ class TestHelixRollRelax(unittest.TestCase):
           012345678901234
         0 [---+[------+
               |       |
-        1 [---+       |
+        1 <---+       |
                       |
         2      <------+
         '''
@@ -8216,7 +8216,7 @@ class TestHelixRollRelax(unittest.TestCase):
           012345678901234567890
         0 [---+[------+[------+
               |       |        \
-        1 [---+       |        |
+        1 <---+       |        |
                       |        /
         2      <------+<------+
         '''
@@ -8227,9 +8227,9 @@ class TestHelixRollRelax(unittest.TestCase):
         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()
+        crossover_addresses_h0 = helices[0].crossover_addresses(design3h.helices)
+        crossover_addresses_h1 = helices[1].crossover_addresses(design3h.helices)
+        crossover_addresses_h2 = helices[2].crossover_addresses(design3h.helices)
 
         self.assertEqual(2, len(crossover_addresses_h0))
         self.assertEqual(1, len(crossover_addresses_h1))
@@ -8248,7 +8248,7 @@ class TestHelixRollRelax(unittest.TestCase):
           012345678901234567890
         0 [---+[------+[------+
               |       |        \
-        1 [---+       |        |
+        1 <---+       |        |
                       |        /
         2      <------+<------+
         '''
@@ -8288,7 +8288,7 @@ class TestHelixRollRelax(unittest.TestCase):
           012345678901234567890123456789012345678901234567890123456789
         0 [---+[--------+[----------+[------+[--------+[--------+
               |         |           |       |         |         |
-        1 [---+<--------+<----------+       |         |         |
+        1 <---+<--------+<----------+       |         |         |
                                             |         |         |
         2                            <------+<--------+<--------+
         '''
@@ -8436,7 +8436,7 @@ class TestHelixRollRelax(unittest.TestCase):
         self.assertAlmostEqual(exp_h1_roll, design2h.helices[1].roll)
 
     def test_helix_crossover_addresses_3_helix_3_strand(self) -> None:
-        xs0 = self.design3helix3strand.helices[0].crossover_addresses()
+        xs0 = self.design3helix3strand.helices[0].crossover_addresses(self.design3helix3strand.helices)
         self.assertEqual(len(xs0), 3)
         h0, o0, f0 = xs0[0]
         h1, o1, f1 = xs0[1]
@@ -8451,7 +8451,7 @@ class TestHelixRollRelax(unittest.TestCase):
         self.assertEqual(f1, True)
         self.assertEqual(f2, True)
 
-        xs1 = self.design3helix3strand.helices[1].crossover_addresses()
+        xs1 = self.design3helix3strand.helices[1].crossover_addresses(self.design3helix3strand.helices)
         self.assertEqual(len(xs1), 2)
         h0, o0, f0 = xs1[0]
         h1, o1, f1 = xs1[1]
@@ -8462,7 +8462,7 @@ class TestHelixRollRelax(unittest.TestCase):
         self.assertEqual(f0, False)
         self.assertEqual(f1, False)
 
-        xs2 = self.design3helix3strand.helices[2].crossover_addresses()
+        xs2 = self.design3helix3strand.helices[2].crossover_addresses(self.design3helix3strand.helices)
         self.assertEqual(len(xs2), 1)
         h0, o0, f0 = xs2[0]
         self.assertEqual(o0, 14)
@@ -8470,7 +8470,7 @@ class TestHelixRollRelax(unittest.TestCase):
         self.assertEqual(h0, False)
 
     def test_helix_crossover_addresses_2_helix_3_strand(self) -> None:
-        xs0 = self.design2h.helices[0].crossover_addresses()
+        xs0 = self.design2h.helices[0].crossover_addresses(self.design2h.helices)
         self.assertEqual(len(xs0), 3)
         h0, o0, f0 = xs0[0]
         h1, o1, f1 = xs0[1]
@@ -8485,7 +8485,7 @@ class TestHelixRollRelax(unittest.TestCase):
         self.assertEqual(h1, True)
         self.assertEqual(h2, True)
 
-        xs1 = self.design2h.helices[1].crossover_addresses()
+        xs1 = self.design2h.helices[1].crossover_addresses(self.design2h.helices)
         self.assertEqual(len(xs1), 3)
         h0, o0, f0 = xs1[0]
         h1, o1, f1 = xs1[1]
@@ -8505,12 +8505,64 @@ class TestHelixRollRelax(unittest.TestCase):
           0123456789012345678901234567890123456789
         0 [---+[--------+[----------+[--+c+-->
               |         |           |
-        1 [---+<--------+<----------+<--+c+--]
+        1 <---+<--------+<----------+<--+c+--]
         '''
         self.design2h.draw_strand(0, 27).move(4).cross(0, 32).move(4)
         self.design2h.draw_strand(1, 36).move(-4).cross(1, 31).move(-4)
 
-        xs0 = self.design2h.helices[0].crossover_addresses(allow_intrahelix=False)
+        xs0 = self.design2h.helices[0].crossover_addresses(self.design2h.helices, allow_intrahelix=False)
+        self.assertEqual(len(xs0), 3)
+        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)
+        self.assertEqual(h0, 1)
+        self.assertEqual(h1, 1)
+        self.assertEqual(h2, 1)
+        self.assertEqual(f0, True)
+        self.assertEqual(f1, True)
+        self.assertEqual(f2, True)
+
+        xs1 = self.design2h.helices[1].crossover_addresses(self.design2h.helices, allow_intrahelix=False)
+        self.assertEqual(len(xs1), 3)
+        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)
+        self.assertEqual(h0, 0)
+        self.assertEqual(h1, 0)
+        self.assertEqual(h2, 0)
+        self.assertEqual(f0, False)
+        self.assertEqual(f1, False)
+        self.assertEqual(f2, False)
+
+    def test_helix_crossover_addresses_2_helix_disallow_intergroup_crossovers(self) -> None:
+        '''         1         2         3
+          0123456789012345678901234567890123456789
+        0 [---+[--------+[----------+[--+
+              |         |           |   |
+        1 <---+<--------+<----------+   |
+                                        |
+        group 2                         |
+        2                            <--+
+        '''
+        group2name = 'group 2'
+        group2 = sc.HelixGroup(position=sc.Position3D(x=0, y=-10, z=0), grid=sc.square,
+                               helices_view_order=[2])
+        self.design2h.groups[group2name] = group2
+        self.design2h.add_helix(2, sc.Helix(max_offset=50, group=group2name))
+        self.design2h.draw_strand(0, 27).move(4).cross(2).move(-4)
+
+        self.assertEqual(2, len(self.design2h.groups))
+
+        exp_group_names = list(self.design2h.groups.keys())
+        self.assertListEqual([sc.default_group_name, group2name], exp_group_names)
+
+        xs0 = self.design2h.helices[0].crossover_addresses(self.design2h.helices, allow_intergroup=False)
         self.assertEqual(len(xs0), 3)
         h0, o0, f0 = xs0[0]
         h1, o1, f1 = xs0[1]
@@ -8525,7 +8577,7 @@ class TestHelixRollRelax(unittest.TestCase):
         self.assertEqual(f1, True)
         self.assertEqual(f2, True)
 
-        xs1 = self.design2h.helices[1].crossover_addresses(allow_intrahelix=False)
+        xs1 = self.design2h.helices[1].crossover_addresses(self.design2h.helices, allow_intergroup=False)
         self.assertEqual(len(xs1), 3)
         h0, o0, f0 = xs1[0]
         h1, o1, f1 = xs1[1]
@@ -8540,6 +8592,9 @@ class TestHelixRollRelax(unittest.TestCase):
         self.assertEqual(f1, False)
         self.assertEqual(f2, False)
 
+        xs2 = self.design2h.helices[2].crossover_addresses(self.design2h.helices, allow_intergroup=False)
+        self.assertEqual(len(xs2), 0)
+
     def test_helix_crossover_addresses_2_helix_allow_intrahelix_crossovers(self) -> None:
         '''         1         2         3
           0123456789012345678901234567890123456789
@@ -8550,7 +8605,7 @@ class TestHelixRollRelax(unittest.TestCase):
         self.design2h.draw_strand(0, 27).move(4).cross(0, 32).move(4)
         self.design2h.draw_strand(1, 36).move(-4).cross(1, 31).move(-4)
 
-        xs0 = self.design2h.helices[0].crossover_addresses(allow_intrahelix=True)
+        xs0 = self.design2h.helices[0].crossover_addresses(self.design2h.helices, allow_intrahelix=True)
         self.assertEqual(len(xs0), 5)
         h0, o0, f0 = xs0[0]
         h1, o1, f1 = xs0[1]
@@ -8573,7 +8628,7 @@ class TestHelixRollRelax(unittest.TestCase):
         self.assertEqual(f3, True)
         self.assertEqual(f4, True)
 
-        xs1 = self.design2h.helices[1].crossover_addresses(allow_intrahelix=True)
+        xs1 = self.design2h.helices[1].crossover_addresses(self.design2h.helices, allow_intrahelix=True)
         self.assertEqual(len(xs1), 5)
         h0, o0, f0 = xs1[0]
         h1, o1, f1 = xs1[1]