#![deny(missing_docs)]
//! This crate (codenano) allows one to design DNA molecules.
//!
//! The following creates a strand zig-zagging between three helices,
//! and outputs it to a file:
//!
//! ```
//! use codenano::*;
//! let mut design: Design<(), ()> = Design::new();
//! for i in 0..3 {
//!    design.add_grid_helix(i, 0);
//! }
//! // Let's first design a strand.
//! design.strand(0, 0).to(31)
//!   .cross(1).to(10)
//!   .cross(2).to(21);
//! // Now its reverse complement:
//! design.strand(2, 21).to(10)
//!    .cross(1).to(31)
//!    .cross(0).to(0);
//! design.write_to("my_design.json").unwrap()
//! ```
//!
//! This library is compatible with the `codenano-server` crate, which
//! watches changes made to a file and shows the result in 3D in a web
//! browser. Additionally, the client part (also called codenano) of the the `codenano-server`crate can run
//! finite element simulations and display the secondary structure.

#[macro_use]
extern crate serde_derive;
extern crate failure;
extern crate serde;

pub use failure::Error;
use std::borrow::Cow;
use std::collections::HashMap;
use std::f64::consts::PI;

/// 3D points and vectors
pub mod geometry;
use crate::geometry::*;

/// Extra tools for sequences (random completions, etc.)
pub mod sequences;

/// Conversion to Cadnano format.
pub mod cadnano;

/// The main type of this crate, describing a DNA design.
#[derive(Serialize, Deserialize)]
pub struct Design<StrandLabel, DomainLabel> {
    /// Version of this format.
    pub version: String,
    /// The vector of all helices used in this design. Helices have a
    /// position and an orientation in 3D.
    pub helices: Vec<Helix>,
    /// The vector of strands.
    pub strands: Vec<Strand<StrandLabel, DomainLabel>>,
    /// Parameters of DNA geometry. This can be skipped (in JSON), or
    /// set to `None` in Rust, in which case a default set of
    /// parameters from the literature is used.
    #[serde(skip_serializing_if = "Option::is_none", default)]
    pub parameters: Option<Parameters>,
}

fn is_false(x: &bool) -> bool {
    !*x
}

fn none<Label>() -> Option<Label> {
    None
}

/// A DNA strand.
#[derive(Debug, Default, Clone, Serialize, Deserialize)]
pub struct Strand<Label, DomainLabel> {
    /// The (ordered) vector of domains, where each domain is a
    /// directed interval of a helix.
    pub domains: Vec<Domain<DomainLabel>>,
    /// The sequence of this strand, if any. If the sequence is longer
    /// than specified by the domains, a prefix is assumed. Can be
    /// skipped in the serialisation.
    #[serde(skip_serializing_if = "Option::is_none", default)]
    pub sequence: Option<Cow<'static, str>>,
    /// Is this sequence cyclic? Can be skipped (and defaults to
    /// `false`) in the serialization.
    #[serde(skip_serializing_if = "is_false", default)]
    pub cyclic: bool,
    /// Colour of this strand. If skipped, a default colour will be
    /// chosen automatically.
    #[serde(skip_serializing_if = "Option::is_none", default)]
    pub color: Option<u32>,
    /// An optional label for the strand. Can be
    /// `serde_json::Value::Null`, and skipped in the serialisation.
    #[serde(skip_serializing_if = "Option::is_none", default = "none")]
    pub label: Option<Label>,
}

impl<Label, DomainLabel> Strand<Label, DomainLabel> {
    /// Find the nucl-th nucleotide of the strand. Output is (h, x, b) so that the nucl-th nucleotide
    /// of the strand is on helix h at position x with orientation b.
    pub fn find_nucl(&self, nucl: isize) -> (isize, isize, bool) {
        let mut current = 0;
        let mut i = 0;
        while current + &self.domains[i].length() <= nucl {
            current += &self.domains[i].length();
            i += 1;
        }
        let position = if self.domains[i].right {
            self.domains[i].start + nucl - current
        } else {
            self.domains[i].end - 1 - nucl + current
        };
        (self.domains[i].helix, position, self.domains[i].right)
    }
}



/// A domain, i.e. an interval of a helix.
#[derive(Debug, Default, Clone, Serialize, Deserialize)]
pub struct Domain<Label> {
    /// Index of the helix in the array of helices. Indices start at
    /// 0.
    pub helix: isize,
    /// Position of the leftmost base of this domain along the helix
    /// (this might be the first or last base of the domain, depending
    /// on the `orientation` parameter below).
    pub start: isize,
    /// Position of the first base after the rightmost base of the
    /// domain, along the helix. Domains must always be such that
    /// `domain.start < domain.end`.
    pub end: isize,
    /// If true, the "5' to 3'" direction of this domain runs in the
    /// same direction as the helix, i.e. "to the right" along the
    /// axis of the helix. Else, the 5' to 3' runs to the left along
    /// the axis.
    pub right: bool,
    /// An optional label that can be attached to strands.
    #[serde(skip_serializing_if = "Option::is_none", default = "none")]
    pub label: Option<Label>,
    /// In addition to the strand-level sequence, individual domains
    /// may have sequences too. The precedence has to be defined by
    /// the user of this library.
    pub sequence: Option<Cow<'static, str>>,
}

