codenano 0.5.1

/// This type represents an abstract OpenGL drawing. After it is built
/// using the `new` and `reset` methods, its `points` field can be
/// sent to OpenGL as a triangle strip for rendering.
pub struct OpenGL {
    /// Points of the vertices to render.
    pub points: Vec<f32>,
    geom: bool,
    /// The size of a single point, in f32 (i.e. the size in bytes is four times this field).
    pub size: usize,
    current: Option<(Vector, Vector, Vector)>,
    /// Number of subdivisions (horizontal and vertical) in a sphere.
    pub sphere_n: usize,
    /// Number of subdivisions in the backbone cylinders.
    pub backbone_n: usize,
    /// Number of subdivisions in the basepair cylinders.
    pub pair_n: usize,
}

impl OpenGL {
    /// Create a new `OpenGL`. This does not load any design, but
    /// merely creates the necessary structures.
    pub fn new(
        use_geometry_shader: bool,
        sphere_n: usize,
        backbone_n: usize,
        pair_n: usize,
    ) -> Self {
        OpenGL {
            size: if use_geometry_shader { 5 } else { 8 },
            geom: use_geometry_shader,
            points: Vec::new(),
            current: None,
            sphere_n,
            backbone_n,
            pair_n,
        }
    }

    /// Reset this type and load a new design. This must also be used
    /// initially after calling the `new` method.
    pub fn reset<A, B>(&mut self, design: &crate::Design<A, B>) {
        self.points.clear();
        self.current = None;
        for i in 0..design.strands.len() {
            self.add_strand(design, self.geom, i);
        }
        let mut nucl = HashMap::new();
        for i in 0..design.strands.len() {
            self.add_pairs(&design, self.geom, &mut nucl, i);
        }
    }


    fn push_p64(&mut self, x: [f64; 3]) {
        self.points.push(x[0] as f32);
        self.points.push(x[1] as f32);
        self.points.push(x[2] as f32);
    }

    fn push_int(&mut self, x: u32) {
        self.points.push(unsafe { std::mem::transmute(x) });
    }

    fn add_cylinder(&mut self, i: usize, point: [f64; 3], color: u32, r: f64, n: usize) {
        assert!(!self.geom);
        let p = Vector([point[0], point[1], point[2]]);
        if let Some((p0, vx0, vy0)) = self.current {
            let vec = (p - p0).normalize();
            let vx = vec.normal();
            let vy = vec.vec_product(vx);
            let _vx0 = if vx0.is_null() { vx } else { vx0 };
            let _vy0 = if vy0.is_null() { vy } else { vy0 };
            let mut t = Vector::NULL;
            let mut normal = Vector::NULL;
            let mut is_first = true;
            for rot in 0..=2 * n + 1 {
                let angle = -((rot / 2) as f64) * 2. * PI / (n as f64);
                normal = r * vx * angle.cos() + r * vy * angle.sin();
                t = if rot % 2 == 0 {
                    p0 + normal
                } else {
                    p + normal
                };
                if is_first {
                    self.push_p64(t.0);
                    self.push_p64(normal.normalize().0);
                    self.push_int(i as u32);
                    self.push_int(color);
                    is_first = false;
                }
                self.push_p64(t.0);
                self.push_p64(normal.normalize().0);
                self.push_int(i as u32);
                self.push_int(color);
            }
            self.push_p64(t.0);
            self.push_p64(normal.normalize().0);
            self.push_int(i as u32);
            self.push_int(color);
            self.current = Some((p, vx, vy));
        } else {
            self.current = Some((p, Vector::NULL, Vector::NULL));
        }
    }

    fn add_sphere(&mut self, k: usize, point: [f64; 3], color: u32, r: f64, n: usize) {
        assert!(!self.geom);
        let p = Vector([point[0], point[1], point[2]]);
        let mut last = p;
        let mut normal = Vector::NULL;
        for i in 0..n {
            let alpha = (i as f64) * PI / (n as f64);
            let alpha_ = ((i + 1) as f64) * PI / (n as f64);
            let r0 = r * alpha.sin();
            let r1 = r * alpha_.sin();
            for j in 0..=2 * n + 1 {
                let beta = -2. * ((j / 2) as f64) * PI / (n as f64);
                normal = if j & 1 == 0 {
                    Vector([r0 * beta.sin(), -r * alpha.cos(), r0 * beta.cos()])
                } else {
                    Vector([r1 * beta.sin(), -r * alpha_.cos(), r1 * beta.cos()])
                };
                last = p + normal;
                if j == 0 {
                    self.push_p64(last.0);
                    self.push_p64(normal.normalize().0);
                    self.push_int(k as u32);
                    self.push_int(color);
                }
                self.push_p64(last.0);
                self.push_p64(normal.normalize().0);
                self.push_int(k as u32);
                self.push_int(color)
            }
            self.push_p64(last.0);
            self.push_p64(normal.normalize().0);
            self.push_int(k as u32);
            self.push_int(color)
        }
    }

