~mht/cra

2fa28e40a0c7fcd01a655929ad159eea85ee897d — Martin Hafskjold Thoresen 1 year, 3 months ago 1910306
Move persistence things to new file

We'd like to keep `main.rs` to be argument parsing and stuff like that.
4 files changed, 533 insertions(+), 450 deletions(-)

M cra/Cargo.toml
M cra/src/main.rs
M cra/src/output.rs
A cra/src/persistence.rs
M cra/Cargo.toml => cra/Cargo.toml +3 -1
@@ 1,7 1,9 @@
[package]
name = "cra"
version = "0.1.0"
authors = ["Martin Hafskjold Thoresen <git@mht.technology>"]
authors = ["Martin Hafskjold Thoresen <git@mht.wtf>"]
description = "blablabla"

[dependencies]
time = "0.1.42"
clap = "2.33"

M cra/src/main.rs => cra/src/main.rs +84 -447
@@ 1,13 1,25 @@
extern crate clap;
// std
use std::cmp::Ordering::*;
use std::fs::File;
use std::io::{BufRead, BufReader, Write};
use std::io::{BufRead, BufReader, Read, Write};
use std::path::Path;

const R_CUTOFF: f64 = 250.0;

// mods
mod output;
mod persistence;

// mods use
use persistence::Variant::*;
use persistence::{compute_alpha_values_2d, reduce, Persistence, Simplex};


// consts
const R_CUTOFF: f64 = 250.0;
const MAX_DIM: usize = 3;