impl<Label> Domain<Label> {
    /// Iterate through the positions of this domain, in 5' to 3'
    /// order (meaning that the values produced by this iterator might
    /// be increasing or decreasing).
    pub fn iter(&self) -> DomainIter {
        DomainIter {
            start: self.start,
            end: self.end,
            right: self.right,
        }
    }
    /// Translate this domain. The first parameter is the translation
    /// along the helix, the second one is a translation across
    /// helices (probably most meaningful for a flat design).
    pub fn translate(self, dx: isize, dy: isize) -> Self {
        use std::convert::TryFrom;
        Domain {
            start: self.start + dx,
            end: self.end + dx,
            helix: usize::try_from(self.helix as isize + dy).unwrap() as isize,
            ..self
        }
    }
    /// Translate this domain along its helix.
    pub fn shift_x(self, dx: isize) -> Self {
        Domain {
            start: self.start + dx,
            end: self.end + dx,
            ..self
        }
    }
    /// Translate this domain to a different helix (probably most
    /// meaningful for a flat design).
    pub fn shift_y(self, dy: isize) -> Self {
        use std::convert::TryFrom;
        Domain {
            helix: usize::try_from(self.helix as isize + dy).unwrap() as isize,
            ..self
        }
    }

    /// Number of Nucleotides on the domain
    pub fn length(&self) -> isize {
        self.end - self.start
    }
}

/// An iterator over all positions of a domain.
pub struct DomainIter {
    start: isize,
    end: isize,
    right: bool,
}

impl Iterator for DomainIter {
    type Item = isize;
    fn next(&mut self) -> Option<Self::Item> {
        if self.start == self.end {
            None
        } else {
            if self.right {
                let s = self.start;
                self.start += 1;
                Some(s)
            } else {
                let s = self.end;
                self.end -= 1;
                Some(s - 1)
            }
        }
    }
}

/// The sequence of the M13 bacteriophage, used for example for DNA
/// Origami.
pub const M13_7249: &'static str = include_str!("m13");
/// A cut version of the M13, which can be obtained with enzymes from
/// the original genome (the `M13_7249` one).
pub const M13_594: &'static str = include_str!("m13_594");

/// DNA geometric parameters.
#[derive(Copy, Clone, Debug, Serialize, Deserialize)]
pub struct Parameters {
    /// Distance between two consecutive bases along the axis of a
    /// helix, in nanometers.
    pub z_step: f64,
    /// Radius of a helix, in nanometers.
    pub helix_radius: f64,
    /// Number of bases per turn in nanometers.
    pub bases_per_turn: f64,
    /// Minor groove angle. DNA helices have a "minor groove" and a
    /// "major groove", meaning that two paired nucleotides are not at
    /// opposite positions around a double helix (i.e. at an angle of
    /// 180°), but instead have a different angle.
    ///
    /// Strands are directed. The "normal" direction is called "5' to
    /// 3'" (named after parts of the nucleotides). This parameter is
    /// the small angle, which is clockwise from the normal strand to
    /// the reverse strand.
    pub groove_angle: f64,

    /// Gap between two neighbouring helices.
    pub inter_helix_gap: f64,
}

impl Default for Parameters {
    fn default() -> Self {
        Parameters {
            // z-step and helix radius from:
            //
            // Single-molecule portrait of DNA and RNA double helices,
            // J. Ricardo Arias-Gonzalez, Integrative Biology, Royal
            // Society of Chemistry, 2014, vol. 6, p.904
            z_step: 0.332,
            helix_radius: 1.,
            // bases per turn from Wu Rothemund (Nature Chemistry).
            bases_per_turn: 10.44,
            groove_angle: -24. * PI / 34.,
            // From Paul's paper.
            inter_helix_gap: 1.,
        }
    }
}

/// A DNA helix. All bases of all strands must be on a helix.
///
/// The three angles are illustrated in the following image, from [the NASA website](https://www.grc.nasa.gov/www/k-12/airplane/rotations.html):
///
/// ![Aircraft angles](https://www.grc.nasa.gov/www/k-12/airplane/Images/rotations.gif)
#[derive(Serialize, Deserialize)]
pub struct Helix {
    /// Position of the origin of the helix axis.
    pub origin: Point<f64>,

    /// Angle around the axis of the helix.
    pub roll: f64,

    /// Horizontal rotation.
    pub yaw: f64,

    /// Vertical rotation.
    pub pitch: f64,
}

impl Helix {
    /// Angle of base number `n` around this helix.
    pub fn theta(&self, n: isize, right: bool, cst: &Parameters) -> f64 {
        let shift = if right {
            cst.groove_angle
        } else {
            0.
        };
        n as f64 * 2. * PI / cst.bases_per_turn + shift + self.roll
    }

    /*
    // Return `R(y,theta)*point` where `R(u,theta)` is the matrix of the rotation of axis u and
    // angle theta and `*` is a matrix,vector multiplication
    fn ry(point: &[f64], theta: f64) -> [f64; 3] {
        [
            point[0] * theta.cos() + point[2] * theta.sin(),
            point[1],
            point[0] * -theta.sin() + point[2] * theta.cos(),
        ]
    }

    // Return `R(z,theta)*point` where `R(u,theta)` is the matrix of the rotation of axis u and
    // angle theta and `*` is a matrix,vector multiplication
    fn rz(point: &[f64], theta: f64) -> [f64; 3] {
        [
            point[0] * theta.cos() + point[1] * -theta.sin(),
            point[0] * theta.sin() + point[1] * theta.cos(),
            point[2],
        ]
    }
    */

