use super::*;
type IntervalVec = Vec<Option<(isize, isize)>>;
type PositionVec = Vec<Option<[f64; 3]>>;
impl Helix {
fn optimise_relative(
&mut self,
other: &mut Helix,
xovers: &Vec<(isize, bool, isize, bool)>,
cst: &Parameters,
id_other: usize,
positions: &PositionVec,
intervals: &IntervalVec,
) {
let (contact_self, contact_other) =
self.optimise_relative_angle(xovers, cst, id_other, positions, intervals);
let delta = (1., 0.);
other.pitch = self.pitch;
other.yaw = if delta.0 > 0. {
self.yaw
} else {
PI + self.yaw
};
let roll_other = PI + contact_self - contact_other + self.roll;
other.roll = if delta.0 > 0. {
roll_other
} else {
PI - roll_other
};
let r = cst.helix_radius + cst.inter_helix_gap / 2.;
let contact = [
0.,
r * (contact_self + self.roll).cos(),
r * (contact_self + self.roll).sin(),
];
let forward = self.rotate_point([delta.1, 0., 0.]);
let contact_point = self.rotate_point(contact);
other.position.x = self.position.x + forward[0] + 2. * contact_point[0];
other.position.y = self.position.y + forward[1] + 2. * contact_point[1];
other.position.z = self.position.z + forward[2] + 2. * contact_point[2];
}
fn optimise_relative_angle(
&self,
xovers: &Vec<(isize, bool, isize, bool)>,
cst: &Parameters,
other_id: usize,
positions: &PositionVec,
intervals: &IntervalVec,
) -> (f64, f64) {
let angles = xovers_to_angle(xovers, cst);
let mut opt_phi = 0.;
let mut opt_theta = 0.;
let mut opt = std::f64::INFINITY;
for i in 0..100 {
let phi = 2. * i as f64 * PI / 100.;
let r = cst.helix_radius + cst.inter_helix_gap / 2.;
let contact = [0., r * (phi + self.roll).cos(), r * (phi + self.roll).sin()];
let contact_point = self.rotate_point(contact);
let other_position = [
self.position.x + 2. * contact_point[0],
self.position.y + 2. * contact_point[1],
self.position.z + 2. * contact_point[2],
];
if collision(other_id, other_position, intervals, positions, 2. * r) {
continue;
}
for j in 0..100 {
let theta = 2. * j as f64 * PI / 100.;
let guess = evaluate_angle(phi, theta, &angles);
if guess < opt {
opt_phi = phi;
opt_theta = theta;
opt = guess;
}
}
}
(opt_phi, opt_theta)
}
}
fn collision(
other_id: usize,
other_position: [f64; 3],
intervals: &IntervalVec,
positions: &PositionVec,
objective: f64,
) -> bool {
if intervals[other_id].is_none() {
return false;
}
let (x1, x2) = intervals[other_id].unwrap();
for i in 0..intervals.len() {
if positions[i].is_some() && intervals[i].is_some() {
let (y1, y2) = intervals[i].unwrap();
if x1 <= y2 && y1 <= x2 {
let vec = Point::from_coord(other_position) - Point::from_coord(positions[i].unwrap());
if vec.norm() < objective {
return true;
}
}
}
}
false
}
fn evaluate_angle(contact_a: f64, contact_b: f64, angles: &Vec<(f64, f64)>) -> f64 {
let convert = |x: f64| PI / 2. - x;
let roll_b = contact_a + PI - contact_b;
angles
.iter()
.map(|(alpha, beta)| {
let y_a = convert(*alpha).sin();
let y_b = convert(beta + roll_b).sin() + 2. * convert(contact_a).sin();
let z_a = convert(*alpha).cos();
let z_b = convert(beta + roll_b).cos() + 2. * convert(contact_a).cos();
let delta_y = y_a - y_b;
let delta_z = z_a - z_b;
let len = (delta_y * delta_y + delta_z * delta_z).sqrt();
len
})
.fold(0., |acc, x| acc.max(x))
}
fn xovers_to_angle(xovers: &Vec<(isize, bool, isize, bool)>, cst: &Parameters) -> Vec<(f64, f64)> {
let theta = |n, right| {
let shift = if right { cst.groove_angle } else { 0. };
n as f64 * 2. * PI / cst.bases_per_turn + shift
};
(0..xovers.len())
.map(|i| {
let (x, b1, y, b2) = xovers[i];
(theta(x, b1), theta(y, b2))
})
.collect()
}
impl<S, T> Design<S, T> {
pub fn separate_helix(&mut self, n: usize) {
println!("helix {}", n);
let mut intervals = Vec::new();
for s in &self.strands {
for d in &s.domains {
if d.length() > 0 && d.helix == n as isize {
intervals.push((d.start, d.end));
}
}
}
if intervals.len() <= 1 {
return;
}
println!("intervals {}\n {:?}", n, intervals);
intervals.sort();
println!("intervals sorted {}\n {:?}", n, intervals);
let mut merged_intervals = Vec::new();
merged_intervals.push(intervals[0]);
for i in 1..intervals.len() {
let interval_1 = merged_intervals.last_mut().expect("last mut");
let interval_2 = intervals[i];
if interval_1.1 >= interval_2.0 {
interval_1.1 = interval_2.1;
} else {
merged_intervals.push(interval_2);
}
}
println!("intervals merged {:?}", merged_intervals);
if merged_intervals.len() <= 1 {
return;
}
let mut new_helix = Vec::new();
new_helix.push(n);
for _ in 1..merged_intervals.len() {
new_helix.push(self.helices.len());
self.helices.push(self.helices[n].clone());
}
for s in &mut self.strands {
for d in &mut s.domains {
if d.length() > 0 && d.helix == n as isize {
println!("start {}", d.start);
println!("merged {:?}", merged_intervals);
let mut h = 0;
while d.start > merged_intervals[h].1 {
h += 1;
}
d.helix = new_helix[h] as isize;
}
}
}
}
pub fn get_interval_helix(&self, n: usize) -> Option<(isize, isize)> {
let mut intervals = Vec::new();
for s in &self.strands {
for d in &s.domains {
if d.length() > 0 && d.helix == n as isize {
intervals.push((d.start, d.end));
}
}
}
if intervals.len() == 0 {
return None;
} else if intervals.len() == 1 {
return Some(intervals[0]);
}
intervals.sort();
let mut merged_intervals = Vec::new();
merged_intervals.push(intervals[0]);
for i in 1..intervals.len() {
let interval_1 = merged_intervals.last_mut().expect("last mut");
let interval_2 = intervals[i];
if interval_1.1 >= interval_2.0 {
interval_1.1 = interval_2.1;
} else {
merged_intervals.push(interval_2);
}
}
return Some((
merged_intervals.first().unwrap().0,
merged_intervals.last().unwrap().1,
));
}
pub fn separate_all_helices(&mut self) {
for i in 0..self.helices.len() {
self.separate_helix(i);
}
}
fn get_parralel_graph(&self) -> Vec<Vec<usize>> {
let nb_helix = self.helices.len();
let mut ret = vec![Vec::new(); nb_helix];
let xovers = self.get_cross_overs();
for i in 0..nb_helix {
for j in (i + 1)..nb_helix {
let nb_cross = xovers
.iter()
.filter(|(h1, _, _, h2, _, _)| {
(*h1 == i as isize && *h2 == j as isize)
|| (*h2 == i as isize && *h1 == j as isize)
})
.count();
if nb_cross >= 2 {
ret[i].push(j);
ret[j].push(i);
}
}
}
ret
}
pub fn guess_bfs(&mut self, separate: bool) {
if separate {
self.separate_all_helices();
}
let graph = self.get_parralel_graph();
let intervals: IntervalVec = (0..self.helices.len())
.map(|i| self.get_interval_helix(i))
.collect();
println!("graph_got");
let nb_helix = self.helices.len();
let mut seen = vec![false; nb_helix];
let mut to_do = VecDeque::new();
let param = &self.parameters.unwrap();
let mut position: PositionVec = vec![None; nb_helix];
for i in 0..nb_helix {
if !seen[i] {
seen[i] = true;
position[i] = {
let position = self.helices[i].position;
Some([position.x, position.y, position.z])
};
to_do.push_back((i, -1));
while to_do.len() > 0 {
let (target, source) = to_do.pop_front().unwrap();
if source >= 0 {
let xovers = self.get_cross_overs_helices(source as isize, target as isize);
let source = source as usize;
if source < target {
let split = self.helices.split_at_mut(target);
split.0[source].optimise_relative(
&mut split.1[0],
&xovers,
param,
target,
&position,
&intervals,
);
} else {
let split = self.helices.split_at_mut(source);
split.1[0].optimise_relative(
&mut split.0[target],
&xovers,
param,
target,
&position,
&intervals,
);
}
position[target] = {
let position = self.helices[target].position;
Some([position.x, position.y, position.z])
};
}
for k in &graph[target] {
let k = *k;
if !seen[k] {
seen[k] = true;
to_do.push_back((k, target as isize));
}
}
}
}
}
}
}