fn f64_eq(a: f64, b: f64) -> bool {
    // TOOD: this is stupid and should be removed.
    let u = 1_000_000.0;
    let a = (a * u).round() as usize;
    let b = (b * u).round() as usize;


@@ 16,8 28,6 @@ fn f64_eq(a: f64, b: f64) -> bool {
    // (a - b).abs() < 10_000.0 * std::f64::EPSILON
}

const MAX_DIM: usize = 3;

/// All the statistics we want to measure for performance reasoning.
#[derive(Debug)]
pub struct Statistics {


@@ 157,85 167,67 @@ impl Point {
    }
}

/// One simplex
#[derive(Debug, Clone)]
pub struct Simplex {
    /// Original index. We need this for bookkeeping.
    j: usize,
    /// Indices of the faces of this simplex
    faces: Vec<usize>,
    /// The ball radius when this simplex is born.
    r_when_born: f64,
}

impl Simplex {
    /// Get the dimension of this simplex.
    pub fn dim(&self) -> isize {
        self.faces.len() as isize - 1
    }
}

/// Data we use for stuff.
pub struct Persistence {
    /// All simplices in the data set
    simplices: Vec<Simplex>,
    /// Point location data, for rendering.
    points: Vec<[f64; 2]>,
}

impl Persistence {
    /// Fill in some stats
    fn register_stats(&self, stats: &mut Statistics) {
        for s in &self.simplices {
            if s.dim() < 0 {
                continue;
            }
            stats.num_of_dimens[s.dim() as usize] += 1;
        }
    }
}

pub fn read_input_stdin2() -> Result<Persistence, Box<dyn std::error::Error>> {
    let stdin = std::io::stdin();
    let mut lines = BufReader::new(stdin.lock()).lines();

    // File format:
    //     number-of-verts
    //     point data of the vertices
    //     empty line
    //     triangle point indices

    let first_line = lines.next().ok_or("empty stdin")??;
/// Read a point triangulation from `reader`.
///
/// ### File format
/// The file format used is as follows
///
///     number-of-verts
///     point data of the vertices
///     empty line
///     triangle point indices
///
/// ### Note
/// This only reads a 2D triangulation.

pub fn read_triangulation<R: Read>(
    r: &mut R,
) -> Result<(Vec<[f64; 2]>, Vec<Simplex>), Box<dyn std::error::Error>> {
    let mut lines = BufReader::new(r).lines();
    let first_line = lines.next().ok_or("reader is empty")??;
    let number_of_vertices: usize = first_line.parse::<usize>()?;

    let points = (&mut lines)
        .take(number_of_vertices)
        .map(|line| {
            let line = line.unwrap();
            // TODO: Errors here aren't caught.
            let line = line.expect("Need to catch errors in here!");
            let mut buf = [0f64; 2];
            for (i, string) in line.split(' ').enumerate() {
                buf[i] = string.parse::<f64>().unwrap();
                assert!(
                    i < 2,
                    "the line should not contain more than two numbers! (2D only!)"
                );
                buf[i] = string
                    .parse::<f64>()
                    .expect("Need to catch errors in here!");
            }
            buf
        })
        .collect::<Vec<[f64; 2]>>();

    assert_eq!(lines.next().unwrap().unwrap(), ""); // the empty line
    assert_eq!(
        lines
            .next()
            .expect("Missing empty line in triangulation format!")?,
        ""
    ); // the empty line

    // Map the triangle index to the indices of the points it has. This goes directly into the
    // `points` array.
    // Next we read in the face indices.
    let mut faces_to_points = Vec::new();
    for line in lines {
        let mut buf = [0usize; 3];
        for (i, string) in line.unwrap().split(' ').enumerate() {
            assert!(
                i < 3,
                "the line should not contain more than three indices (2D only!)"
            );
            buf[i] = string.parse::<usize>().unwrap();
        }
        faces_to_points.push(buf);
    }

    // Now we have a vec of euclidian point data, and a vec of indices to the point data which
    // forms triangles. Next up we construct the simplices.

    // TOOD: clean this up: it is bad.
    let simplices: Vec<Vec<usize>> = {
        let mut simplices: Vec<Vec<usize>> =
            Vec::with_capacity(points.len() + faces_to_points.len() + 1);


@@ 284,7 276,7 @@ pub fn read_input_stdin2() -> Result<Persistence, Box<dyn std::error::Error>> {
        simplices
    };

    let mut simplices = simplices
    let simplices = simplices
        .into_iter()
        .enumerate()
        .map(|(j, v)| Simplex {


@@ 294,402 286,47 @@ pub fn read_input_stdin2() -> Result<Persistence, Box<dyn std::error::Error>> {
        })
        .collect::<Vec<_>>();

    let simplices_ptr = simplices.as_ptr();
    for (current_simplex_i, simplex) in simplices.iter_mut().enumerate() {
        let simpl = |i: usize| unsafe {
            // We really only require `current_simplex_i != i` to avoid aliasing,
            // but this is a logic error if not.
            assert!(current_simplex_i > i);
            let target: *mut Simplex = simplices_ptr.offset(i as isize) as *mut _;
            &mut *target
        };
        match simplex.dim() {
            -1 | 0 => simplex.r_when_born = 0.0,
            1 => {
                let p0 = Point(points[simplex.faces[0] - 1]);
                let p1 = Point(points[simplex.faces[1] - 1]);
                let diff = p1 - p0;
                let r = 0.5 * diff.len();
                simplex.r_when_born = r;
            }
            2 => {
                let xy_i = simplex.faces[0];
                let xy = simpl(xy_i);
                let yz_i = simplex.faces[1];
                let _yz = simpl(yz_i);
                let zx_i = simplex.faces[2];
                let zx = simpl(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();

                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 ajust the `r_when_born` on
                    // this edge.
                    if am == ax {
                        // debug_assert_eq!(yz / 2.0, simpl(yz_i).r_when_born);
                        simpl(yz_i).r_when_born = circumradius;
                    } else if am == ay {
                        // debug_assert_eq!(xz / 2.0, simpl(zx_i).r_when_born);
                        simpl(zx_i).r_when_born = circumradius;
                    } else {
                        // debug_assert_eq!(xy / 2.0, simpl(xy_i).r_when_born);
                        simpl(xy_i).r_when_born = circumradius;
                    }
                }
                simplex.r_when_born = circumradius;
            }
            more => panic!("Only support 2d simplices! ({} > 2)", more),
        }
    }
    // The times in which a simplex dies is determined by the pairing: it will die when a simplex
    // of dimension d+1 is born.

    Ok(Persistence { simplices, points })
}

/// The result from reducing.
#[derive(Debug, PartialEq)]
pub struct Reduction {
    /// Two simplices are paired when a simplex has been fully reduced; it is then paired with its highest boundary element.
    pairings: Vec<(usize, usize)>,
    ///
    j_stored_at: Vec<usize>,
    /// This maps an index to the reduced simplex which maximum boundary element is that index.
    simplex_with_low: Vec<IndList>,
    /// Here we track which elements are added to which columns.
    /// This was probably used for the output?
    additions: Vec<(usize, usize)>,
}

impl Reduction {
    /// Write out the reduced thing to the `writer`. Useful for debugging.
    #[allow(dead_code)]
    fn write_out<W: Write>(
        &self,
        writer: &mut W,
        simplices: &[Simplex],
    ) -> Result<(), std::io::Error> {
        for (&j, il) in self.j_stored_at.iter().zip(self.simplex_with_low.iter()) {
            let simplex_i = simplices.iter().position(|s| s.j == j).unwrap();
            let simplex = &simplices[simplex_i];
            writeln!(
                writer,
                "j={}, faces={:?} => {:?}",
                simplex.j, simplex.faces, il
            )?;
        }
        Ok(())
    }

    /// Write out a table containing information of the reduction.
    fn write_table<W: Write>(
        &self,
        w: &mut W,
        simplices: &[Simplex],
    ) -> Result<(), std::io::Error> {
        let n = self.simplex_with_low.len();

        let additions_of = {
            let mut v = vec![vec![]; n];
            for &(from, to) in &self.additions {
                v[to].push(from);
            }
            for ve in v.iter_mut() {
                ve.sort();
            }
            v
        };

        let i_of_j = {
            let mut v = vec![None; n];
            for (i, &j) in self.j_stored_at.iter().enumerate() {
                // simplex `j` is stored at index `i` in `self.simplex_with_low`.
                if i > 0 {
                    v[j] = Some(i);
                }
            }
            v
        };

        writeln!(w, "j;| faces;| adds;| reduced")?;
        for j in 1..n {
            let s = &simplices[j];
            write!(w, "{};| {:?};| {:?};| ", s.j, s.faces, additions_of[j])?;
            if let Some(_i) = i_of_j[j] {
                // writeln!(w, "{:?}", self.simplex_with_low[i].get())?;
                writeln!(w, " +")?;
            } else {
                writeln!(w)?;
            }
        }
        Ok(())
    }
}

/// The three variants of the algorithm.
#[derive(PartialEq)]
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,
    /// It's not clear why this is interesting?
    ExhaustiveWhenD1,
    Ok((points, simplices))
}
use Variant::*;

/// 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, stats: &mut Statistics) -> Reduction {
    let t0 = time::precise_time_ns();
    let n = p.simplices.len();
    // `simplex_with_low` maps the highest index `m` to the simplex containing it.
    let mut simplex_with_low: Vec<IndList> = vec![Null; n];
    // `j_stored_at[m]` contains the index of the simplex stored at `simplex_with_low[m]`.
    let mut j_stored_at: Vec<usize> = vec![0; n];

    let mut pairings = Vec::new();
    let mut additions = Vec::new();

	// We're reducing each simplex.
    for simplex in &p.simplices {
        let j = simplex.j;
        // With no boundary we're done.
        if simplex.faces.last().is_none() {
            continue;
        }
        let s_dim = simplex.dim().max(0) as usize;
        // If this simplex is born too late we also skip it.
        if simplex.r_when_born >= R_CUTOFF {
            stats.skip(s_dim);
            continue;
        }
        // TODO: this probably doens't need to be cloned.
        let mut current_simplex = simplex.faces.clone();
        current_simplex.sort();

        let mut r_iter = 0;
        loop {
            r_iter += 1;
            let low: usize = *current_simplex.last().unwrap();
            let cand = &simplex_with_low[low];
            stats.swl();

			if cand.is_list() {
                additions.push((j_stored_at[low], j));
                column_add(&mut current_simplex, cand.get(), s_dim, stats);
                if current_simplex.is_empty() {
                    stats.zeroed(s_dim);
                    break;
                }
			} else {
                debug_assert!(low <= j);
                pairings.push((low, j));

                stats.placed_full(s_dim);
                simplex_with_low[low] = List(current_simplex);
                j_stored_at[low] = j;

                if variant == Exhaustive || (variant == ExhaustiveWhenD1 && simplex.dim() == 1) {
                    // We want to reduce the column `j` with `low(j) = low`. Try to remove all indices in this
                    // column.
                    let mut iter = 0;
                    'search: loop {
                        let this = simplex_with_low[low].get();
                        let list_len = this.len();
                        for (iter_i, face_i) in (0..list_len - 1).rev().enumerate() {
                            let this_index = this[face_i];
                            let other = &simplex_with_low[this_index];
                            stats.swl();
                            if other.is_null() {
                                // nothing to see here
                                continue;
                            }
                            let this_mut: &mut Vec<usize> = unsafe {
                                // We are indexing in twice into the same array, once mutably, but this is fine,
                                // since the indices are not the same.
                                debug_assert!(this_index != low);
                                let ptr: *const IndList =
                                    simplex_with_low.as_ptr().offset(low as isize);
                                let mptr: *mut IndList = ptr as *mut IndList;
                                (&mut *mptr).get_mut()
                            };

                            additions.push((j_stored_at[this_index], j));
                            column_add(this_mut, other.get(), s_dim, stats);
                            iter += 1;
                            stats.ex_search(iter_i + 1);
                            continue 'search;
                        }
                        stats.ex_search(list_len);
                        stats.ex_reductions(iter);
                        break;
                    }
                }
                break;
            }
        }
        stats.reductions(r_iter);
    }

    // Find the essential cycles: these are entries that did not give death, and are `null` in the
    // `simplex_with_low` (in matrix terms, this means that row `i` has no `low`s in it).

    // for &cyc in zeroed.difference(&killed) {
    //     let s = p.simplices.iter().find(|&s| s.j == cyc).unwrap();
    //     eprintln!("[reduce] Essential cycle: {:?}", s);
    // }

    // for (i, ptr) in simplex_with_low.iter().enumerate() {
    //     eprint!("{:2}: ", i);
    //     if ptr.is_null() {
    //         eprintln!("null");
    //     } else if ptr.is_implicit() {
    //         eprintln!("{}", i);
    //     } else {
    //         for i in ptr.get().0.iter() {
    //             eprint!("{} ", i);
    //         }
    //         eprintln!("");
    //     }
    // }

    let t1 = time::precise_time_ns();
    stats.time = t1 - t0;
    Reduction {
        pairings,
        j_stored_at,
        simplex_with_low,
        additions,
    }
}

/// Represent a list of indices for the boundary of a simplex.
/// This is really just an `Option<Vec<usize>>`, but for some reason we decided to make our own wrapper.
#[derive(Debug, Clone, PartialEq, Eq)]
enum IndList {
    Null,
    List(Vec<usize>),
/// Read a 2D triangulation from `stdin` and compute the alpha values for the simplices.
pub fn read_input_stdin2() -> Result<Persistence, Box<dyn std::error::Error>> {
    let mut stdin = std::io::stdin();
    let (points, mut simplices) = read_triangulation(&mut stdin)?;
    compute_alpha_values_2d(&points, &mut simplices);
    Ok(Persistence { simplices, points })
}

use IndList::*;

impl IndList {
    fn is_null(&self) -> bool {
        self == &IndList::Null
    }
fn main() {
    use clap::*;
    let matches = clap_app!(
        myapp =>
        (name: crate_name!())
        (about: crate_description!())
        (version: crate_version!())
        (author: crate_authors!())
        (@arg triangulation: -t --triangulation<FILE> "Input a 2D triangulation in text form.")
        (@arg exhaustive: -e --exhaustive "Run the exhaustive variant of the algorithm.")
        (@arg svg: --svg<FILE> "Draw the compex to an .svg.")
        (@arg alpha: -a --alpha "Ignore the ordering of the simplices given and calculate and use the alpha ordering.")
        (@arg stats: -s --stats "Print out statistics.")
        (@arg graphviz: --graphiviz<FILE> "Use `graphviz` to graph out which simplices are added together.")
    ).get_matches();

    fn is_list(&self) -> bool {
        match self {
            IndList::List(ref _v) => true,
            _ => false,
        }
    if let Some(path) = matches.value_of("triangulation") {
    	// TODO: input triangulation thing here
        println!("{:?}", path);
    }

    #[allow(dead_code)]
    fn list(&self) -> Option<&Vec<usize>> {
        match self {
            List(ref v) => Some(v),
            _ => None,
        }
    }
	if let Some(path) = matches.value_of("svg") {
    	// TODO: `output::svg()` here
	}

    fn get(&self) -> &Vec<usize> {
        match self {
            List(ref v) => v,
            _ => panic!("IndList::get on non-List value!"),
        }
    }

    fn get_mut(&mut self) -> &mut Vec<usize> {
        match self {
            List(ref mut v) => v,
            _ => panic!("IndList::get on non-List value!"),
        }
    }
}

/// Add the `other` column into the `this` column. Since we're in (mod 2) arithmetic, this `xor`s
/// together the two vectors.
fn column_add(
    this: &mut Vec<usize>,
    other: &Vec<usize>,
    simplex_dim: usize,
    stats: &mut Statistics,
) {
    // TODO: hehe
    static mut BUFFER: Option<Vec<usize>> = None;
    let buffer: &mut Vec<usize> = unsafe {
        if BUFFER.is_none() {
            BUFFER = Some(Vec::with_capacity(1_000));
        }
        BUFFER.as_mut().unwrap()
    };
    debug_assert_eq!(buffer.len(), 0);
    stats.col_add(simplex_dim);
    stats.add_sizes(this.len(), other.len());

    let mut our_i = 0;
    let mut their_i = 0;
    let our_last = this.len();
    let their_last = other.len();

    while our_i < our_last && their_i < their_last {
        if this[our_i] < other[their_i] {
            buffer.push(this[our_i]);
            our_i += 1;
        } else if this[our_i] == other[their_i] {
            our_i += 1;
            their_i += 1;
        } else {
            buffer.push(other[their_i]);
            their_i += 1;
        }
    }

    buffer.extend(&other[their_i..]);
    buffer.extend(&this[our_i..]);
    this.clear();
    this.extend(&*buffer);
    buffer.clear();
}
    return;

fn main() {
    let mut args = std::env::args();
    let name = match args.nth(1) {
        Some(name) => name,

M cra/src/output.rs => cra/src/output.rs +3 -2
@@ 5,7 5,8 @@ use std::iter;
use std::path::Path;
use std::process::{Command, Stdio};

use crate::{Persistence, Point, Reduction, Simplex, Statistics, MAX_DIM, R_CUTOFF};
use crate::{Point, Statistics, MAX_DIM, R_CUTOFF};
use persistence::{Persistence, Reduction, Simplex};

/// Format the given `usize` to a `String` with underscores in between each third digit.
///


@@ 43,7 44,7 @@ pub fn svg(persistence: &Persistence, std: &Reduction, exh: &Reduction, out_path
    const DRAW_TEXT: bool = false;
    let mut f = std::io::BufWriter::new(std::fs::File::create(out_path).unwrap());

	// The padding we use for the .svg, so that we know text will be in the image.
    // The padding we use for the .svg, so that we know text will be in the image.
    const PADDING: f64 = 50f64;

    let mut min_x = persistence.points[0][0];

A cra/src/persistence.rs => cra/src/persistence.rs +443 -0
@@ 0,0 1,443 @@
use crate::{Point, Statistics, R_CUTOFF};
use std::io::Write;

/// The three variants of the algorithm.
#[derive(PartialEq)]
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,
    /// It's not clear why this is interesting?
    ExhaustiveWhenD1,
}
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,
}

impl Simplex {
    /// Get the dimension of this simplex.
    pub fn dim(&self) -> isize {
        self.faces.len() as isize - 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]>,
}

impl Persistence {
    /// Fill in some stats
    pub fn register_stats(&self, stats: &mut Statistics) {
        for s in &self.simplices {
            if s.dim() < 0 {
                continue;
            }
            stats.num_of_dimens[s.dim() as usize] += 1;
        }
    }
}

/// The result from reducing.
#[derive(Debug, PartialEq)]
pub struct Reduction {
    /// Two simplices are paired when a simplex has been fully reduced; it is then paired with its highest boundary element.
    pub pairings: Vec<(usize, usize)>,
    ///
    pub j_stored_at: Vec<usize>,
    /// This maps an index to the reduced simplex which maximum boundary element is that index.
    pub simplex_with_low: Vec<IndList>,
    /// Here we track which elements are added to which columns.
    /// This was probably used for the output?
    pub additions: Vec<(usize, usize)>,
}

impl Reduction {
    /// Write out the reduced thing to the `writer`. Useful for debugging.
    #[allow(dead_code)]
    pub fn write_out<W: Write>(
        &self,
        writer: &mut W,
        simplices: &[Simplex],
    ) -> Result<(), std::io::Error> {
        for (&j, il) in self.j_stored_at.iter().zip(self.simplex_with_low.iter()) {
            let simplex_i = simplices.iter().position(|s| s.j == j).unwrap();
            let simplex = &simplices[simplex_i];
            writeln!(
                writer,
                "j={}, faces={:?} => {:?}",
                simplex.j, simplex.faces, il
            )?;
        }
        Ok(())
    }

    /// Write out a table containing information of the reduction.
    pub fn write_table<W: Write>(
        &self,
        w: &mut W,
        simplices: &[Simplex],
    ) -> Result<(), std::io::Error> {
        let n = self.simplex_with_low.len();

        let additions_of = {
            let mut v = vec![vec![]; n];
            for &(from, to) in &self.additions {
                v[to].push(from);
            }
            for ve in v.iter_mut() {
                ve.sort();
            }
            v
        };

        let i_of_j = {
            let mut v = vec![None; n];
            for (i, &j) in self.j_stored_at.iter().enumerate() {
                // simplex `j` is stored at index `i` in `self.simplex_with_low`.
                if i > 0 {
                    v[j] = Some(i);
                }
            }
            v
        };

        writeln!(w, "j;| faces;| adds;| reduced")?;
        for j in 1..n {
            let s = &simplices[j];
            write!(w, "{};| {:?};| {:?};| ", s.j, s.faces, additions_of[j])?;
            if let Some(_i) = i_of_j[j] {
                // writeln!(w, "{:?}", self.simplex_with_low[i].get())?;
                writeln!(w, " +")?;
            } else {
                writeln!(w)?;
            }
        }
        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, stats: &mut Statistics) -> Reduction {
    let t0 = time::precise_time_ns();
    let n = p.simplices.len();
    // `simplex_with_low` maps the highest index `m` to the simplex containing it.
    let mut simplex_with_low: Vec<IndList> = vec![Null; n];
    // `j_stored_at[m]` contains the index of the simplex stored at `simplex_with_low[m]`.
    let mut j_stored_at: Vec<usize> = vec![0; n];

    let mut pairings = Vec::new();
    let mut additions = Vec::new();

    // We're reducing each simplex.
    for simplex in &p.simplices {
        let j = simplex.j;
        // With no boundary we're done.
        if simplex.faces.last().is_none() {
            continue;
        }
        let s_dim = simplex.dim().max(0) as usize;
        // If this simplex is born too late we also skip it.
        if simplex.r_when_born >= R_CUTOFF {
            stats.skip(s_dim);
            continue;
        }
        // TODO: this probably doens't need to be cloned.
        let mut current_simplex = simplex.faces.clone();
        current_simplex.sort();

        let mut r_iter = 0;
        loop {
            r_iter += 1;
            let low: usize = *current_simplex.last().unwrap();
            let cand = &simplex_with_low[low];
            stats.swl();

            if cand.is_list() {
                additions.push((j_stored_at[low], j));
                column_add(&mut current_simplex, cand.get(), s_dim, stats);
                if current_simplex.is_empty() {
                    stats.zeroed(s_dim);
                    break;
                }
            } else {
                debug_assert!(low <= j);
                pairings.push((low, j));

                stats.placed_full(s_dim);
                simplex_with_low[low] = List(current_simplex);
                j_stored_at[low] = j;

                if variant == Exhaustive || (variant == ExhaustiveWhenD1 && simplex.dim() == 1) {
                    // We want to reduce the column `j` with `low(j) = low`. Try to remove all indices in this
                    // column.
                    let mut iter = 0;
                    'search: loop {
                        let this = simplex_with_low[low].get();
                        let list_len = this.len();
                        for (iter_i, face_i) in (0..list_len - 1).rev().enumerate() {
                            let this_index = this[face_i];
                            let other = &simplex_with_low[this_index];
                            stats.swl();
                            if other.is_null() {
                                // nothing to see here
                                continue;
                            }
                            let this_mut: &mut Vec<usize> = unsafe {
                                // We are indexing in twice into the same array, once mutably, but this is fine,
                                // since the indices are not the same.
                                debug_assert!(this_index != low);
                                let ptr: *const IndList =
                                    simplex_with_low.as_ptr().offset(low as isize);
                                let mptr: *mut IndList = ptr as *mut IndList;
                                (&mut *mptr).get_mut()
                            };

                            additions.push((j_stored_at[this_index], j));
                            column_add(this_mut, other.get(), s_dim, stats);
                            iter += 1;
                            stats.ex_search(iter_i + 1);
                            continue 'search;
                        }
                        stats.ex_search(list_len);
                        stats.ex_reductions(iter);
                        break;
                    }
                }
                break;
            }
        }
        stats.reductions(r_iter);
    }

    // Find the essential cycles: these are entries that did not give death, and are `null` in the
    // `simplex_with_low` (in matrix terms, this means that row `i` has no `low`s in it).

    // for &cyc in zeroed.difference(&killed) {
    //     let s = p.simplices.iter().find(|&s| s.j == cyc).unwrap();
    //     eprintln!("[reduce] Essential cycle: {:?}", s);
    // }

    // for (i, ptr) in simplex_with_low.iter().enumerate() {
    //     eprint!("{:2}: ", i);
    //     if ptr.is_null() {
    //         eprintln!("null");
    //     } else if ptr.is_implicit() {
    //         eprintln!("{}", i);
    //     } else {
    //         for i in ptr.get().0.iter() {
    //             eprint!("{} ", i);
    //         }
    //         eprintln!("");
    //     }
    // }

    let t1 = time::precise_time_ns();
    stats.time = t1 - t0;
    Reduction {
        pairings,
        j_stored_at,
        simplex_with_low,
        additions,
    }
}

/// Represent a list of indices for the boundary of a simplex.
/// This is really just an `Option<Vec<usize>>`, but for some reason we decided to make our own wrapper.
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum IndList {
    Null,
    List(Vec<usize>),
}

use self::IndList::*;

impl IndList {
    pub fn is_null(&self) -> bool {
        self == &IndList::Null
    }

    pub fn is_list(&self) -> bool {
        match self {
            IndList::List(ref _v) => true,
            _ => false,
        }
    }

    #[allow(dead_code)]
    pub fn list(&self) -> Option<&Vec<usize>> {
        match self {
            List(ref v) => Some(v),
            _ => None,
        }
    }

    pub fn get(&self) -> &Vec<usize> {
        match self {
            List(ref v) => v,
            _ => panic!("IndList::get on non-List value!"),
        }
    }

    pub fn get_mut(&mut self) -> &mut Vec<usize> {
        match self {
            List(ref mut v) => v,
            _ => panic!("IndList::get on non-List value!"),
        }
    }
}

/// Add the `other` column into the `this` column. Since we're in (mod 2) arithmetic, this `xor`s
/// together the two vectors.
fn column_add(
    this: &mut Vec<usize>,
    other: &Vec<usize>,
    simplex_dim: usize,
    stats: &mut Statistics,
) {
    // TODO: hehe
    static mut BUFFER: Option<Vec<usize>> = None;
    let buffer: &mut Vec<usize> = unsafe {
        if BUFFER.is_none() {
            BUFFER = Some(Vec::with_capacity(1_000));
        }
        BUFFER.as_mut().unwrap()
    };
    debug_assert_eq!(buffer.len(), 0);
    stats.col_add(simplex_dim);
    stats.add_sizes(this.len(), other.len());

    let mut our_i = 0;
    let mut their_i = 0;
    let our_last = this.len();
    let their_last = other.len();

    while our_i < our_last && their_i < their_last {
        if this[our_i] < other[their_i] {
            buffer.push(this[our_i]);
            our_i += 1;
        } else if this[our_i] == other[their_i] {
            our_i += 1;
            their_i += 1;
        } else {
            buffer.push(other[their_i]);
            their_i += 1;
        }
    }

    buffer.extend(&other[their_i..]);
    buffer.extend(&this[our_i..]);
    this.clear();
    this.extend(&*buffer);
    buffer.clear();
}

/// 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: &mut [Simplex]) {
    fn alpha_1d(simplex: &Simplex, points: &[[f64; 2]]) -> f64 {
        let p0 = Point(points[simplex.faces[0] - 1]);
        let p1 = Point(points[simplex.faces[1] - 1]);
        let diff = p1 - p0;
        0.5 * diff.len()
    }

    fn alpha_2d(simplex: &Simplex, points: &[[f64; 2]], simplices: &[Simplex]) -> f64 {
        let simpl = |i| unsafe { simplex_mut(i, simplices) };

        let xy_i = simplex.faces[0];
        let xy = simpl(xy_i);
        let yz_i = simplex.faces[1];
        let _yz = simpl(yz_i);
        let zx_i = simplex.faces[2];
        let zx = simpl(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();

        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 {
                // debug_assert_eq!(yz / 2.0, simpl(yz_i).r_when_born);
                simpl(yz_i).r_when_born = circumradius;
            } else if am == ay {
                // debug_assert_eq!(xz / 2.0, simpl(zx_i).r_when_born);
                simpl(zx_i).r_when_born = circumradius;
            } else {
                // debug_assert_eq!(xy / 2.0, simpl(xy_i).r_when_born);
                simpl(xy_i).r_when_born = circumradius;
            }
        }
        circumradius
    }

    unsafe fn simplex_mut(i: usize, simplices: &[Simplex]) -> &mut Simplex {
        let simplices_ptr = simplices.as_ptr();
        let target: *mut Simplex = simplices_ptr.offset(i as isize) as *mut _;
        &mut *target
    }

    for current_simplex_i in 0..simplices.len() {
        let (former, latter) = simplices.split_at_mut(current_simplex_i);
        let simplex = &mut latter[0];
        match simplex.dim() {
            -1 | 0 => simplex.r_when_born = 0.0,
            1 => {
                simplex.r_when_born = alpha_1d(&simplex, &points);
            }
            2 => {
                simplex.r_when_born = alpha_2d(&simplex, &points, former);
            }
            more => panic!("Only support 2d simplices! ({} > 2)", more),
        }
    }
}