    /// 3D position of a nucleotide on this helix. `n` is the position along the axis, and `right` is true iff the 5' to 3' direction of the strand containing that nucleotide runs in the same direction as the axis of the helix.
    pub fn space_pos(&self, p: &Parameters, n: isize, right: bool) -> [f64; 3] {
        let theta = self.theta(n, right, p);
        let mut ret = [
            n as f64 * p.z_step,
            theta.cos() * p.helix_radius,
            theta.sin() * p.helix_radius,
        ];
        
        let forward = [self.yaw.cos() * self.pitch.cos(),
                       self.pitch.sin(),
                       -self.yaw.sin() * self.pitch.cos()];
        let right = [self.yaw.sin(),
                     0.,
                     self.yaw.cos()];
        let up = [ right[1] * forward[2] - right[2] * forward[1],
                   right[2] * forward[0] - right[0] * forward[2],
                   right[0] * forward[1] - right[1] * forward[0]];

        ret = [ret[0] * forward[0] + ret[1] * up[0] + ret[2] * right[0],
         ret[0] * forward[1] + ret[1] * up[1] + ret[2] * right[1],
         ret[0] * forward[2] + ret[1] * up[2] + ret[2] * right[2],
        ];
        /*
        ret = Helix::ry(&ret, self.yaw);
        ret = Helix::rz(&ret, self.pitch);
        */
        ret[0] += self.origin.x;
        ret[1] += self.origin.y;
        ret[2] += self.origin.z;
        ret
    }

    /// 3D position of the projection of the nucleotide on its helix.
    /// `n` is the position along the axis.
    pub fn axis_pos(&self, p: &Parameters, n:isize) -> [f64; 3] {
        let mut ret = [
            n as f64 * p.z_step,
            0.,
            0.,
        ];
        
        let forward = [self.yaw.cos() * self.pitch.cos(),
                       self.pitch.sin(),
                       -self.yaw.sin() * self.pitch.cos()];
        let right = [self.yaw.sin(),
                     0.,
                     self.yaw.cos()];
        let up = [ right[1] * forward[2] - right[2] * forward[1],
                   right[2] * forward[0] - right[0] * forward[2],
                   right[0] * forward[1] - right[1] * forward[0]];

        ret = [ret[0] * forward[0] + ret[1] * up[0] + ret[2] * right[0],
         ret[0] * forward[1] + ret[1] * up[1] + ret[2] * right[1],
         ret[0] * forward[2] + ret[1] * up[2] + ret[2] * right[2],
        ];

        ret[0] += self.origin.x;
        ret[1] += self.origin.y;
        ret[2] += self.origin.z;
        ret


    }
}

/// Identity of a helix, useful for referring to it.
#[derive(Debug, Copy, Clone, PartialEq, Eq, PartialOrd, Ord, Hash)]
pub struct HelixId(pub isize);

impl From<isize> for HelixId {
    fn from(e: isize) -> Self {
        HelixId(e)
    }
}

impl From<i32> for HelixId {
    fn from(e: i32) -> Self {
        HelixId(e as isize)
    }
}

impl From<usize> for HelixId {
    fn from(e: usize) -> Self {
        HelixId(e as isize)
    }
}

/// Identity of a strand, useful for referring to it.
pub struct StrandId(pub usize);

/// This is just a communication type between the server and the
/// client, serialisable in JSON. No need to really expose it.
#[doc(hidden)]
#[derive(Serialize, Deserialize)]
pub struct Res {
    pub out: String,
    pub err: String,
}

/// A pointer to a strand being built, containing a current domain and
/// position.
pub struct StrandRef<'a, StrandLabel, DomainLabel> {
    builder: &'a mut Design<StrandLabel, DomainLabel>,
    strand_id: usize,
    x: Option<isize>,
    x0: Option<isize>,
}

