Merge pull request #901 from UC-Davis-molecular-computing/dev

Dev

lib/src/actions/actions.dart

@@ -1183,6 +1183,27 @@ abstract class HelixRollSetAtOther
   }
 }
 
+/////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// relax helix rolls
+
+abstract class RelaxHelixRolls
+    with BuiltJsonSerializable, UndoableAction
+    implements Built<RelaxHelixRolls, RelaxHelixRollsBuilder> {
+  bool get only_selected;
+
+  /************************ begin BuiltValue boilerplate ************************/
+  factory RelaxHelixRolls({bool only_selected}) = _$RelaxHelixRolls._;
+
+  RelaxHelixRolls._();
+
+  static Serializer<RelaxHelixRolls> get serializer => _$relaxHelixRollsSerializer;
+
+  @override
+  String short_description() {
+    return "set helix rolls to unstrain crossovers";
+  }
+}
+
 /////////////////////////////////////////////////////////////////////////////////////////////////////////////
 // Error message
 

lib/src/constants.dart

@@ -8,12 +8,12 @@ import 'state/grid.dart';
 
 // WARNING: Do not modify line below, except for the version string
 //          (and also add new version string to scadnano_versions_to_link).
-const String CURRENT_VERSION = "0.18.3";
+const String CURRENT_VERSION = "0.18.4";
 const String INITIAL_VERSION = "0.1.0";
 
 // scadnano versions that we deploy so that older versions can be used.
 final scadnano_older_versions_to_link = [
-  "0.18.2",
+  "0.18.3",
   "0.17.14",
   // "0.17.13",
   // "0.17.12",
@@ -91,7 +91,7 @@ final default_geometry = Geometry();
 const scadnano_css_stylesheet_name_no_ext = r'scadnano-styles';
 const scadnano_css_stylesheet_name = '${scadnano_css_stylesheet_name_no_ext}.css';
 
-const NUM_DIGITS_PRECISION_POSITION_DISPLAYED = 2;
+const NUM_DIGITS_PRECISION_POSITION_DISPLAYED = 1;
 
 /// DISTANCE_BETWEEN_HELICES_SVG is set to (BASE_WIDTH_SVG * 2.5/0.332) based on the following calculation,
 /// to attempt to make the DNA appear to scale in 2D drawings:

lib/src/middleware/edit_select_mode_change.dart

@@ -17,8 +17,7 @@ const selectable_css_style_non_domain_or_end = {
 };
 
 const selectable_css_style_domain = {
-  'stroke': 'hotpink',
-  'stroke-width': '5pt', // makes thicker when selected so easier to see
+  'filter': 'url("#shadow")',
 };
 
 const selectable_css_style_end = {

lib/src/reducers/helices_reducer.dart

@@ -32,6 +32,7 @@ Reducer<BuiltMap<int, Helix>> helices_local_reducer = combineReducers([
 ]);
 
 GlobalReducer<BuiltMap<int, Helix>, AppState> helices_global_reducer = combineGlobalReducers([
+  TypedGlobalReducer<BuiltMap<int, Helix>, AppState, actions.RelaxHelixRolls>(relax_helix_rolls_reducer),
   TypedGlobalReducer<BuiltMap<int, Helix>, AppState, actions.GroupChange>(helix_group_change_reducer),
   TypedGlobalReducer<BuiltMap<int, Helix>, AppState, actions.GridChange>(helix_grid_change_reducer),
   TypedGlobalReducer<BuiltMap<int, Helix>, AppState, actions.HelixGridPositionSet>(
@@ -478,6 +479,22 @@ BuiltMap<int, Helix> helix_grid_change_reducer(
   return new_helices.build();
 }
 
+BuiltMap<int, Helix> relax_helix_rolls_reducer(
+    BuiltMap<int, Helix> helices, AppState state, actions.RelaxHelixRolls action) {
+  var helix_idxs_to_relax =
+      action.only_selected ? state.ui_state.side_selected_helix_idxs : state.design.helix_idxs;
+
+  var new_helices_map = helices.toMap();
+  for (var helix_idx in helix_idxs_to_relax) {
+    var helix = new_helices_map[helix_idx];
+    var crossover_addresses = state.design.helix_to_crossover_addresses[helix_idx];
+    var helix_relaxed = helix.relax_roll(helices, crossover_addresses);
+    new_helices_map[helix_idx] = helix_relaxed;
+  }
+
+  return new_helices_map.build();
+}
+
 BuiltMap<int, Helix> helix_group_change_reducer(
     BuiltMap<int, Helix> helices, AppState state, actions.GroupChange action) {
   return helices.map_values((idx, helix) =>

lib/src/reducers/load_dna_file_reducer.dart

@@ -1,4 +1,5 @@
 import 'dart:convert';
+import 'dart:html';
 
 import 'package:built_collection/built_collection.dart';
 import 'package:scadnano/src/dna_file_type.dart';
@@ -38,11 +39,18 @@ The design has the following problem:
 ${error.cause}
 ${util.stack_trace_message_bug_report(stack_trace)}''';
   } catch (error, stack_trace) {
+    window.alert(
+        'I was unable to process that file. Only scadnano .sc files are supported for opening via the menu File-->Open or dragging onto the browser. If you are trying to import a cadnano file (ending. in .json), use the menu option File-->Import cadnano v2. Here is the full error message: ');
     error_message = "I encountered an error while reading the file ${action.filename}:"
         '\n\n$hline'
         '\n* error type:    ${error.runtimeType}'
         '\n* error message: ${error.toString()}'
         '\n$hline'
+        '\nIf the file is imported from cadnano, use the option of importing "cadnano v2".'
+        '\nTo do this, go to the menu options, then perform Select File --> Import "cadnano v2".'
+        '\nImporting a .json file is not allowed when selecting the "Open" option. Use a file with the .sc extension.'
+        '\nIt is the same thing for dragging a file into the browser. A .json file will also cause an error, even if it is dragged in.'
+        '\nEssentially, for .sc files, you can use the open or drag and drop feature. For .json files, you have to import "cadnano v2".'
         '\n\nThat file\'s contents are printed below.'
         '${util.stack_trace_message_bug_report(stack_trace)}'
         '\n\nThe file ${action.filename} has this content:\n\n${action.content}';

lib/src/serializers.dart

@@ -61,6 +61,7 @@ import 'state/domain_name_mismatch.dart';
 part 'serializers.g.dart';
 
 @SerializersFor([
+  RelaxHelixRolls,
   CopySelectedStandsToClipboardImage,
   LoadingDialogHide,
   LoadingDialogShow,

lib/src/state/design.dart

@@ -2353,6 +2353,61 @@ 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];
+      if (crossover_addresses.isNotEmpty) {
+        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/dialog.dart

@@ -149,6 +149,7 @@ abstract class Dialog with BuiltJsonSerializable implements Built<Dialog, Dialog
   @BuiltValueField(serialize: false, compare: false)
   ProcessCallback get process_saved_response;
 
+  // TODO: document this
   bool get use_saved_response;
 
   BuiltList<DialogItem> get items;

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';
@@ -394,4 +395,41 @@ abstract class Helix with BuiltJsonSerializable, UnusedFields implements Built<H
     }
     return ticks.build();
   }
+
+  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
+    // :return:
+    //     backbone angle at *offset* for the strand in the direction given by *forward*.
+    var degrees_per_base = 360 / geometry.bases_per_turn;
+    var angle = this.roll + offset * degrees_per_base;
+    if (!forward) {
+      angle += geometry.minor_groove_angle;
+    }
+    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;
+}