implemented logic for relaxing helix rolls based on crossover locatio…

…ns (not in UI yet)

lib/src/state/design.dart

@@ -2353,6 +2353,59 @@ abstract class Design with UnusedFields implements Built<Design, DesignBuilder>,
     }
     return strand_idx;
   }
+
+  /// maps helix indices to list of addresses of crossovers on that helix (helix_idx of Address is idx
+  /// of *other* helix, so the Addresses are not interpretable as Addresses in the Design, since the
+  /// other two parts (offset and forward) refer to the current helix (the key in the returned Map)
+  @memoized
+  BuiltMap<int, BuiltList<Address>> get helix_to_crossover_addresses {
+    Map<int, List<Address>> ret = {for (int helix_idx in this.helix_idxs) helix_idx: []};
+    for (int helix_idx in this.helix_idxs) {
+      var domains = this.domains_on_helix(helix_idx);
+      for (Domain domain in domains) {
+        var strand = this.substrand_to_strand[domain];
+        var domains_on_strand = strand.domains;
+        var num_domains = domains_on_strand.length;
+        var domain_idx = domains_on_strand.indexOf(domain);
+
+        // if not first domain, then there is a crossover to the previous domain
+        if (domain_idx > 0) {
+          var offset = domain.offset_5p;
+          var other_domain = domains_on_strand[domain_idx - 1];
+          var other_helix_idx = other_domain.helix;
+          ret[helix_idx].add(Address(helix_idx: other_helix_idx, offset: offset, forward: domain.forward));
+        }
+
+        // if not last domain, then there is a crossover to the next domain
+        if (domain_idx < num_domains - 1) {
+          var offset = domain.offset_3p;
+          var other_domain = domains_on_strand[domain_idx + 1];
+          var other_helix_idx = other_domain.helix;
+          ret[helix_idx].add(Address(helix_idx: other_helix_idx, offset: offset, forward: domain.forward));
+        }
+      }
+    }
+
+    // convert to built map of built lists
+    Map<int, BuiltList<Address>> ret_built_lists = {
+      for (int helix_idx in this.helix_idxs) helix_idx: ret[helix_idx].build()
+    };
+    return ret_built_lists.build();
+  }
+
+  /// returns design with all helix rolls relaxed (based on crossover locations)
+  Design relax_helix_rolls() {
+    Map<int, Helix> helices_relaxed = this.helices.toMap();
+
+    for (var helix_idx in helices_relaxed.keys) {
+      var helix = helices_relaxed[helix_idx];
+      var crossover_addresses = this.helix_to_crossover_addresses[helix_idx];
+      helix = helix.relax_roll(this.helices, crossover_addresses);
+      helices_relaxed[helix_idx] = helix;
+    }
+
+    return this.rebuild((b) => b..helices.replace(helices_relaxed));
+  }
 }
 
 Map<String, HelixGroup> _calculate_groups_from_helix_builder(

lib/src/state/helix.dart

@@ -2,6 +2,7 @@ import 'dart:math';
 
 import 'package:built_value/serializer.dart';
 import 'package:built_collection/built_collection.dart';
+import 'package:tuple/tuple.dart';
 import '../state/design.dart';
 import '../state/position3d.dart';
 import '../state/unused_fields.dart';
@@ -395,15 +396,13 @@ abstract class Helix with BuiltJsonSerializable, UnusedFields implements Built<H
     return ticks.build();
   }
 