impl<StrandLabel: serde::Serialize, DomainLabel: serde::Serialize>
    Design<StrandLabel, DomainLabel>
{
    /// Initiates a design.
    pub fn new() -> Self {
        Design {
            version: env!("CARGO_PKG_VERSION").to_string(),
            helices: Vec::new(),
            strands: Vec::new(),
            parameters: Some(Parameters::default()),
        }
    }

    /// Creates a flat helix (in the `z = 0` plane).
    pub fn square_grid_helix(&self, h: isize, v: isize) -> Helix {
        Helix {
            origin: Point {
                x: 0.,
                y: h as f64
                    * (self.parameters.unwrap().helix_radius * 2.
                        + self.parameters.unwrap().inter_helix_gap),
                z: v as f64
                    * (self.parameters.unwrap().helix_radius * 2.
                        + self.parameters.unwrap().inter_helix_gap),
            },
            roll: 0.,
            pitch: 0.,
            yaw: 0.,
        }
    }

    /// Creates a flat helix in the `z = 0` plane with adjusted roll.
    pub fn adjusted_grid_helix(&self, h: isize, v: isize) -> Helix {
        let delta_roll = self.parameters.unwrap().groove_angle - PI;
        Helix {
            origin: Point {
                x: 0.,
                y: h as f64
                    * (self.parameters.unwrap().helix_radius * 2.
                        + self.parameters.unwrap().inter_helix_gap),
                z: v as f64
                    * (self.parameters.unwrap().helix_radius * 2.
                        + self.parameters.unwrap().inter_helix_gap),
            },
            roll: ((h + v) as f64) * delta_roll,
            pitch: 0.,
            yaw: 0.,
        }
    }

    /// Initiates a design with the given geometric parameters.
    pub fn with_parameters(parameters: Parameters) -> Self {
        let mut b = Design::new();
        b.parameters = Some(parameters);
        b
    }

    /// Add a helix created by `Self::square_grid_helix`.
    pub fn add_grid_helix(&mut self, h: isize, v: isize) -> HelixId {
        self.helices.push(self.square_grid_helix(h, v));
        HelixId((self.helices.len() - 1) as isize)
    }

    /// Add a helix created by `Self::add_grid_helix`.
    pub fn add_adjusted_grid_helix(&mut self, h: isize, v: isize) -> HelixId {
        self.helices.push(self.adjusted_grid_helix(h, v));
        HelixId((self.helices.len() - 1) as isize)
    }


    /// Add an arbitrary helix. Parameters x,y,z are in nanometers and are the 3D coordinates
    /// of the centre of the helix. roll,pitch,yaw are in radians.  
    pub fn add_helix(
        &mut self,
        x: f64,
        y: f64,
        z: f64,
        roll: f64,
        pitch: f64,
        yaw: f64,
    ) -> HelixId {
        self.helices.push(Helix {
            origin: Point { x, y, z },
            roll,
            pitch,
            yaw,
        });
        HelixId((self.helices.len() - 1) as isize)
    }

    /// Start a strand, placing the first base at the position given
    /// by `helix` and `start`.
    pub fn strand<H: Into<HelixId>>(
        &mut self,
        helix: H,
        start: isize,
    ) -> StrandRef<StrandLabel, DomainLabel> {
        let strand_id = self.strands.len();
        self.strands.push(Strand {
            domains: vec![Domain {
                helix: helix.into().0,
                start,
                end: start,
                right: true,
                label: None,
                sequence: None,
            }],
            label: None,
            sequence: None,
            cyclic: false,
            color: None,
        });
        StrandRef {
            strand_id,
            builder: self,
            x: Some(start),
            x0: Some(start),
        }
    }

    /// Create a strand with no initial domain. Useful for adding
    /// custom domains, for example domains obtained by geometric
    /// operations from other domains.
    pub fn empty_strand(&mut self) -> StrandRef<StrandLabel, DomainLabel> {
        let strand_id = self.strands.len();
        self.strands.push(Strand {
            domains: Vec::new(),
            cyclic: false,
            color: None,
            sequence: None,
            label: None,
        });
        StrandRef {
            strand_id,
            builder: self,
            x: None,
            x0: None,
        }
    }

    /// Return a StrandRef that can build the strand `strand_id`
    pub fn get_strand_ref(&mut self, strand_id: usize) -> StrandRef<StrandLabel, DomainLabel> {
        let strand = &self.strands[strand_id];
        if strand.cyclic {
            panic!("Attempt to keep building on a cyclic strand");
        }
        let domain = &strand.domains[strand.domains.len() - 1];
        let x0 = if domain.right {
            domain.start
        } else {
            domain.end - 1
        };
        let x = if domain.right {
            domain.end - 1
        } else {
            domain.start
        };
        StrandRef {
            strand_id,
            builder: self,
            x: Some(x),
            x0: Some(x0)
        }
    }

    /// Write this design to the standard output.
    pub fn write(&self) -> Result<(), Error> {
        serde_json::to_writer(std::io::stdout(), self)?;
        Ok(())
    }

    /// Write this design to a file. If the `REMOTE` environment
    /// variable is set, write it to the standard output instead.
    pub fn write_to<P: AsRef<std::path::Path>>(&self, p: P) -> Result<(), Error> {
        if std::env::var("REMOTE").is_ok() {
            serde_json::to_writer_pretty(std::io::stdout(), &self)?
        } else {
            serde_json::to_writer_pretty(std::fs::File::create(p)?, &self)?
        }
        Ok(())
    }

    /// Tests whether all bases have different positions, and returns
    /// an offending strand else. The return format is `(helix,
    /// position, orientation, i, j)`, where `i` and `j` are the indices of
    /// two strands sharing this position.
    /// Bases that are the start/end of a cycle are not detected and do not cause their
    /// strand to be offending
    pub fn unique_positions(&self) -> Option<(isize, isize, bool, usize, usize)> {
        let mut positions = HashMap::new();
        for (i, s) in self.strands.iter().enumerate() {
            for dom in s.domains.iter() {
                for pos in dom.iter() {
                    if let Some(prev) = positions.insert((dom.helix, pos, dom.right), i) {
                        if s.cyclic {
                            // unwrap OK because we know s.domain[0] has at least one nucl
                            let j = positions.get(&(s.domains[0].helix, s.domains[0].iter().next().unwrap(),
                                s.domains[0].right));
                            if let Some(x) = j {
                                if *x == prev {
                                    continue;
                                }
                            }
                        }
                            return Some((dom.helix, pos, dom.right, i, prev));
                    }
                }
            }
        }
        None
    }


    /// Tests whether all bases have different positions, and warns the user otherwise.
    /// Write this design to a file. If the `REMOTE` environment
    /// variable is set, write it to the standard output instead.
    pub fn check_and_write_to<P: AsRef<std::path::Path>>(&self, p: P) -> Result<(), Error> {
        if let Some((helix, pos, right, i, prev)) = self.unique_positions() {
            let sens = if right { "Sens" } else { "AntiSens" };
            println!("WARNING: SOME STRANDS ARE CONFLICTING");
            println!("First conflict found: On helix {}, at position {}, on the {} strand",
                helix,
                pos,
                sens);
            println!("Strands number {} and {} are in conflict", i, prev);
        }
        self.write_to(p)
    }

    fn get_cross_overs(&self) -> Vec<(isize, isize, bool, isize, isize, bool)> {
        let mut ret = Vec::new();
        for s in &self.strands {
            for i in 0..(s.domains.len() - 1) {
                if s.domains[i].helix != s.domains[i + 1].helix {
                    let origin = if s.domains[i].right {
                        s.domains[i].end - 1
                    } else {
                        s.domains[i].start
                    };
                    let destination = if s.domains[i + 1].right {
                        s.domains[i + 1].start
                    } else {
                        s.domains[i + 1].end - 1
                    };
                    ret.push((s.domains[i].helix,
                              origin,
                              s.domains[i].right,
                              s.domains[i + 1].helix,
                              destination,
                              s.domains[i + 1].right));
                }
            }
        }
        ret
    }

    fn get_cross_overs_helices(&self, h1: isize, h2:isize) -> Vec<(isize, bool, isize, bool)> {
        let cross_overs = self.get_cross_overs();
        let mut ret = Vec::new();
        for (g1, x1, b1, g2, x2, b2) in cross_overs {
            if g1 == h1 && g2 == h2 {
                ret.push((x1, b1, x2, b2));
            } else if g2 == h1 && g1 == h2 {
                ret.push((x2, b2, x1, b1));
            }
        }
        ret
    }

    fn get_total_dist_helices(&self, h1: isize, h2:isize, theta: f64, phi: f64) -> f64 {
        let cross_overs = self.get_cross_overs_helices(h1, h2);
        let helix1 = Helix {
            roll: theta,
            ..self.helices[h1 as usize]
        };
        let helix2 = Helix {
            roll: phi,
            ..self.helices[h2 as usize]
        };

        let mut ret = 0.;
        //println!("cross_overs: {:?}", cross_overs);

        for (x1, b1, x2, b2) in cross_overs {
            let pos_1 = helix1.space_pos(&self.parameters.unwrap(), x1, b1);
            let pos_2 = helix2.space_pos(&self.parameters.unwrap(), x2, b2);
            ret += (pos_1[0] - pos_2[0]) * (pos_1[0] - pos_2[0]);
            ret += (pos_1[1] - pos_2[1]) * (pos_1[1] - pos_2[1]);
            ret += (pos_1[2] - pos_2[2]) * (pos_1[2] - pos_2[2]);
        }
        ret
    }

    /// Adjust the roll of two helices to minimize the distance at crossing overs
    pub fn adjust_helix(&mut self, h1: usize, h2: usize) {
        let mut theta = self.helices[h1].roll;
        let mut phi = self.helices[h2].roll;

        loop {
            let best  = self.get_total_dist_helices(h1 as isize, h2 as isize, theta, phi);
            //println!("{}, {}, {}", best, theta, phi);
            if self.get_total_dist_helices(h1 as isize, h2 as isize, theta - 0.1, phi - 0.1) < best {
                theta -= 0.1;
                phi -= 0.1;
            } else if self.get_total_dist_helices(h1 as isize, h2 as isize, theta + 0.1, phi - 0.1) < best {
                theta += 0.1;
                phi -= 0.1;
            } else if self.get_total_dist_helices(h1 as isize, h2 as isize, theta - 0.1, phi + 0.1) < best {
                theta -= 0.1;
                phi += 0.1;
            } else if self.get_total_dist_helices(h1 as isize, h2 as isize, theta + 0.1, phi + 0.1) < best {
                theta += 0.1;
                phi += 0.1;
            } else if self.get_total_dist_helices(h1 as isize, h2 as isize, theta, phi + 0.1) < best {
                phi += 0.1;
            } else if self.get_total_dist_helices(h1 as isize, h2 as isize, theta, phi - 0.1) < best {
                phi -= 0.1;
            } else if self.get_total_dist_helices(h1 as isize, h2 as isize, theta - 0.1, phi) < best {
                theta -= 0.1;
            } else if self.get_total_dist_helices(h1 as isize, h2 as isize, theta + 0.1, phi) < best {
                theta += 0.1;
            } else {
                self.helices[h1].roll = theta;
                self.helices[h2].roll = phi;
                return;
            }
        }
    }

}

