use crate::{log, Point};
use std::io::Write;
/// The three variants of the algorithm.
#[derive(PartialEq, Debug, Clone, Copy)]
pub enum Variant {
/// The starndard version terminates a simplex reduction when no other reduced simplex shares it's
/// highest boundary simplex.
Standard,
/// The exhaustive variant tries to lower the indices of the reduces simplcides after the max
/// has been found.
Exhaustive,
}
use self::Variant::*;
/// One simplex
#[derive(Debug, Clone)]
pub struct Simplex {
/// Original index. We need this for bookkeeping.
pub j: usize,
/// Indices of the faces of this simplex
pub faces: Vec<usize>,
// The ball radius when this simplex is born.
// pub r_when_born: f64,
pub dim: isize,
}
impl Simplex {
/// Get the dimension of this simplex.
pub fn dim(&self) -> isize {
self.dim
}
/// Make the empty simplex
pub fn empty() -> Self {
Simplex {
j: 0,
faces: Vec::new(),
dim: -1,
}
}
}
/// Data we use for stuff.
pub struct Persistence {
/// All simplices in the data set
pub simplices: Vec<Simplex>,
/// Point location data, for rendering.
pub points: Vec<[f64; 2]>,
/// Alpha values for all simplices.
/// TODO: reconsider having this here; what if we don't have any alpha?
pub alphas: Vec<f64>,
}
/// The result from reducing.
pub struct Reduction<'a> {
/// When we reduce a simplex we might finding that it's `low` is unique (that `simplex_with_low[low] == None`).
/// We then store the reduced simplex in that slot. This corresponds to pairing `low` with `j`.
pub pairings: Vec<(usize, usize)>,
/// This maps an index to the reduced simplex which maximum boundary element is that index.
pub simplex_with_low: Vec<Option<Simplex>>,
/// Here we track which elements are added to which columns.
/// This was probably used for the output?
pub additions: Vec<(usize, usize)>,
pub log: log::Reduction<'a>,
}
impl<'a> Reduction<'a> {
/// Write debug output for this reduction to `w`. This includes how each simplex was reduced.
pub fn debug_output<W: Write>(
&self,
w: &mut W,
simplices: &[Simplex],
) -> Result<(), std::io::Error> {
use std::collections::HashMap;
let mut added_to = HashMap::<usize, Vec<usize>>::new();
for &(left, right) in &self.additions {
added_to.entry(left).or_insert(Vec::new()).push(right);
}
let mut killed = HashMap::new();
for &(victim, killer) in &self.pairings {
killed.insert(killer, victim);
}
writeln!(w, "j;faces;adds;reduced")?;
for s in simplices {
let empty_vec = Vec::new();
let j = s.j;
let mut initial_faces = s.faces.clone();
initial_faces.sort();
let reduced_faces = if let Some(&place) = killed.get(&j) {
let v = self.simplex_with_low[place]
.as_ref()
.expect("Simplex was killing, but wasn't stored?");
&v.faces
} else {
&empty_vec
};
let adds = added_to.get(&s.j).unwrap_or(&empty_vec);
writeln!(
w,
"{};{:?};{:?};{:?}",
j, initial_faces, adds, reduced_faces
)?;
}
Ok(())
}
}
/// Run the reduction algorithm on the given `Persistence`. The `variant` is ran,
/// and stats are recorded in `stats`.
pub fn reduce(p: &Persistence, variant: Variant) -> Reduction {
let mut pairings = Vec::new();
let mut additions = Vec::new();
let mut reduced_with_low: Vec<Option<Simplex>> = vec![None; p.simplices.len()];
let mut reduction_log = log::Reduction::new(p);
// We're reducing each simplex.
for simplex in &p.simplices {
// we'll store an index here, if we don't fully reduce the simplex
let mut place_index = None;
// With no boundary we're done.
if simplex.faces.last().is_none() {
continue;
}
let mut current_reduced = simplex.clone();
current_reduced.faces.sort();
let mut simplex_log = log::SimplexReduction::new(&simplex);
while current_reduced.faces.len() > 0 {
let low = *current_reduced.faces.last().unwrap();
if let Some(ref other) = reduced_with_low[low] {
additions.push((current_reduced.j, other.j));
simplex_log.add(log::Add::regular(¤t_reduced, other));
simplex_reduce(&mut current_reduced, other);
} else {
place_index = Some(low);
pairings.push((low, simplex.j));
simplex_log.place_index = place_index;
break;
}
}
if variant == Exhaustive {
// Reduce the rest of the simplex. Go backwards. We can skip the last index, since
// that's us. Then we count the numbers of boundary simplices we know we can't reduce
// any further. Since we're going bottom up and looking at `reduced_with_low`, these
// cannot change with later redutions.
let mut safe_simplices = 1;
let mut n = current_reduced.faces.len();
while safe_simplices < n {
let cand_i = current_reduced.faces[n - 1 - safe_simplices];
if let Some(ref cand) = reduced_with_low[cand_i] {
additions.push((current_reduced.j, cand.j));
simplex_log.add(log::Add::exhaustive(¤t_reduced, cand));
simplex_reduce(&mut current_reduced, cand);
n = current_reduced.faces.len();
} else {
safe_simplices += 1;
}
}
}
if current_reduced.faces.len() != 0 {
let i = place_index.expect(
"Reduced simplex isn't empty, yet we haven't set the index in which to store it",
);
reduced_with_low[i] = Some(current_reduced);
}
reduction_log.add(simplex_log);
}
Reduction {
pairings,
simplex_with_low: reduced_with_low,
additions,
log: reduction_log,
}
}
/// Reduce `this` by joining it with `other`.
fn simplex_reduce(this: &mut Simplex, other: &Simplex) {
// TODO: This is terrible!
let mut slow_vec = Vec::new();
let mut _buffer = [0usize; 32];
let num_inds = this.faces.len() + other.faces.len();
let fastpath = num_inds <= 32;
let mut buf: &mut [usize] = &mut _buffer;
if !fastpath {
// We might need more space. Allocate.
// TODO: this is awkward
slow_vec = vec![0usize; num_inds];
buf = &mut slow_vec;
}
let mut our_i = 0;
let mut their_i = 0;
let mut buffer_i = 0;
let our_last = this.faces.len();
let their_last = other.faces.len();
while our_i < our_last && their_i < their_last {
if this.faces[our_i] < other.faces[their_i] {
buf[buffer_i] = this.faces[our_i];
our_i += 1;
buffer_i += 1;
} else if this.faces[our_i] == other.faces[their_i] {
our_i += 1;
their_i += 1;
} else {
buf[buffer_i] = other.faces[their_i];
their_i += 1;
buffer_i += 1;
}
}
while our_i < our_last {
buf[buffer_i] = this.faces[our_i];
our_i += 1;
buffer_i += 1;
}
while their_i < their_last {
buf[buffer_i] = this.faces[their_i];
their_i += 1;
buffer_i += 1;
}
if fastpath {
this.faces.reserve(num_inds);
this.faces.clear();
for &n in &buf[..buffer_i] {
this.faces.push(n);
}
} else {
std::mem::swap(&mut this.faces, &mut slow_vec);
}
}
/// Compute the alpha values for the simplices, given the point location.
/// As the name suggests, this is only for dimension 2 and less.
pub fn compute_alpha_values_2d(
points: &[[f64; 2]],
simplices: &[Simplex],
) -> (Vec<Simplex>, Vec<f64>) {
fn alpha_1d(simplex: &Simplex, points: &[[f64; 2]]) -> f64 {
let p0 = Point(points[simplex.faces[0]]);
let p1 = Point(points[simplex.faces[1]]);
let diff = p1 - p0;
0.5 * diff.len()
}
fn alpha_2d(
simplex: &Simplex,
points: &[[f64; 2]],
simplices: &[Simplex],
) -> (f64, Option<usize>) {
let xy_i = simplex.faces[0];
let xy = &simplices[xy_i];
let yz_i = simplex.faces[1];
let _yz = &simplices[yz_i];
let zx_i = simplex.faces[2];
let zx = &simplices[zx_i];
// The labels here may be messed up due to the sorting when making the edges.
// Get it right.
let (x_i, y_i, z_i) = if xy.faces[0] == zx.faces[0] {
// xy, xz
(xy.faces[0], xy.faces[1], zx.faces[1])
} else if xy.faces[0] == zx.faces[1] {
// xy, zx
(xy.faces[0], xy.faces[1], zx.faces[0])
} else if xy.faces[1] == zx.faces[0] {
// yx, xz
(xy.faces[1], xy.faces[0], zx.faces[1])
} else {
// yx, zx
(xy.faces[1], xy.faces[0], zx.faces[0])
};
let x = Point(points[x_i - 1]);
let y = Point(points[y_i - 1]);
let z = Point(points[z_i - 1]);
let xy = (y - x).len();
let xz = (z - x).len();
let yz = (z - y).len();
let ax = (y - x).angle(z - x);
let ay = (z - y).angle(x - y);
let az = (x - z).angle(y - z);
let am = ax.max(ay).max(az);
// The circumcenter is the point in which the balls from each vertex meets, so
// at `r=circumradius` the face is filled in.
let circumradius = (xy * yz * xz)
/ ((xy + yz + xz) * (xy + yz - xz) * (xy - yz + xz) * (-xy + yz + xz)).sqrt();
let mut edit = None;
if am * 180.0 / std::f64::consts::PI >= 90.0 {
// When one angle is > 90.0, the circumcenter of the triangle is outside. Then
// the balls from the edges that's closest to the circumcenter doesn't meet
// before they are at the circumcenter. Thus we must adjust the `r_when_born` on
// this edge.
if am == ax {
edit = Some(yz_i);
} else if am == ay {
edit = Some(zx_i);
} else {
edit = Some(xy_i);
}
}
(circumradius, edit)
}
#[derive(PartialEq)]
struct CmpSimplex {
dim: usize,
alpha: f64,
index: usize,
}
impl Eq for CmpSimplex {}
impl PartialOrd for CmpSimplex {
fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
Some(self.cmp(other))
}
}
impl Ord for CmpSimplex {
fn cmp(&self, other: &Self) -> std::cmp::Ordering {
self.dim.cmp(&other.dim).then(
self.alpha
.partial_cmp(&other.alpha)
.expect("alpha value was NaN or something!"),
)
}
}
let mut hacked_alphas = std::collections::HashMap::new();
let mut v = Vec::with_capacity(simplices.len());
for current_simplex_i in 0..simplices.len() {
let simplex = &simplices[current_simplex_i];
let alpha = match simplex.dim() {
-1 | 0 => 0.0,
1 => alpha_1d(&simplex, &points),
2 => {
let (a, edit) = alpha_2d(&simplex, &points, simplices);
if let Some(i) = edit {
hacked_alphas.insert(i, a);
}
a
}
more => panic!("Only support 2d simplices! ({} > 2)", more),
};
v.push(CmpSimplex {
dim: simplex.faces.len(),
alpha,
index: current_simplex_i,
});
}
v.sort();
let mut ordered: Vec<Simplex> = Vec::with_capacity(simplices.len());
let mut alphas: Vec<f64> = Vec::with_capacity(simplices.len());
for cs in v.into_iter() {
ordered.push(simplices[cs.index].clone());
if let Some(&a) = hacked_alphas.get(&cs.index) {
alphas.push(a);
} else {
alphas.push(cs.alpha);
}
}
// Now we have simplices in a new order. Next we need to rewrite the faces of all simplices,
// since the indices are in the old ordering. To do this we build an index mapping from old to new
// index.
let mut old_to_new = vec![0; simplices.len()];
for (i, s) in ordered.iter().enumerate() {
old_to_new[s.j] = i;
}
// Now we rewrite all faces.
for (i, s) in ordered.iter_mut().enumerate() {
s.j = i;
for j in s.faces.iter_mut() {
*j = old_to_new[*j];
}
}
(ordered, alphas)
}