Merge pull request #264 from UC-Davis-molecular-computing/257-automat…

…ically-set-helix-rolls-based-on-crossover-locations-relax-the-rolls

257 automatically set helix rolls based on crossover locations relax the rolls

examples/output_designs/relax_helix_rolls.sc

@@ -0,0 +1,40 @@
+{
+  "version": "0.18.1",
+  "grid": "square",
+  "helices": [
+    {
+      "max_offset": 60,
+      "grid_position": [0, 0],
+      "roll": 40.71428571400003,
+      "major_ticks": [0, 5, 13]
+    },
+    {
+      "max_offset": 60,
+      "grid_position": [0, 1],
+      "roll": 72.85714285699999,
+      "major_ticks": [0, 5, 13]
+    },
+    {
+      "max_offset": 60,
+      "grid_position": [1, 0],
+      "roll": 68.57142857100001,
+      "major_ticks": [0, 5, 13]
+    }
+  ],
+  "strands": [
+    {
+      "color": "#f74308",
+      "domains": [
+        {"helix": 0, "forward": true, "start": 0, "end": 5},
+        {"helix": 1, "forward": false, "start": 0, "end": 5}
+      ]
+    },
+    {
+      "color": "#57bb00",
+      "domains": [
+        {"helix": 0, "forward": true, "start": 5, "end": 13},
+        {"helix": 2, "forward": false, "start": 5, "end": 13}
+      ]
+    }
+  ]
+}

examples/relax_helix_rolls.py

@@ -0,0 +1,78 @@
+import scadnano as sc
+import modifications as mod
+import dataclasses
+
+
+def create_design() -> sc.Design:
+    # '''
+    #   0123456789012345678901234567890123456789
+    # 0 [---+[--------+[----------+
+    #       |         |           |
+    # 1 [---+<--------+<----------+
+    #
+    # angle (fraction of 360)
+    #      5/10.5
+    #              (15-10.5)/10.5 = 4.5/10.5
+    #                        (27-21)/10.5 = 6/10.5
+    # '''
+    # design2h = sc.Design(helices=[sc.Helix(max_offset=50) for _ in range(2)], grid=sc.square)
+    # # helix 0 forward
+    # design2h.draw_strand(0, 0).move(5).cross(1).move(-5)
+    # design2h.draw_strand(0, 5).move(10).cross(1).move(-10)
+    # design2h.draw_strand(0, 15).move(12).cross(1).move(-12)
+    #
+    # for helix in design2h.helices.values():
+    #     helix.major_ticks = [0, 5, 15, 27]
+    #
+    # design2h.relax_helix_rolls()
+
+    '''
+      0         1         2         3         4         5         6
+      012345678901234567890123456789012345678901234567890123456789
+    0 [---+[--------+[----------+[------+[--------+[--------+
+          |         |           |       |         |         |
+    1 [---+<--------+<----------+       |         |         |
+                                        |         |         |
+    2                            <------+<--------+<--------+
+
+    angle (fraction of 360)
+         5/10.5
+                 (15-10.5)/10.5 = 4.5/10.5
+                           (27-21)/10.5 = 6/10.5
+    '''
+    helices = [sc.Helix(max_offset=60) for _ in range(3)]
+    helices[2].grid_position = (1, 0)
+    design3h = sc.Design(helices=helices, grid=sc.square)
+
+    # helix 0 forward
+    design3h.draw_strand(0, 0).move(5).cross(1).move(-5)
+    design3h.draw_strand(0, 5).move(10).cross(1).move(-10)
+    design3h.draw_strand(0, 15).move(12).cross(1).move(-12)
+    design3h.draw_strand(0, 27).move(7).cross(2).move(-7)
+    design3h.draw_strand(0, 34).move(10).cross(2).move(-10)
+    design3h.draw_strand(0, 44).move(10).cross(2).move(-10)
+
+    for helix in design3h.helices.values():
+        helix.major_ticks = [0, 5, 15, 27, 34, 44, 54]
+
+    design3h.relax_helix_rolls()
+
+    # return design3h
+
+    helices = [sc.Helix(max_offset=60) for _ in range(3)]
+    helices[2].grid_position = (1, 0)
+    design3h2 = sc.Design(helices=helices, grid=sc.square)
+    design3h2.draw_strand(0, 0).move(5).cross(1).move(-5)
+    design3h2.draw_strand(0, 5).move(8).cross(2).move(-8)
+
+    for helix in design3h2.helices.values():
+        helix.major_ticks = [0, 5, 13]
+
+    design3h2.relax_helix_rolls()
+
+    return design3h2
+
+
+if __name__ == '__main__':
+    d = create_design()
+    d.write_scadnano_file(directory='output_designs')