// http://eleanormaclure.files.wordpress.com/2011/03/colour-coding.pdf
const KELLY: [u32; 19] = [
    // 0xF2F3F4, // White
    // 0x222222, // Black,
    0xF3C300, 0x875692, // 0xF38400, // Orange, too close to others
    0xA1CAF1, 0xBE0032, 0xC2B280, 0x848482, 0x008856, 0xE68FAC, 0x0067A5, 0xF99379, 0x604E97,
    0xF6A600, 0xB3446C, 0xDCD300, 0x882D17, 0x8DB600, 0x654522, 0xE25822, 0x2B3D26,
];

impl<StrandLabel, DomainLabel> Strand<StrandLabel, DomainLabel> {
    /// Provide a default color to the strand.
    pub fn default_color(&self) -> u32 {
        for domain in self.domains.iter() {
            let x1 = if domain.right {
                domain.end - 1
            } else {
                domain.start
            };
            let h = domain.helix as isize;
            let x = x1 + (x1 % 11) + 5 * h;
            let n = KELLY.len() as isize;
            return KELLY[(((x % n) + n) % n) as usize];
        }
        0
    }

    /// Set the color of this strand to be Kelly's nth color.
    pub fn set_kelly_color(&mut self, n: usize) {
        self.color = Some(KELLY[n % KELLY.len()]);
    }
}