-  double backbone_angle_at_offset(int offset, bool forward, Geometry geometry) {
+  double backbone_angle_at_offset(int offset, bool forward) {
     // Computes the backbone angle at *offset* for the strand in the direction given by *forward*.
     //
     // :param offset:
     //     offset on this helix
     // :param forward:
     //     whether to compute angle for the forward or reverse strand
-    // :param geometry:
-    //     :any:`Geometry` parameters to determine bases per turn
     // :return:
     //     backbone angle at *offset* for the strand in the direction given by *forward*.
     var degrees_per_base = 360 / geometry.bases_per_turn;
@@ -414,4 +413,23 @@ abstract class Helix with BuiltJsonSerializable, UnusedFields implements Built<H
     angle %= 360;
     return angle;
   }
+
+  Helix relax_roll(BuiltMap<int, Helix> helices, BuiltList<Address> crossover_addresses) {
+    var roll = compute_relaxed_roll(helices, crossover_addresses);
+    var new_helix = this.rebuild((b) => b..roll = roll);
+    return new_helix;
+  }
+
+  double compute_relaxed_roll(BuiltMap<int, Helix> helices, BuiltList<Address> crossover_addresses) {
+    List<Tuple2<double, double>> angles = [];
+    for (var address in crossover_addresses) {
+      var other_helix = helices[address.helix_idx];
+      var angle_of_other_helix = util.angle_from_helix_to_helix(this, other_helix);
+      var crossover_angle = this.backbone_angle_at_offset(address.offset, address.forward);
+      var relative_angle = Tuple2<double, double>(crossover_angle, angle_of_other_helix);
+      angles.add(relative_angle);
+    }
+    var angle = util.minimum_strain_angle(angles);
+    return angle;
+  }
 }

lib/src/util.dart

@@ -1793,3 +1793,138 @@ Map<int, int> invert_helices_view_order(Iterable<int> helices_view_order) {
   }
   return view_order_inverse;
 }
+
+double degrees(double rad) => rad * 180 / pi;
+
+double radians(double deg) => deg * pi / 180;
+
+// 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.
+double angle_from_helix_to_helix(Helix helix, Helix other_helix) {
+  var p1 = helix.position;
+  var p2 = other_helix.position;
+
+  // negate y_diff because y increases going down in the main view
+  var y_diff = -(p2.y - p1.y);
+  var x_diff = p2.x - p1.x;
+
+  var angle = degrees(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;
+}
+
+// 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.
+double minimum_strain_angle(List<Tuple2<double, double>> relative_angles) {
+  var adjusted_angles = [for (var angle in relative_angles) angle.item1 - angle.item2];
+  var ave_angle = average_angle(adjusted_angles);
+  var min_strain_angle = -ave_angle;
+  min_strain_angle %= 360;
+  return min_strain_angle;
+}
+
+// :param x: angle in degrees
+// :param y: angle in degrees
+// :return: signed difference between angles `x` and `y`, in degrees, in range [-180, 180]
+double angle_distance(double x, double y) {
+  var a = (x - y) % 360;
+  var b = (y - x) % 360;
+  var diff = a < b ? -a : b;
+  return diff;
+}
+
+// :param angles: list of angles in degrees
+// :param angle: angle in degrees
+// :return: sum of squared distances from each angle in `angles` to `angle`
+double sum_squared_angle_distances(List<double> angles, double angle) {
+  double sum = 0.0;
+  for (var a in angles) {
+    var dist = angle_distance(angle, a);
+    sum += dist * dist;
+  }
+  return sum;
+}
+
+// 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.
+double average_angle(List<double> angles) {
+  int num_angles = angles.length;
+  double mean_angle = 0.0;
+  if (num_angles > 0) {
+    mean_angle = angles.reduce((a, b) => a + b) / num_angles;
+  }
+
+  double min_dist = double.infinity;
+  double optimal_angle = 0;
+  for (int n = 0; n < num_angles; n++) {
+    var candidate_angle = mean_angle + 360.0 * n / num_angles;
+    var 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 ((360.0 - optimal_angle).abs() < 0.000000001) {
+    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;
+}
+
+double round(double x, int precision) {
+  var x_big = x * pow(10, precision);
+  int x_big_int = x_big.round();
+  x = x_big_int / pow(10, precision);
+  return x;
+}

lib/src/view/design_side_helix.dart

@@ -67,10 +67,8 @@ class DesignSideHelixComponent extends UiComponent2<DesignSideHelixProps> with P
       grid_position_str = position_str.replaceAll(' ', '');
     }
 
-    var forward_angle =
-        props.helix.backbone_angle_at_offset(props.slice_bar_offset, true, props.helix.geometry);
-    var reverse_angle =
-        props.helix.backbone_angle_at_offset(props.slice_bar_offset, false, props.helix.geometry);
+    var forward_angle = props.helix.backbone_angle_at_offset(props.slice_bar_offset, true);
+    var reverse_angle = props.helix.backbone_angle_at_offset(props.slice_bar_offset, false);
 
     var tooltip = '''\
 position:  ${position_str}