scadnano/scadnano.py

@@ -62,6 +62,7 @@
 import enum
 import itertools
 import re
+import math
 from builtins import ValueError
 from dataclasses import dataclass, field, InitVar, replace
 from typing import Iterator, Tuple, List, Sequence, Iterable, Set, Dict, Union, Optional, Type, cast, Any, \
@@ -335,9 +336,6 @@ class Grid(str, enum.Enum):
     Represents default patterns for laying out helices in the side view.
     Each :any:`Grid` except :py:data:`Grid.none` has an interpretation of a "grid position",
     which is a 2D integer coordinate (`h`, `v`).
-    (scadnano also allows a 3rd coordinate (`h`, `v`, `b`) specifying a "base offset" at which to position
-    the start of the :any:`Helix`, which is not relevant for the side view but will eventually be
-    supported to adjust the main view.)
     """
 
     square = "square"
@@ -1202,7 +1200,7 @@ def from_json(json_map: Dict[str, Any]) -> Position3D:
             z = json_map[position_z_key]
         return Position3D(x=x, y=y, z=z)
 
-    def __add__(self, other: 'Position3D') -> 'Position3D':
+    def __add__(self, other: Position3D) -> Position3D:
         """
         :param other:
             other position to add to this one
@@ -1211,7 +1209,7 @@ def __add__(self, other: 'Position3D') -> 'Position3D':
         """
         return Position3D(self.x + other.x, self.y + other.y, self.z + other.z)
 
-    def clone(self) -> 'Position3D':
+    def clone(self) -> Position3D:
         """
         :return:
             copy of this :any:`Position3D`
@@ -1699,6 +1697,26 @@ def calculate_major_ticks(self, grid: Grid) -> List[int]:
             distance = default_major_tick_distance(grid)
             return list(range(self.major_tick_start, self.max_offset + 1, distance))
 