impl<'a:'b, 'b, StrandLabel, DomainLabel> StrandRef<'a, StrandLabel, DomainLabel> {
    /// Retrieve the identity of the strand currently being built.
    pub fn id(&self) -> StrandId {
        StrandId(self.strand_id)
    }

    /// Return self.builder
    pub fn design(&mut self) -> &mut Design<StrandLabel, DomainLabel> {
        self.builder
    }

    /// Total (current) length of the strand being created. Adding
    /// more points will increase the length.
    pub fn len(&self) -> usize {
        self.domains()
            .iter()
            .map(|d| (d.end - d.start) as usize)
            .sum()
    }

    /// Choose a color for this strand.
    pub fn with_color(self, color: u32) -> Self {
        self.builder.strands[self.strand_id].color = Some(color);
        self
    }

    /// Choose a color for this strand.
    pub fn with_kelly_color(self, color: usize) -> Self {
        self.builder.strands[self.strand_id].color = Some(KELLY[color]);
        self
    }

    /// List the current domains of this strand (more domains can be
    /// added later).
    pub fn domains(&self) -> &[Domain<DomainLabel>] {
        &self.builder.strands[self.strand_id].domains
    }

    /// Gives a mutable pointer to the vector of domains.
    pub fn domains_mut(&mut self) -> &mut Vec<Domain<DomainLabel>> {
        &mut self.builder.strands[self.strand_id].domains
    }

    /// Push a domain onto this strand.
    pub fn push_domain(&mut self, domain: Domain<DomainLabel>) {
        self.x = Some(if domain.right {
            domain.end - 1
        } else {
            domain.start
        });
        self.builder.strands[self.strand_id].domains.push(domain);
    }

    /// Extend this strand with multiple domains.
    pub fn extend<I: Iterator<Item = Domain<DomainLabel>>>(&mut self, i: I) {
        for i in i {
            self.push_domain(i)
        }
    }

    /// Sets the sequence of this strand.
    pub fn with_sequence<I: Into<Cow<'static, str>>>(self, sequence: I) -> Self {
        self.builder.strands[self.strand_id].sequence = Some(sequence.into());
        self
    }

    /// Sets the sequence of this strand.
    pub fn with_domain_sequence<I: Into<Cow<'static, str>>>(self, sequence: I) -> Self {
        self.builder.strands[self.strand_id]
            .domains
            .last_mut()
            .unwrap()
            .sequence = Some(sequence.into());
        self
    }

    /// Move along the helix. This method can be called multiple times
    /// consecutively, but the direction of the domain is set only on
    /// the first call, subsequent calls only change the 3' end of the
    /// domain (which is the rightmost point if `self.right` is true,
    /// and the leftmost point else).
    pub fn to(mut self, x: isize) -> Self {
        let strand = self.builder.strands[self.strand_id]
            .domains
            .last_mut()
            .unwrap();
        self.x = Some(x);
        let x0 = self.x0.unwrap();
        if x > x0 {
            strand.right = true;
            strand.start = x0;
            strand.end = x + 1;
        } else {
            strand.right = false;
            strand.start = x;
            strand.end = x0 + 1;
        }
        self
    }

    /// Remove the last domain of the strand.
    pub fn pop(mut self) -> Self {
        self.builder.strands[self.strand_id].domains.pop();
        if let Some(last) = self.builder.strands[self.strand_id].domains.last() {
            self.x = Some(if last.right {
                last.end - 1
            } else {
                last.start
            })
        } else {
            self.x = None
        }
        self
    }

    /// Change the initial point of the first domain of this strand.
    pub fn rev_to(mut self, x: isize) -> Self {
        if self.builder.strands[self.strand_id].domains.is_empty() {
            return self;
        }
        let ref mut strand = self.builder.strands[self.strand_id].domains[0];
        self.x = Some(x);
        let x0 = self.x0.unwrap();
        if x > x0 {
            strand.right = true;
            strand.start = x0;
            strand.end = x + 1;
        } else {
            strand.right = false;
            strand.start = x;
            strand.end = x0 + 1;
        }
        self
    }

    /// Give a label to this strand.
    pub fn with_label<L: Into<StrandLabel>>(self, label: L) -> Self {
        self.builder.strands[self.strand_id].label = Some(label.into());
        self
    }

    /// Give a label to the most recently produced domain on this strand.
    pub fn with_domain_label<L: Into<DomainLabel>>(self, label: L) -> Self {
        self.builder.strands[self.strand_id]
            .domains
            .last_mut()
            .unwrap()
            .label = Some(label.into());
        self
    }

    /// Finish a domain, and start a new one on the same helix and at
    /// the same position.
    pub fn next_domain_to(mut self, to: isize) -> Self {
        let helix = self.builder.strands[self.strand_id]
            .domains
            .last()
            .unwrap()
            .helix;
        // If the previous domain is on the same helix.
        let (start, end, right) = if to < self.x.unwrap() {
            let n = (to, self.x.unwrap(), false);
            self.x = Some(to);
            n
        } else {
            let n = (self.x.unwrap() + 1, to + 1, true);
            self.x = Some(to);
            n
        };
        let domain = Domain {
            helix,
            start,
            end,
            right,
            label: None,
            sequence: None,
        };
        self.builder.strands[self.strand_id].domains.push(domain);
        self
    }

    /// Move the cursor to a different helix, keeping the same current
    /// position along the axis of the helix (most meaningful when the helices are in the same plane).
    pub fn cross<H: Into<HelixId>>(self, h: H) -> Self {
        let x = self.x.unwrap();
        self.cross_to(h, x)
    }

    /// Move the cursor to a different helix, at the specified
    /// position.
    pub fn cross_to<H: Into<HelixId>>(mut self, h: H, start: isize) -> Self {
        let helix = h.into().0;
        if helix >= self.builder.helices.len() as isize {
            panic!("Crossing to an undefined helix")
        }
        self.builder.strands[self.strand_id].domains.push(Domain {
            helix,
            start,
            end: start,
            right: true,
            label: None,
            sequence: None,
        });
        self.x0 = Some(start);
        self.x = Some(start);
        self
    }

    /// Compute the horizontal symmetry of this strand with respect to
    /// x-coordinate `x`.
    pub fn hflip(self, x: isize) -> Self {
        for dom in self.builder.strands[self.strand_id].domains.iter_mut() {
            let start = dom.start;
            dom.start = x - dom.end;
            dom.end = x - start;
        }
        self.builder.strands[self.strand_id].domains.reverse();
        self
    }

    /// Compute the vertical symmetry of this strand with respect to
    /// y-coordinate `y`. Most meaningful in flat designs.
    pub fn vflip<H: Into<HelixId>>(self, y: H) -> Self {
        let y = y.into();
        for dom in self.builder.strands[self.strand_id].domains.iter_mut() {
            assert!(y.0 >= dom.helix);
            dom.helix = 2 * y.0 - dom.helix
        }
        self
    }

    /// Translate this strand (all domains) to a different
    /// position. `x` is the horizontal translation (along the axis of
    /// the helices), and `y` moves between helices. Most meaningful
    /// in a flat, on-grid design.
    pub fn translate(self, x: isize, y: isize) -> Self {
        use std::convert::TryFrom;
        for dom in self.builder.strands[self.strand_id].domains.iter_mut() {
            dom.helix = usize::try_from(dom.helix as isize + y).unwrap() as isize;
            dom.start += x;
            dom.end += x;
        }
        self
    }

    /// Ends a cyclic strands.
    /// Set `self.cycle` to `true` and link the end and the start of the strand.
    pub fn cycle(self) -> () {
        let helix = self.builder.strands[self.strand_id].domains[0].helix;
        let start = self.builder.strands[self.strand_id].domains[0].start;
        let id = self.strand_id;
        self.builder.strands[id].cyclic = true;
        self.cross(helix).to(start);
    }
}