    fn add_strand<A, B>(&mut self, design: &crate::Design<A, B>, geom: bool, i: usize) {
        let ref s = design.strands[i];
        self.current = None;
        let color = if let Some(ref color) = s.color {
            color.as_int()
        } else {
            0x888888
        };
        for d in s.domains.iter() {
            if d.helix < 0 {
                continue;
            }
            let h = d.helix as usize;
            if d.forward {
                for x in d.start..d.end {
                    let point = design.helices[h].space_pos(
                        design
                            .parameters
                            .as_ref()
                            .unwrap_or(&crate::Parameters::DEFAULT),
                        x,
                        true,
                    );
                    if geom {
                        self.push_p64(point);
                        self.push_int(i as u32);
                        self.push_int(color);
                    } else {
                        self.add_cylinder(i, point, color, 0.1, self.backbone_n);
                        self.add_sphere(i, point, color, 0.2, self.sphere_n);
                    }
                }
            } else {
                for x in (d.start..d.end).rev() {
                    let point = design.helices[h].space_pos(
                        design
                            .parameters
                            .as_ref()
                            .unwrap_or(&crate::Parameters::DEFAULT),
                        x,
                        false,
                    );
                    if geom {
                        self.push_p64(point);
                        self.push_int(i as u32);
                        self.push_int(color);
                    } else {
                        self.add_cylinder(i, point, color, 0.1, self.backbone_n);
                        self.add_sphere(i, point, color, 0.2, self.sphere_n);
                    }
                }
            }
        }
        // Make an infinite last point to handle the last segment of
        // the linestrip.
        if geom {
            self.push_p64([1. / 0., 0., 0.]);
            self.push_int(i as u32);
            self.push_int(color);
        }
    }

    fn add_pairs<A, B>(
        &mut self,
        design: &crate::Design<A, B>,
        geom: bool,
        nucl: &mut HashMap<(isize, isize, bool), usize>,
        i: usize,
    ) {
        let ref s = design.strands[i];
        for dom in s.domains.iter() {
            if dom.helix < 0 {
                continue;
            }
            for x in dom.start..dom.end {
                nucl.insert((dom.helix, x, dom.forward), i);
                if let Some(&j) = nucl.get(&(dom.helix, x, !dom.forward)) {
                    let a = design.helices[dom.helix as usize].space_pos(
                        design
                            .parameters
                            .as_ref()
                            .unwrap_or(&crate::Parameters::DEFAULT),
                        x,
                        dom.forward,
                    );
                    let b = design.helices[dom.helix as usize].space_pos(
                        design
                            .parameters
                            .as_ref()
                            .unwrap_or(&crate::Parameters::DEFAULT),
                        x,
                        !dom.forward,
                    );
                    let mut a = Vector([a[0], a[1], a[2]]);
                    let mut b = Vector([b[0], b[1], b[2]]);
                    let v = (b - a).normalize();
                    a = a + v * 0.19;
                    b = b - v * 0.19;
                    if geom {
                        let k = (i << 16) | j;
                        self.push_p64(a.0);
                        self.push_int(k as u32);
                        self.push_int(0xaaaaaa);
                        self.push_p64(b.0);
                        self.push_int(k as u32);
                        self.push_int(0xaaaaaa);

                        self.push_p64([1. / 0., 0., 0.]);
                        self.push_int(k as u32);
                        self.push_int(0);
                    } else {
                        self.current = Some((a, Vector::NULL, Vector::NULL));
                        self.add_cylinder((i << 16) | j, b.0, 0xaaaaaa, 0.05, self.pair_n);
                    }
                }
            }
        }
    }
}

use std::f64::consts::PI;
use std::collections::HashMap;

#[derive(Debug, Clone, Copy, PartialEq)]
struct Vector(pub [f64; 3]);

impl Vector {
    /// Return one vector orthogonal to `self`.
    fn normal(self) -> Self {
        // Find the smallest coordinate in absolute value (and avoid
        // it, for numerical stability).
        let mut arg = 0;
        if self.0[1].abs() < self.0[0].abs() {
            arg = 1;
        }
        if self.0[2].abs() < self.0[arg] {
            arg = 2;
        }
        (Vector(if arg == 0 {
            [0., self.0[2], -self.0[1]]
        } else if arg == 1 {
            [self.0[2], 0., -self.0[0]]
        } else {
            [self.0[1], -self.0[0], 0.]
        }))
        .normalize()
    }

    /// Normalize `self`, i.e. divide it by its norm.
    fn normalize(&self) -> Self {
        let norm = (self.0[0] * self.0[0] + self.0[1] * self.0[1] + self.0[2] * self.0[2]).sqrt();
        Vector([self.0[0] / norm, self.0[1] / norm, self.0[2] / norm])
    }

    /// Vector product of two vectors.
    fn vec_product(&self, b: Self) -> Self {
        Vector([
            self.0[1] * b.0[2] - self.0[2] * b.0[1],
            self.0[2] * b.0[0] - self.0[0] * b.0[2],
            self.0[0] * b.0[1] - self.0[1] * b.0[0],
        ])
    }

    const NULL: Vector = Vector([0., 0., 0.]);

    /// Is this the null vector?
    fn is_null(&self) -> bool {
        self.0 == Self::NULL.0
    }
}

impl std::ops::Sub<Vector> for Vector {
    type Output = Vector;
    fn sub(self, b: Vector) -> Self::Output {
        Vector([self.0[0] - b.0[0], self.0[1] - b.0[1], self.0[2] - b.0[2]])
    }
}

impl std::ops::Add<Vector> for Vector {
    type Output = Vector;
    fn add(self, b: Vector) -> Self::Output {
        Vector([self.0[0] + b.0[0], self.0[1] + b.0[1], self.0[2] + b.0[2]])
    }
}

impl std::ops::Mul<f64> for Vector {
    type Output = Vector;
    fn mul(self, b: f64) -> Self::Output {
        Vector([self.0[0] * b, self.0[1] * b, self.0[2] * b])
    }
}

impl std::ops::Mul<Vector> for f64 {
    type Output = Vector;
    fn mul(self, b: Vector) -> Self::Output {
        Vector([b.0[0] * self, b.0[1] * self, b.0[2] * self])
    }
}