+    def calculate_position(self, grid: Grid, geometry: Optional[Geometry] = None) -> Position3D:
+        """
+        :param grid:
+            :any:`Grid` of this :any:`Helix` used to interpret the field :data:`Helix.grid_position`.
+            Must be None if :data:`Helix.grid_position` is None.
+        :param geometry:
+            :any:`Geometry` parameters to determine distance between helices in a grid.
+            Must be None if :data:`Helix.grid_position` is None.
+        :return:
+            Position of this :any:`Helix` in 3D space, based on its :data:`Helix.grid_position`
+            if it is not None, or its :data:`Helix.position` otherwise.
+        """
+        if self.grid_position is None:
+            assert grid == Grid.none
+            return self.position
+        else:
+            assert grid is not None
+            assert geometry is not None
+            return grid_position_to_position(self.grid_position, grid, geometry)
+
     @property
     def domains(self) -> List[Domain]:
         """
@@ -1742,8 +1760,184 @@ def backbone_angle_at_offset(self, offset: int, forward: bool, geometry: Geometr
         angle = self.roll + offset * degrees_per_base
         if not forward:
             angle += geometry.minor_groove_angle
+        angle %= 360
         return angle
 
+    def crossovers(self) -> List[Tuple[int, int, bool]]:
+        """
+        :return:
+            list of triples (`offset`, `helix_idx`, `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]] = []
+        for domain in self.domains:
+            strand = domain.strand()
+            domains = strand.bound_domains()
+            num_domains = len(domains)
+            domain_idx = 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[domain_idx - 1]
+                other_helix_idx = other_domain.helix
+                crossovers.append((offset, other_helix_idx, 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[domain_idx + 1]
+                other_helix_idx = other_domain.helix
+                crossovers.append((offset, other_helix_idx, domain.forward))
+
+        return crossovers
+
+    def relax_roll(self, helices: Dict[int, Helix], grid: Grid, geometry: Geometry) -> None:
+        """
+        Like :meth:`Design.relax_helix_rolls`, but only for this :any:`Helix`.
+        """
+        angle = self.compute_relaxed_roll(helices, grid, geometry)
+        self.roll = angle
+
+    def compute_relaxed_roll(self, helices: Dict[int, Helix], grid: Grid, geometry: Geometry) -> float:
+        """
+        Like :meth:`Helix.relax_roll`, but just returns the new roll without altering this :any:`Helix`,
+        rather than changing the field :data:`Helix.roll`.
+        """
+        angles = []
+        for offset, helix_idx, forward in self.crossovers():
+            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)
+            relative_angle = (crossover_angle, angle_of_other_helix)
+            angles.append(relative_angle)
+        angle = minimum_strain_angle(angles)
+        return angle
+
+
+def angle_from_helix_to_helix(helix: Helix, other_helix: Helix,
+                              grid: Optional[Grid] = None, geometry: Optional[Geometry] = None) -> float:
+    """
+    Computes angle between `helix` and `other_helix` in degrees.
+
+    :param helix:
+        first helix
+    :param other_helix:
+        second helix
+    :param grid:
+        :any:`Grid` to use when calculating Helix positions
+    :param geometry:
+        :any:`Geometry` to use when calculating Helix positions
+    :return:
+        angle between `helix` and `other_helix` in degrees.
+    """
+    p1 = helix.calculate_position(grid, geometry)
+    p2 = other_helix.calculate_position(grid, geometry)
+
+    # negate y_diff because y increases going down in the main view
+    y_diff = -(p2.y - p1.y)
+    x_diff = p2.x - p1.x
+
+    angle = math.degrees(math.atan2(y_diff, x_diff))
+
+    # negate angle because we rotate clockwise
+    angle = -angle
+
+    # subtract 90 since we define 0 angle to be up instead of right
+    angle += 90
+
+    # normalize to be in range [0, 360)
+    angle %= 360
+
+    return angle
+
+
+def minimum_strain_angle(relative_angles: List[Tuple[float, float]]) -> float:
+    r"""
+    Computes the angle that minimizes the "strain" of all relative angles in the given list.
+
+    A "relative angle" is a pair :math:`(\theta, \mu)`. The strain is set to 0 by setting
+    :math:`\theta = \mu`; more generally the strain is :math:`(\theta-\mu)^2`, where :math:`\theta-\mu`
+    is the "angular difference" (e.g., 10-350 is 20 since 350 is also -10 mod 360).
+
+    The constraint is that in the list
+    [:math:`(\theta_1, \mu_1)`, :math:`(\theta_2, \mu_2)`, ..., :math:`(\theta_n, \mu_n)`],
+    we can rotate all angles :math:`\theta_i` by the same amount :math:`\theta`.
+    So this calculates the angle :math:`\theta` that minimizes
+    :math:`\sum_i [(\theta + \theta_i) - \mu_i]^2`
+
+    :param relative_angles:
+        List of :math:`(\theta_i, \mu_i)` pairs, where :math:`\theta_i = \mu_i` means 0 strain, and angles
+        are in units of degrees.
+    :return:
+        angle :math:`\theta` by which to rotate all angles :math:`\theta_i`
+        (but not changing any "zero angle" :math:`\mu_i`)
+        such that :math:`\sum_i [(\theta + \theta_i) - \mu_i]^2` is minimized.
+    """
+    adjusted_angles = [angle - zero_angle for angle, zero_angle in relative_angles]
+    ave_angle = average_angle(adjusted_angles)
+    min_strain_angle = -ave_angle
+    min_strain_angle %= 360
+    return min_strain_angle
+
+
+def angle_distance(x: float, y: float) -> float:
+    """
+    :param x: angle in degrees
+    :param y: angle in degrees
+    :return: signed difference between angles `x` and `y`, in degrees, in range [-180, 180]
+    """
+    a = (x - y) % 360
+    b = (y - x) % 360
+    diff = -a if a < b else b
+    return diff
+
+
+def sum_squared_angle_distances(angles: List[float], angle: float) -> float:
+    """
+    :param angles: list of angles in degrees
+    :param angle: angle in degrees
+    :return: sum of squared distances from each angle in `angles` to `angle`
+    """
+    return sum(angle_distance(angle, a) ** 2 for a in angles)
+
+
+def average_angle(angles: List[float]) -> float:
+    """
+    Calculate the "circular mean" of the angles in `angles`. Note this coincides with the arithemtic mean
+    for certain lists of angles, e.g., [0, 10, 50], in a way that the circular mean calculated via
+    interpreting angles as unit vectors (https://en.wikipedia.org/wiki/Circular_mean) does not.
+
+    This algorithm is due to Julian Panetta. (https://julianpanetta.com/)
+
+    :param angles:
+        List of angles in degrees.
+    :return:
+        average angle of the list of angles, normalized to be between 0 and 360.
+    """
+    num_angles = len(angles)
+    mean_angle = sum(angles) / num_angles
+    min_dist = float('inf')
+    optimal_angle = 0
+    for n in range(num_angles):
+        candidate_angle = mean_angle + 360.0 * n / num_angles
+        candidate_dist = sum_squared_angle_distances(angles, candidate_angle)
+        if min_dist > candidate_dist:
+            min_dist = candidate_dist
+            optimal_angle = candidate_angle
+
+    optimal_angle %= 360.0
+
+    # taking mod 360 sometimes results in 360.0 instead of 0.0. This is hacky but fixes it.
+    if abs(360.0 - optimal_angle) < 10 ** (-8):
+        optimal_angle = 0.0
+
+    # in case it's a nice round number, let's get rid of the floating-point error artifacts here
+    optimal_angle = round(optimal_angle, 9)
+
+    return optimal_angle
+
 
 def _is_close(x1: float, x2: float) -> bool:
     return abs(x1 - x2) < _floating_point_tolerance
@@ -8174,9 +8368,6 @@ def strand_with_name(self, name: str) -> Optional[Strand]:
                 return strand
         return None
 
-    def grid_of_helix(self, helix):
-        pass
-
     def add_helix(self, idx: int, helix: Helix) -> None:
         """
         Adds `helix` as a new :any:`Helix` with index `idx` to this Design.
@@ -8202,6 +8393,18 @@ def add_helix(self, idx: int, helix: Helix) -> None:
         self._set_helices_grid_positions_or_positions()
         self._assign_default_helices_view_orders_to_groups()
 
+    def relax_helix_rolls(self) -> None:
+        """
+        Sets all helix rolls to "relax" them based on their crossovers.
+
+        This calculates the "strain" of each crossover c as the absolute value d_c of the distance between
+        the angle to the helix to which it is connected and the angle of that crossover given the
+        current helix roll. It minimizes sum_c d_c^2, i.e., minimize the sum of the squares of the strains.
+        """
+        for helix in self.helices.values():
+            helix_group = self.groups[helix.group]
+            helix.relax_roll(self.helices, helix_group.grid, self.geometry)
+
 
 def _find_index_pair_same_object(elts: Union[List, Dict]) -> Optional[Tuple]:
     # return pair of indices representing same object in elts, or None if they do not exist