fn sequence<StrandLabel, DomainLabel>(
    positions: &HashMap<(usize, isize, bool), char>,
    strand: &Strand<StrandLabel, DomainLabel>,
) -> String {
    let mut seq = String::new();
    // println!("label = {:?}", strand.label);
    for dom in strand.domains.iter() {
        if let Some(ref s) = dom.sequence {
            seq.push_str(s);
            continue;
        }
        assert!(dom.helix >= 0);
        for pos in dom.iter() {
            seq.push(sequences::complement(
                *positions
                    .get(&(dom.helix as usize, pos, !dom.right))
                    .unwrap(),
            ))
        }
    }
    seq
}

/// Common geometric transforms on domains (most useful when working
/// on grids or partial grids).
pub mod transforms {
    use super::Domain;

    /// Rotate all the domains provided as input by 180°. This can be
    /// used on domains obtained with
    /// [`StrandRef::domains()`](../struct.StrandRef.html#method.domains)
    /// method.
    pub fn rotate<Label>(domains: &mut [Domain<Label>]) {
        let max_helix = domains.iter().map(|x| x.helix).max().unwrap();
        for dom in domains.iter_mut() {
            let start = dom.start;
            dom.start = -dom.end;
            dom.end = -start;
            dom.helix = max_helix - dom.helix;
            dom.right = !dom.right;
        }
    }

    /// Flip all the domains provided as input vertically, i.e. relative
    /// to a horizontal axis. This can be used on domains obtained
    /// with
    /// [`StrandRef::domains()`](../struct.StrandRef.html#method.domains)
    /// method.
    pub fn v_flip<Label>(domains: &mut [Domain<Label>], max_helix: isize) {
        for dom in domains.iter_mut() {
            dom.helix = max_helix - dom.helix;
        }
    }

    /// Reverse the direction of the list of domains provided as
    /// input, starting from the very last base back to the first
    /// one. This can be used on domains obtained with
    /// [`StrandRef::domains()`](../struct.StrandRef.html#method.domains)
    /// method.
    pub fn reverse<Label>(domains: &mut [Domain<Label>]) {
        for dom in domains.iter_mut() {
            dom.right = !dom.right;
        }
    }

    /// Flip all the domains provided as input horizontally,
    /// i.e. relative to a vertical axis. This can be used on domains
    /// obtained with
    /// [`StrandRef::domains()`](../struct.StrandRef.html#method.domains)
    /// method.
    pub fn h_flip<Label>(domains: &mut [Domain<Label>]) {
        for dom in domains.iter_mut() {
            let start = dom.start;
            dom.start = -dom.end;
            dom.end = -start;
        }
    }

