~mht/cra

deef01a7309fbc9e6f49ec728ed7888f9bbb2c1f — Martin Hafskjold Thoresen 1 year, 3 months ago fc08a98
Did a lot of things.

This is not tested at all, since CGAL is stupid.
4 files changed, 301 insertions(+), 270 deletions(-)

M cra/src/main.rs
M cra/src/output.rs
M cra/src/persistence.rs
M cra/triangulate2.cpp
M cra/src/main.rs => cra/src/main.rs +115 -53
@@ 164,6 164,84 @@ impl Point {
    }
}

/// Read a boundary file. The first line should be the number of vertices,
/// and consequent lines should be whitespace separated indices for the boundary
/// of that simplex, using 1-indexing.
/// Everything after a '#' in a line is ignored, and empty lines are ignored.
///
/// For now, we only allow the boundary to be the minimal possible for that
/// dimension. Eg, the dimension of a simplex is the size of its boundary minus 1.
///
/// ## Example
///
/// The following file would define the boundary of a tetrahedron.
///
/// 	4			# 4 vertices
/// 	1 2			# simplex #5 is the first edge
/// 	1 3			#  6
/// 	1 4			#  7
/// 	2 3			#  8
/// 	2 4			#  9
/// 	3 4			# 10
/// 	5 6 8   	# 11 - the first triangle 123
/// 	5 7 9		# 12 (124)
/// 	6 7 10  	# 13 (134)
/// 	8 9 10      # 14 (234)
///     11 12 13 14 # the tetrahedron
///
pub fn read_boundary<R: Read>(r: R) -> Result<Vec<Simplex>, Box<dyn std::error::Error>> {
    let mut simplices = vec![Simplex::empty()];

    let mut lines = BufReader::new(r).lines();
    let num_vertices = lines
        .next()
        .ok_or("The first line should contain the number of vertices, but was empty")??
        .parse::<usize>()?;
    for i in 0..num_vertices {
        simplices.push(Simplex {
            j: i,
            faces: vec![0],
        });
    }

    let mut j = simplices.len();
    for line in lines.flat_map(|id| id) {
        let line: &str = if let Some(i) = line.find('#') {
            &line[..i]
        } else {
            &line
        }
        .trim();

        if line.is_empty() {
            continue;
        }

        let mut v = Vec::new();
        for f in line.split_whitespace() {
            v.push(f.parse::<usize>()?);
        }
        simplices.push(Simplex { j, faces: v });
        j += 1;
    }

    Ok(simplices)
}

/// Read vertex 2D positions, one for each line.
pub fn read_2d_points<R: Read>(r: R) -> Result<Vec<[f64; 2]>, Box<dyn std::error::Error>> {
    let lines = BufReader::new(r).lines();
    let mut points = Vec::new();
    for line in lines {
        let line = line?;
        let mut sp = line.split_whitespace();
        let x = sp.next().ok_or("Missing x in line")?.parse::<f64>()?;
        let y = sp.next().ok_or("Missing x in line")?.parse::<f64>()?;
        points.push([x, y]);
    }
    Ok(points)
}

/// Read a point triangulation from `reader`.
///
/// ### File format


@@ 284,6 362,17 @@ pub fn read_triangulation<R: Read>(
    Ok((points, simplices))
}

/// Open the file given, or if the filename is "-", open `stdin`.
fn open_file_or_stdin(
    filename: &str,
) -> Result<Box<dyn std::io::Read>, Box<dyn std::error::Error>> {
    if filename == "-" {
        Ok(Box::new(std::io::stdin()))
    } else {
        Ok(Box::new(File::open(Path::new(filename))?))
    }
}