    /// Trnsalate all the domains provided as input by the specified
    /// number of helices (parameter `dh`) and bases (parameter
    /// `dx`). This can be used on domains obtained with
    /// [`StrandRef::domains()`](../struct.StrandRef.html#method.domains)
    /// method.
    pub fn translate<Label>(domains: &mut [Domain<Label>], dh: isize, dx: isize) {
        use std::convert::TryFrom;
        for dom in domains.iter_mut() {
            dom.helix = usize::try_from(dom.helix as isize + dh).unwrap() as isize;
            dom.start += dx;
            dom.end += dx;
        }
    }
}

impl<StrandLabel, DomainLabel> Design<StrandLabel, DomainLabel> {
    fn positions(&self) -> HashMap<(usize, isize, bool), char> {
        // First add sequences to all positions
        let mut positions = HashMap::new();
        for s in self.strands.iter() {
            if let Some(ref seq) = s.sequence {
                let mut seq = seq.chars();
                for dom in s.domains.iter() {
                    if dom.helix < 0 {
                        continue;
                    }
                    for pos in dom.iter() {
                        let c = seq.next().unwrap();
                        positions.insert((dom.helix as usize, pos, dom.right), c);
                    }
                }
            }
        }
        positions
    }

    /// Output all the sequences of the strands whose sequence has not
    /// been forced.
    pub fn sequences(&self) -> Vec<String> {
        let positions = self.positions();
        let mut sequences = Vec::new();
        for s in self.strands.iter() {
            if s.sequence.is_none() {
                sequences.push(sequence(&positions, s))
            }
        }
        sequences
    }
}

#[cfg(feature = "excel")]
impl<StrandLabel: Ord + std::fmt::Display, DomainLabel>
    Design<StrandLabel, DomainLabel>
{
    /// Output plates with the sequences, ordered by label, as an
    /// Excel spreadsheet that can be ordered from [the IDT
    /// website](https://idtdna.com).
    ///
    /// All strands `s` such that `filter(s).is_some()` are included,
    /// and only those. Moreover, if `filter(s) == Some(true)`, a new
    /// plate is started if the current plate is nonempty.
    ///
    /// Moreover, strands are ran through the `filter` function in the
    /// same order as they are in `self`.
    pub fn make_plates_96<
            H: FnMut(&Strand<StrandLabel, DomainLabel>) -> usize,
        F: FnMut(&Strand<StrandLabel, DomainLabel>) -> Option<bool>,
        G: FnMut(usize) -> String,
    >(
        &mut self,
        file: &str,
        mut seq_class: H,
        mut filter: F,
        mut format: G,
    ) {
        self.strands.sort_by(|a, b| seq_class(&a).cmp(&seq_class(b)));
        let mut current_plate = 0;
        let mut current_class = 0;

        let mut plate = Plate96::new();
        let mut plates = vec![Vec::new()];
        let positions = self.positions();
        for s in self.strands.iter() {
            let class = seq_class(s);
            if class != current_class && plate.row > 0 {
                plate.incr_col();
            }

            if plate.n != current_plate {
                plates.push(Vec::new());
                current_plate = plate.n;
            }
            if let Some(x) = filter(s) {
                // println!("{:?}", plate);
                if x && (plate.row != 0 || plate.column != 0) {
                    plate.incr_plate();
                    plates.push(Vec::new());
                    current_plate = plate.n;
                }
                plates.last_mut().unwrap().push((
                    plate.column,
                    plate.row,
                    &s.label,
                    sequence(&positions, s),
                ));
                plate.incr()
            }
            current_class = class;
        }
        use simple_excel_writer::*;
        let mut wb = Workbook::create(file);
        for (i, seqs) in plates.iter().enumerate() {
            let mut sheet = wb.create_sheet(&format(i));
            sheet.add_column(Column { width: 30.0 });
            sheet.add_column(Column { width: 30.0 });
            sheet.add_column(Column { width: 30.0 });
            wb.write_sheet(&mut sheet, |sw| {
                sw.append_row(row!["Well Position", "Name", "Sequence"])?;
                for &(col, row, label, ref seq) in seqs.iter() {
                    let row = (b'A' + row as u8) as char;
                    let label = if let Some(ref label) = label {
                        format!("{}", label)
                    } else {
                        format!("{}{}", row, col + 1)
                    };
                    sw.append_row(row![
                        format!("{}{}", row, col + 1).as_str(),
                        label.as_str(),
                        seq.as_str()
                    ])?
                }
                Ok(())
            })
            .unwrap()
        }
        wb.close().unwrap();
    }
}

#[cfg(feature = "excel")]
#[derive(Debug, Copy, Clone)]
/// Representation of a 96-well plate being filled. Plates are filled
/// from top to bottom (i.e. rows are incremented first), then from
/// left to right (columns are incremented).
pub struct Plate96 {
    /// Plate number.
    pub n: usize,
    /// Column number.
    pub column: usize,
    /// Row number.
    pub row: usize,
}

#[cfg(feature = "excel")]
impl Plate96 {
    /// Create a plate
    pub fn new() -> Self {
        Plate96 {
            n: 0,
            column: 0,
            row: 0,
        }
    }
    /// Move to the next well.
    pub fn incr(&mut self) {
        self.row += 1;
        if self.row >= 8 {
            self.incr_col()
        }
    }
    /// Move to the next column.
    pub fn incr_col(&mut self) {
        self.column += 1;
        self.row = 0;
        if self.column >= 12 {
            self.incr_plate()
        }
    }
    /// Move to the next plate.
    pub fn incr_plate(&mut self) {
        self.column = 0;
        self.row = 0;
        self.n += 1;
    }
}