fn main() {
    use clap::*;
    let matches = clap_app!(


@@ 303,7 392,7 @@ fn main() {
        (@arg points: --points[FILE] "Input the location data for the 1-simplices.")
        (@arg debug: -d "print debug information to `stdout`")
    ).arg(Arg::from_usage(
            "[boundary] -b, --boundary <FILE> 'The boundary of all simplices in the complex.'",
            "<boundary> -b, --boundary <FILE> 'The boundary of all simplices in the complex.'",
        ).long_help("\
This file should contain the boundaries of your simplices.
The boundary of the empty complex and the vertices are implied, and so


@@ 311,63 400,36 @@ the first boundaries in the list should be that of the 1-simplices.
"),
    ).get_matches();

    let mut persistence = if let Some(path) = matches.value_of("triangulation") {
        // TODO: input triangulation thing here
        let mut reader: Box<dyn std::io::Read> = if path == "-" {
            Box::new(std::io::stdin())
        } else {
            Box::new(File::open(path).expect("No such file"))
        };
        let (points, simplices) =
            read_triangulation(&mut reader).expect("failed to read triangulation");
        let (simplices, alphas) = compute_alpha_values_2d(&points, &simplices);
        Persistence {
            simplices,
            points,
            alphas,
        }
    } else {
        // Need to take some other format, then!
        unimplemented!();
    let debug = matches.is_present("debug");

    let mut simplices = {
        let mut reader = open_file_or_stdin(matches.value_of("boundary").unwrap())
            .expect("failed to open given file");
        read_boundary(&mut reader).expect("Couldn't read boundary")
    };

    // TODO: this should be done somewhere else!
    {
        // Map j to sorted index.
        let sorted_index_of_j = {
            let mut v = (0..persistence.simplices.len()).collect::<Vec<_>>();
            for (i, s) in persistence.simplices.iter().enumerate() {
                v[s.j] = i;
            }
            v
        };
    let mut alphas = None;
    let mut points = None;

        // Change all `r-values` to be in the sorted format, such that `simplex[a].j == a`.
        for (j, s) in persistence.simplices.iter_mut().enumerate() {
            for face in s.faces.iter_mut() {
                *face = sorted_index_of_j[*face];
            }
            s.j = j;
        }
    if matches.is_present("alpha") {
        let pts = {
            let mut reader = open_file_or_stdin(matches.value_of("points").unwrap())
                .expect("failed to open given file");
            read_2d_points(&mut reader).expect("Couldn't read 2d points")
        };

        // Since our sorting is not really numerically stable at all, we still need this assertion.
        for (j, s) in persistence.simplices.iter().enumerate() {
            if let Some(&face) = s.faces.iter().max() {
                if face >= j {
                    eprintln!("face is after simplex! {} >= {}", face, j);
                    eprintln!("{:?}", persistence.simplices[face]);
                    eprintln!("{:?}", persistence.simplices[j]);
                    eprintln!(
                        "{:#?}",
                        &persistence.simplices
                            [(j.saturating_sub(3)..((face + 4).min(persistence.simplices.len())))]
                    );
                    panic!("The ordering is not right!");
                }
            }
        }
        let (ordered_simplices, alphs) = compute_alpha_values_2d(&pts, &simplices);
        alphas = Some(alphs);
        points = Some(pts);
        simplices = ordered_simplices;
    }

    let mut persistence = Persistence {
        simplices,
        points: points.unwrap_or(Vec::new()),
        alphas: alphas.unwrap_or(Vec::new()),
    };

    let mut e_stats = Statistics::new();
    let mut r_stats = Statistics::new();



@@ 378,7 440,7 @@ the first boundaries in the list should be that of the 1-simplices.
    let exhaustive_reduction = reduce(&persistence, Exhaustive, &mut e_stats);
    assert_eq!(standard_reduction.pairings, exhaustive_reduction.pairings);

    if false {
    if debug {
        let mut out = std::io::stderr();
        standard_reduction
            .write_table(&mut out, &persistence.simplices)

M cra/src/output.rs => cra/src/output.rs +105 -121
@@ 42,11 42,12 @@ fn thousand_split(mut num: usize) -> String {
///
pub fn svg(persistence: &Persistence, std: &Reduction, exh: &Reduction, out_path: &Path) {
    // Whether to draw text on all simplices outputted.
    const DRAW_TEXT: bool = false;
    const DRAW_TEXT: bool = true;
    // The padding we use for the .svg, so that we know text will be in the image.
    const PADDING: f64 = 50f64;

    let mut f = std::io::BufWriter::new(std::fs::File::create(out_path).unwrap());
    let mut _f = std::io::BufWriter::new(std::fs::File::create(out_path).unwrap());
    let f = &mut _f;

    let min_x = persistence
        .points


@@ 79,6 80,7 @@ pub fn svg(persistence: &Persistence, std: &Reduction, exh: &Reduction, out_path

    let xo = PADDING - min_x;
    let yo = PADDING - min_y;
    let off = Point([xo, yo]);

    let alpha = &persistence.alphas;



@@ 110,76 112,64 @@ pub fn svg(persistence: &Persistence, std: &Reduction, exh: &Reduction, out_path
            let p1 = Point(persistence.points[y_i - 1]);
            let p2 = Point(persistence.points[z_i - 1]);

            let draw_line = |f: &mut std::io::BufWriter<File>, x0, y0, x1, y1| {
                writeln!(
                    f,
                    r#"<line x1="{}" y1="{}" x2="{}" y2="{}"
                                 style="stroke:black;stroke-width:1;" />"#,
                    x0 + xo,
                    y0 + yo,
                    x1 + xo,
                    y1 + yo,
                )
            };

            let draw_i = |f: &mut std::io::BufWriter<File>, x, y, i| {
                writeln!(
                    f,
                    "<text x=\"{}\" y=\"{}\" style=\"fill:green;font-size:8px\">{}</text>",
                    x + xo,
                    y + yo,
                    i
                )
            };

            let draw_if = |f: &mut std::io::BufWriter<File>, x, y, i, ff| {
                writeln!(
                    f,
                    "<text x=\"{}\" y=\"{}\" style=\"fill:green;font-size:8px\">{} ({:4.1})</text>",
                    x + xo,
                    y + yo,
                    i,
                    ff
                )
            };

            // See whether we can fill in the triangle or not.
            if alpha[simplex_i] >= R_CUTOFF {
                // Cannot fill triangle. Draw edges, maybe?
                if alpha[xy_i] < R_CUTOFF {
                    write!(
                        f,
                        r#"<line x1="{}" y1="{}" x2="{}" y2="{}"
                                 style="stroke:black;stroke-width:1;" />\n"#,
                        p0.0[0] + xo,
                        p0.0[1] + yo,
                        p1.0[0] + xo,
                        p1.0[1] + yo
                    )
                    .unwrap();
                    draw_line(f, p0.0[0], p0.0[1], p1.0[0], p1.0[1]).unwrap();
                    // with text
                    if DRAW_TEXT {
                        let xy_middle = (p1 - p0) / 2.0 + p0;
                        writeln!(
                            f,
                            "<text x=\"{}\" y=\"{}\" style=\"fill:red;font-size:8px\">{}</text>",
                            xy_middle.0[0] + xo,
                            xy_middle.0[1] + yo,
                            xy.j
                        )
                        .unwrap();
                        draw_i(f, xy_middle.0[0], xy_middle.0[1], xy.j).unwrap();
                    }
                }

                if alpha[zx_i] < R_CUTOFF {
                    write!(
                        f,
                        r#"<line x1="{}" y1="{}" x2="{}" y2="{}"
                                 style="stroke:black;stroke-width:1;" />\n"#,
                        p0.0[0] + xo,
                        p0.0[1] + yo,
                        p2.0[0] + xo,
                        p2.0[1] + yo
                    )
                    .unwrap();
                    draw_line(f, p0.0[0], p0.0[1], p2.0[0], p2.0[1]).unwrap();
                    if DRAW_TEXT {
                        let zx_middle = (p0 - p2) / 2.0 + p2;
                        writeln!(
                            f,
                            "<text x=\"{}\" y=\"{}\" style=\"fill:red;font-size:8px\">{}</text>",
                            zx_middle.0[0] + xo,
                            zx_middle.0[1] + yo,
                            zx.j
                        )
                        .unwrap();
                        draw_i(f, zx_middle.0[0], zx_middle.0[1], zx.j).unwrap();
                    }
                }

                if alpha[yz_i] < R_CUTOFF {
                    write!(
                        f,
                        r#"<line x1="{}" y1="{}" x2="{}" y2="{}"
                                 style="stroke:black;stroke-width:1;" />\n"#,
                        p1.0[0] + xo,
                        p1.0[1] + yo,
                        p2.0[0] + xo,
                        p2.0[1] + yo
                    )
                    .unwrap();
                    draw_line(f, p1.0[0], p1.0[1], p2.0[0], p2.0[1]).unwrap();
                    if DRAW_TEXT {
                        let yz_middle = (p2 - p1) / 2.0 + p1;
                        writeln!(
                            f,
                            "<text x=\"{}\" y=\"{}\" style=\"fill:red;font-size:8px\">{}</text>",
                            yz_middle.0[0] + xo,
                            yz_middle.0[1] + yo,
                            yz.j
                        )
                        .unwrap();
                        draw_i(f, yz_middle.0[0], yz_middle.0[1], yz.j).unwrap();
                    }
                }
            } else {


@@ 192,47 182,41 @@ pub fn svg(persistence: &Persistence, std: &Reduction, exh: &Reduction, out_path
                    p2.0[0] + xo,
                    p2.0[1] + yo
                ).unwrap();
                let center = (p0 + p1 + p2) / 3.0 + off;

                writeln!(
                    f,
                    "<text x=\"{}\" y=\"{}\" style=\"fill:black;font-size:8px\">{} ({:4.1})</text>",
                    center.0[0], center.0[1], simplex_i, alpha[simplex_i]
                )
                .unwrap();

                // Draw all edges unconditionally.

                if DRAW_TEXT {
                    let xy_middle = (p1 - p0) / 2.0 + p0;
                    writeln!(
                        f,
                        "<text x=\"{}\" y=\"{}\" style=\"fill:red;font-size:8px\">{}</text>",
                        xy_middle.0[0] + xo,
                        xy_middle.0[1] + yo,
                        xy.j
                    )
                    .unwrap();
                    draw_if(f, xy_middle.0[0], xy_middle.0[1], xy.j, alpha[xy_i]).unwrap();
                    let yz_middle = (p2 - p1) / 2.0 + p1;
                    writeln!(
                        f,
                        "<text x=\"{}\" y=\"{}\" style=\"fill:red;font-size:8px\">{}</text>",
                        yz_middle.0[0] + xo,
                        yz_middle.0[1] + yo,
                        yz.j
                    )
                    .unwrap();
                    draw_if(f, yz_middle.0[0], yz_middle.0[1], yz.j, alpha[yz_i]).unwrap();
                    let zx_middle = (p0 - p2) / 2.0 + p2;
                    writeln!(
                        f,
                        "<text x=\"{}\" y=\"{}\" style=\"fill:red;font-size:8px\">{}</text>",
                        zx_middle.0[0] + xo,
                        zx_middle.0[1] + yo,
                        zx.j
                    )
                    .unwrap();
                    draw_if(f, zx_middle.0[0], zx_middle.0[1], zx.j, alpha[zx_i]).unwrap();
                }
            }
        } else if simplex.dim() == 1 {
            // Edges
        } else if simplex.dim() == 0 {
            // points
        }
    }

    // TODO: what is this?
    let n = persistence.simplices.len();
    let mut std_jumps = vec![0; n];
    let mut exh_jumps = vec![0; n];
    for i in 0..n {
        if std.simplex_with_low[i].is_list() {
            let s = std.simplex_with_low[i].get();
            let e = exh.simplex_with_low[i].get();
        if std.simplex_with_low[i].is_some() {
            let s = std.simplex_with_low[i].as_ref().unwrap();
            let e = exh.simplex_with_low[i].as_ref().unwrap();
            if s.len() == 2 {
                std_jumps[i] = s[0];
                exh_jumps[i] = e[0];


@@ 299,12 283,12 @@ plot '{}' using 1:2 pt 7 ps 0.25 title "pairs", x"#,
//     }
//     freq
// }
// 
//
// pub fn histogram(data: &[usize], label: &str, out_file: &Path) {
//     if data.len() == 0 {
//         panic!("output_histogram: `data` cannot be empty!");
//     }
// 
//
//     let mut frequency = make_histogram(data);
//     // Remove single outliers, and don't give out trailing zeroes to gnuplot.
//     let mut last_i = 0;


@@ 316,13 300,13 @@ plot '{}' using 1:2 pt 7 ps 0.25 title "pairs", x"#,
//             last_i = i;
//         }
//     }
// 
//
//     let file_path = Path::new("tmp-histogram");
//     let mut file = File::create(&file_path).unwrap();
//     for (ind, freq) in frequency.iter().enumerate().take(last_i + 1) {
//         write!(file, "{} {}\n", ind, freq).unwrap();
//     }
// 
//
//     let plot_script = format!(
//         r#"
// stats '{file_path}' nooutput


@@ 341,16 325,16 @@ plot '{}' using 1:2 pt 7 ps 0.25 title "pairs", x"#,
//         label,
//         file_path = file_path.display()
//     );
// 
//
//     run_gnuplot(&plot_script, out_file);
//     std::fs::remove_file(file_path).unwrap();
// }
// 
//
// pub fn histogram2(data_a: &[usize], data_b: &[usize], label: &str, out_file: &Path) {
//     if data_a.len() == 0 || data_b.len() == 0 {
//         panic!("output_histogram: `data` cannot be empty!");
//     }
// 
//
//     let mut a_hist = make_histogram(data_a);
//     let mut b_hist = make_histogram(data_b);
//     if a_hist.len() < b_hist.len() {


@@ 362,10 346,10 @@ plot '{}' using 1:2 pt 7 ps 0.25 title "pairs", x"#,
//             b_hist.push(0);
//         }
//     }
// 
//
//     let file_path = Path::new("tmp-2histogram");
//     let mut file = File::create(&file_path).unwrap();
// 
//
//     for (a, b) in a_hist
//         .iter()
//         .zip(b_hist.iter())


@@ 373,7 357,7 @@ plot '{}' using 1:2 pt 7 ps 0.25 title "pairs", x"#,
//     {
//         write!(file, "{} {}\n", a, b).unwrap();
//     }
// 
//
//     let plot_script = format!(
//         r#"
// stats '{file_path}' nooutput


@@ 400,7 384,7 @@ plot '{}' using 1:2 pt 7 ps 0.25 title "pairs", x"#,
//         label,
//         file_path = file_path.display()
//     );
// 
//
//     run_gnuplot(&plot_script, out_file);
//     std::fs::remove_file(file_path).unwrap();
// }


@@ 602,6 586,30 @@ pub fn tex_table<W: Write>(
    Ok(())
}

// Get the three vertices of a triangle. This isn't super trivial, since the triangle has
// three edges, and so we have to eliminate the duplicates. Still pretty straight forward.
fn get_triangle_xyz(simplex: &Simplex, simplices: &[Simplex]) -> (usize, usize, usize) {
    let xy_i = simplex.faces[0];
    let xy = &simplices[xy_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.
    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])
    }
}

/// Write out a dotfile with the directed edges given by `simplex_with_low`
pub fn graphviz<W: Write>(
    p: &Persistence,


@@ 609,30 617,6 @@ pub fn graphviz<W: Write>(
    exh: &Reduction,
    w: &mut W,
) -> Result<(), std::io::Error> {
    // Get the three vertices of a triangle. This isn't super trivial, since the triangle has
    // three edges, and so we have to eliminate the duplicates. Still pretty straight forward.
    fn get_triangle_xyz(simplex: &Simplex, simplices: &[Simplex]) -> (usize, usize, usize) {
        let xy_i = simplex.faces[0];
        let xy = &simplices[xy_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.
        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])
        }
    }

    const SCALE: f64 = 5.0;
    writeln!(w, "digraph {{")?;
    writeln!(w, "  rankdir=LR;")?;


@@ 681,9 665,9 @@ pub fn graphviz<W: Write>(

    // TODO: Rethink what we want here
    for i in 0..p.simplices.len() {
        if std.simplex_with_low[i].is_list() {
            let s = std.simplex_with_low[i].get();
            let e = exh.simplex_with_low[i].get();
        if std.simplex_with_low[i].is_some() {
            let s = std.simplex_with_low[i].as_ref().unwrap();
            let e = exh.simplex_with_low[i].as_ref().unwrap();
            assert_eq!(std.j_stored_at[i], exh.j_stored_at[i]);
            if s.len() == 2 {
                writeln!(w, "  {} -> {} [color=\"red\"]", s[1], s[0])?;

M cra/src/persistence.rs => cra/src/persistence.rs +47 -68
@@ 31,6 31,13 @@ impl Simplex {
    pub fn dim(&self) -> isize {
        self.faces.len() as isize - 1
    }
    /// Make the empty simplex
    pub fn empty() -> Self {
        Simplex {
            j: 0,
            faces: Vec::new(),
        }
    }
}

/// Data we use for stuff.


@@ 64,7 71,7 @@ pub struct Reduction {
    ///
    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>,
    pub simplex_with_low: Vec<Option<Vec<usize>>>,
    /// Here we track which elements are added to which columns.
    /// This was probably used for the output?
    pub additions: Vec<(usize, usize)>,


@@ 142,7 149,7 @@ pub fn reduce(p: &Persistence, variant: Variant, stats: &mut Statistics) -> Redu
    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];
    let mut simplex_with_low: Vec<Option<Vec<usize>>> = vec![None; 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];



@@ 173,9 180,9 @@ pub fn reduce(p: &Persistence, variant: Variant, stats: &mut Statistics) -> Redu
            let cand = &simplex_with_low[low];
            stats.swl();

            if cand.is_list() {
            if cand.is_some() {
                additions.push((j_stored_at[low], j));
                column_add(&mut current_simplex, cand.get(), s_dim, stats);
                column_add(&mut current_simplex, cand.as_ref().unwrap(), s_dim, stats);
                if current_simplex.is_empty() {
                    stats.zeroed(s_dim);
                    break;


@@ 185,7 192,7 @@ pub fn reduce(p: &Persistence, variant: Variant, stats: &mut Statistics) -> Redu
                pairings.push((low, j));

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

                if variant == Exhaustive || (variant == ExhaustiveWhenD1 && simplex.dim() == 1) {


@@ 193,13 200,13 @@ pub fn reduce(p: &Persistence, variant: Variant, stats: &mut Statistics) -> Redu
                    // column.
                    let mut iter = 0;
                    'search: loop {
                        let this = simplex_with_low[low].get();
                        let this = simplex_with_low[low].as_ref().unwrap();
                        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() {
                            if other.is_none() {
                                // nothing to see here
                                continue;
                            }


@@ 207,14 214,14 @@ pub fn reduce(p: &Persistence, variant: Variant, stats: &mut Statistics) -> Redu
                                // 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 =
                                let ptr: *const Option<Vec<usize>> =
                                    simplex_with_low.as_ptr().offset(low as isize);
                                let mptr: *mut IndList = ptr as *mut IndList;
                                (&mut *mptr).get_mut()
                                let mptr = ptr as *mut Option<Vec<usize>>;
                                (*mptr).as_mut().unwrap()
                            };

                            additions.push((j_stored_at[this_index], j));
                            column_add(this_mut, other.get(), s_dim, stats);
                            column_add(this_mut, other.as_ref().unwrap(), s_dim, stats);
                            iter += 1;
                            stats.ex_search(iter_i + 1);
                            continue 'search;


@@ 262,51 269,6 @@ pub fn reduce(p: &Persistence, variant: Variant, stats: &mut Statistics) -> Redu
    }
}

/// 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(


@@ 359,8 321,8 @@ pub fn compute_alpha_values_2d(
    simplices: &[Simplex],
) -> (Vec<Simplex>, Vec<f64>) {
    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 p0 = Point(points[simplex.faces[0]]);
        let p1 = Point(points[simplex.faces[1]]);
        let diff = p1 - p0;
        0.5 * diff.len()
    }


@@ 410,9 372,8 @@ pub fn compute_alpha_values_2d(
        let circumradius = (xy * yz * xz)
            / ((xy + yz + xz) * (xy + yz - xz) * (xy - yz + xz) * (-xy + yz + xz)).sqrt();

		let mut edit = None;
		// TODO: check this out again
        if am * 180.0 / std::f64::consts::PI >= 120.0 {
        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


@@ 428,14 389,14 @@ pub fn compute_alpha_values_2d(
        (circumradius, edit)
    }

	#[derive(PartialEq)]
    #[derive(PartialEq)]
    struct CmpSimplex {
        dim: usize,
        alpha: f64,
        index: usize,
    }

    impl Eq for CmpSimplex { }
    impl Eq for CmpSimplex {}

    impl PartialOrd for CmpSimplex {
        fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {


@@ 445,12 406,15 @@ pub fn compute_alpha_values_2d(

    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!"))
            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 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];


@@ 463,13 427,13 @@ pub fn compute_alpha_values_2d(
                    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
            index: current_simplex_i,
        });
    }



@@ 485,5 449,20 @@ pub fn compute_alpha_values_2d(
        }
    }

    // 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)
}

M cra/triangulate2.cpp => cra/triangulate2.cpp +34 -28
@@ 1,53 1,59 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

#include <vector>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Delaunay_triangulation_2<K>      Triangulation;
typedef Triangulation::Point          Point;
typedef CGAL::Delaunay_triangulation_2<K> Triangulation;
typedef Triangulation::Point Point;
typedef Triangulation::Vertex Vertex;
typedef Triangulation::Edge Edge;
typedef Triangulation::Face Face;

void triangulate_stdin() {
  // Read in the points into a vec
  float x, y;
  std::vector<Point> points;
  while (1) {
    if (!(std::cin >> x)) break;
    if (!(std::cin >> y)) break;
    if (!(std::cin >> x))
      break;
    if (!(std::cin >> y))
      break;
    points.push_back(Point(x, y));
  }

  // Output number of points
  printf("%zu\n", points.size());
  // output all points
  for (auto p : points) {
    printf("%f %f\n", p.x(), p.y());
  }

  // blank line divisor
  printf("\n");

  std::map<Point, int> index_of_vertex;
  int i = 0;
  for (Point p : points) {
    index_of_vertex[p] = i++;
  }
  
  // here you go, cgal
  Triangulation triangulation;
  triangulation.insert(points.begin(), points.end());

  for (auto fi = triangulation.finite_faces_begin(); 
       fi != triangulation.finite_faces_end();
       fi++) {
    auto p0 = fi->vertex(0)->point(),
         p1 = fi->vertex(1)->point(),
         p2 = fi->vertex(2)->point();
    auto i0 = index_of_vertex[p0],
         i1 = index_of_vertex[p1],
         i2 = index_of_vertex[p2];
    printf("%d %d %d\n", i0, i1, i2);
  int j = 0;
  // std::map<Vertex, int> vert_to_j;
  for (auto fi = triangulation.finite_vertices_begin();
       fi != triangulation.finite_vertices_end(); fi++) {
    // vert_to_j[*fi] = ++j;
    // Skip outputting the points, since they all only have the empty set.
  }

  // std::map<Edge, int> edge_to_j;
  for (auto fi = triangulation.finite_edges_begin();
       fi != triangulation.finite_edges_end(); fi++) {
    printf("edge %d\n", fi->second);
    // edge_to_j[*fi] = ++j;
    // int v0 = fi->first->vertex(0);
	// int v1 = fi->second;
    // auto i0 = vert_to_j[v0], i1 = vert_to_j[v1];
    // printf("%d %d\n", i0, i1);
  }

//   for (auto fi = triangulation.finite_faces_begin();
//        fi != triangulation.finite_faces_end(); fi++) {
//     auto e0 = fi->edge(0), e1 = fi->edge(1), e2 = fi->edge(2);
//     auto i0 = edge_to_j[e0], i1 = edge_to_j[e1], i2 = edge_to_j[e2];
//     printf("%d %d %d\n", i0, i1, i2);
//   }
}

int main() {