~mht/cra

92ecc9ca15537a2719f06e58368c89f4249ff7cd — Martin Hafskjold Thoresen 1 year, 3 months ago bd8a9e3
Move output things to `output.rs`, and fix warnings
3 files changed, 761 insertions(+), 741 deletions(-)

D cra/rust-toolchain
M cra/src/main.rs
A cra/src/output.rs
D cra/rust-toolchain => cra/rust-toolchain +0 -1
@@ 1,1 0,0 @@
nightly

M cra/src/main.rs => cra/src/main.rs +18 -740
@@ 1,15 1,12 @@
#![feature(nll)]

use std::cmp::Ordering::*;
use std::fmt::Write as fmtWrite;
use std::fs::File;
use std::io::{BufRead, BufReader, Write};
use std::iter;
use std::path::Path;
use std::process::{Command, Stdio};

const R_CUTOFF: f64 = 250.0;

mod output;

fn f64_eq(a: f64, b: f64) -> bool {
    let u = 1_000_000.0;
    let a = (a * u).round() as usize;


@@ 119,121 116,6 @@ impl Statistics {
    pub fn swl(&mut self) {
        self.simplex_with_low_queries += 1;
    }

    pub fn write_stats<W: Write>(&self, writer: &mut W) -> Result<(), std::io::Error> {
        let mut print_pairs = Vec::new();

        print_pairs.push(("Column additions", thousand_split(self.col_adds)));

        let adds = self.add_size_sum.iter().sum::<usize>() as usize;
        print_pairs.push(("Estimate total 'column_add' cost", thousand_split(adds)));

        print_pairs.push((
            "`simplex_with_low` queries",
            thousand_split(self.simplex_with_low_queries),
        ));

        print_pairs.push(("Time spent (ns)", thousand_split(self.time as usize)));

        // ~10 cycles per L2 hit seems reasonable?
        let instr_cost_estimate = self.simplex_with_low_queries * 10 + adds;
        print_pairs.push((
            "Instruction cost estimate",
            thousand_split(instr_cost_estimate),
        ));

        let add_size0 =
            self.add_size.iter().map(|(a, _)| a).sum::<usize>() as f64 / self.add_size.len() as f64;
        print_pairs.push(("Average size of first operand", format!("{}", add_size0)));

        let add_size1 =
            self.add_size.iter().map(|(_, b)| b).sum::<usize>() as f64 / self.add_size.len() as f64;
        print_pairs.push(("Average size of second operand", format!("{}", add_size1)));

        let num_iters = self.num_iters.iter().sum::<usize>() as f64 / self.num_iters.len() as f64;
        print_pairs.push((
            "Average #iterations before reducing a column",
            format!("{}", num_iters),
        ));

        let column_add_cost =
            self.add_size_sum.iter().sum::<usize>() as f64 / self.add_size_sum.len() as f64;
        print_pairs.push((
            "Average cost of one column addition",
            format!("{}", column_add_cost),
        ));

        if self.ex_reductions.len() > 0 {
            let ex_reductions =
                self.ex_reductions.iter().sum::<usize>() as f64 / self.ex_reductions.len() as f64;
            print_pairs.push((
                "Average iterations to exh. reduce a col",
                format!("{}", ex_reductions),
            ));
        }

        if self.ex_searches.len() > 0 {
            let ex_searches =
                self.ex_searches.iter().sum::<usize>() as f64 / self.ex_searches.len() as f64;
            print_pairs.push((
                "Average searches to find `j=low(k)`",
                format!("{}", ex_searches),
            ));
        }

        print_pairs.push(("Number of zeroed columns", thousand_split(self.zeroed)));
        print_pairs.push((
            "Number of single columns",
            thousand_split(self.placed_single),
        ));
        print_pairs.push((
            "Number of remaining columns",
            thousand_split(self.placed_full),
        ));

        let mut last_pos = 0;
        for (i, &n) in self.num_of_dimens.iter().enumerate() {
            if n > 0 {
                last_pos = i;
            }
        }
        let strs = (0..last_pos + 1)
            .map(|d| format!("Simplices with d={}", d))
            .collect::<Vec<_>>();
        for d in 0..last_pos + 1 {
            print_pairs.push((&strs[d], thousand_split(self.num_of_dimens[d])));
        }

        let max_str = print_pairs.iter().map(|(s, _)| s.len()).max().unwrap();
        for (left, right) in print_pairs {
            write!(writer, "{:>2$}: {}\n", left, right, max_str)?;
        }
        Ok(())
    }
}

fn thousand_split(mut num: usize) -> String {
    let mut out = String::new();
    let mut nums = vec![];
    while num > 0 {
        nums.push(num % 1_000);
        num /= 1_000;
    }
    for i in (1..nums.len()).rev() {
        if i == nums.len() - 1 {
            write!(out, "{},", nums[i]).unwrap();
        } else {
            write!(out, "{:03},", nums[i]).unwrap();
        }
    }
    if nums.len() == 0 {
        write!(out, "0").unwrap();
    } else if nums.len() == 1 {
        write!(out, "{}", nums[0]).unwrap();
    } else {
        write!(out, "{:03}", nums[0]).unwrap();
    }
    out
}

#[derive(Copy, Clone, Debug)]


@@ 313,7 195,7 @@ impl Persistence {
    }
}

pub fn read_input_stdin2() -> Result<Persistence, Box<std::error::Error>> {
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();



@@ 508,6 390,8 @@ pub struct Reduction {
}

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();


@@ 548,7 432,7 @@ impl Reduction {
            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, " +");
                writeln!(w, " +")?;
            } else {
                writeln!(w)?;
            }


@@ 708,6 592,7 @@ impl IndList {
        }
    }

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


@@ 770,613 655,6 @@ fn column_add(this: &mut Vec<usize>, other: &Vec<usize>, simplex_dim: usize, sta
    buffer.clear();
}

fn output_svg(persistence: &Persistence, std: &Reduction, exh: &Reduction, out_path: &Path) {
    const DRAW_TEXT: bool = false;
    let mut f = std::io::BufWriter::new(std::fs::File::create(out_path).unwrap());

    const PADDING: f64 = 50f64;

    let mut min_x = persistence.points[0][0];
    let mut max_x = persistence.points[0][0];
    let mut min_y = persistence.points[0][1];
    let mut max_y = persistence.points[0][1];
    for &point in &persistence.points {
        min_x = min_x.min(point[0]);
        min_y = min_y.min(point[1]);
        max_x = max_x.max(point[0]);
        max_y = max_y.max(point[1]);
    }
    let xr = max_x - min_x;
    let yr = max_y - min_y;

    write!(
        f,
        r#"<svg xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink"
                height="{}" width="{}">\n"#,
        yr + PADDING * 2.0,
        xr + PADDING * 2.0
    )
    .unwrap();

    let xo = PADDING - min_x;
    let yo = PADDING - min_y;

    for simplex in persistence.simplices.iter() {
        if simplex.dim() == 2 {
            let xy_i = simplex.faces[0];
            let xy = &persistence.simplices[xy_i];
            let yz_i = simplex.faces[1];
            let yz = &persistence.simplices[yz_i];
            let zx_i = simplex.faces[2];
            let zx = &persistence.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 p0 = Point(persistence.points[x_i - 1]);
            let p1 = Point(persistence.points[y_i - 1]);
            let p2 = Point(persistence.points[z_i - 1]);

            if simplex.r_when_born >= R_CUTOFF {
                if xy.r_when_born < 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();
                    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();
                    }
                }

                if zx.r_when_born < 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();
                    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();
                    }
                }

                if yz.r_when_born < 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();
                    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();
                    }
                }
            } else {
                writeln!(f,
                    r#"<polygon points="{},{} {},{} {},{}" style="fill:#e0a0f0;stroke:black;stroke-width:1;" />"#,
                    p0.0[0] + xo,
                    p0.0[1] + yo,
                    p1.0[0] + xo,
                    p1.0[1] + yo,
                    p2.0[0] + xo,
                    p2.0[1] + yo
                ).unwrap();

                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();
                    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();
                    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();
                }
            }
        }
    }

    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 s.len() == 2 {
                std_jumps[i] = s[0];
                exh_jumps[i] = e[0];
            }
        }
    }

    for (i, v) in persistence.points.iter().enumerate() {
        write!(
            f,
            r#"<circle cx="{}" cy="{}" r="2" fill="red"/>\n"#,
            v[0] + PADDING - min_x,
            v[1] + PADDING - min_y
        )
        .unwrap();
        if DRAW_TEXT {
            write!(f,
                r#"<text style="fill:blue;font-size:10px" x="{}" y="{}" >{}[{}/{}]</text>\n"#,
                v[0] + PADDING - min_x,
                v[1] + PADDING - min_y,
                i + 1,
                std_jumps[i+1],
                exh_jumps[i+1],
            )
            .unwrap();
        }
    }
    write!(f, "</svg>\n").unwrap();
}

// Run `gnuplot` with the given script and out file for output.
fn run_gnuplot(script: &str, out_file: &Path) {
    let mut gnuplot = Command::new("gnuplot")
        .stdin(Stdio::piped())
        .stdout(Stdio::piped())
        .spawn()
        .expect("Failed to run `gnuplot`");
    {
        let stdin = gnuplot
            .stdin
            .as_mut()
            .expect("Failed to open stdin for `gnuplot`");
        stdin
            .write_all(script.as_bytes())
            .expect("Failed to write to stdin");
    }
    let output = gnuplot
        .wait_with_output()
        .expect("`wait_for_output` failed!");
    match output.status.code() {
        Some(0) => {}
        other => panic!("`gnuplot` process returned error: {:?}", other),
    }

    let mut f = File::create(out_file).unwrap();
    f.write_all(&output.stdout).unwrap();
}

fn histogram(data: &[usize]) -> Vec<usize> {
    if data.len() == 0 {
        panic!("Trying to create histogram of empty data set");
    }
    let max = data.iter().max().unwrap();
    let mut freq = iter::repeat(0).take(max + 1).collect::<Vec<_>>();
    for &d in data.iter() {
        freq[d] += 1;
    }
    freq
}

fn output_persistence(pairs: &[(f64, f64)], out_file: &Path) {
    let in_file_path = "temp-persistence";
    let mut file = File::create(in_file_path).unwrap();

    for (a, b) in pairs {
        write!(file, "{} {}\n", a, b).unwrap();
    }

    let plot_script = format!(
        r#"
set terminal pdf size 12cm,10cm
set xlabel "Birth"
set ylabel "Death"
set key outside center top horizontal
plot '{}' using 1:2 pt 7 ps 0.25 title "pairs", x"#,
        in_file_path
    );

    run_gnuplot(&plot_script, out_file);
    std::fs::remove_file(in_file_path).unwrap();
}

fn output_histogram(data: &[usize], label: &str, out_file: &Path) {
    if data.len() == 0 {
        panic!("output_histogram: `data` cannot be empty!");
    }

    let mut frequency = histogram(data);
    // Remove single outliers, and don't give out trailing zeroes to gnuplot.
    let mut last_i = 0;
    for (i, f) in frequency.iter_mut().enumerate() {
        if *f == 1 {
            *f = 0;
        }
        if *f > 1 {
            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
set grid
set terminal pdf size 18cm,10cm
set logscale y
set xtics out
set ytics out
set xrange [STATS_min_x - 0.5:]
set xlabel "{}"
set ylabel "Frequency"
set style fill solid 1.0 noborder
set boxwidth 0.5
plot '{file_path}' using ($1):($2) with boxes lc rgb"gray40" notitle
"#,
        label,
        file_path = file_path.display()
    );

    run_gnuplot(&plot_script, out_file);
    std::fs::remove_file(file_path).unwrap();
}

fn output_2histogram(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 = histogram(data_a);
    let mut b_hist = histogram(data_b);
    if a_hist.len() < b_hist.len() {
        for _ in 0..(b_hist.len() - a_hist.len()) {
            a_hist.push(0);
        }
    } else if a_hist.len() > b_hist.len() {
        for _ in 0..(a_hist.len() - b_hist.len()) {
            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())
        .skip_while(|(&a, &b)| a == 0 && b == 0)
    {
        write!(file, "{} {}\n", a, b).unwrap();
    }

    let plot_script = format!(
        r#"
stats '{file_path}' nooutput
set terminal pdf size 18cm,10cm
set logscale y
set xtics out
set ytics out
set xrange [STATS_min_x - 0.5:]
set xlabel "{}"
set ylabel "Frequency"
set boxwidth 0.5
set style fill solid 1.0 noborder
set key outside center top horizontal
a = 0
cum_a(x, i) = (a=a + x * i, a)
b = 0
cum_b(x, i) = (b=b + x * i, b)
set grid
plot '{file_path}' using ($0):(cum_a($1, $0)) with lines   lc rgb"gray20" title "Regular",\
     '{file_path}' using ($0):(cum_b($2, $0)) with lines   lc rgb"0x88bb88" title "Exhaustive",\
     '{file_path}' using ($0 - 0.25):($1+1) with boxes lc rgb"gray20" title "Regular",\
     '{file_path}' using ($0 + 0.25):($2+1) with boxes lc rgb"0x88bb88" title "Exhaustive"
"#,
        label,
        file_path = file_path.display()
    );

    run_gnuplot(&plot_script, out_file);
    std::fs::remove_file(file_path).unwrap();
}

fn tex_table<W: Write>(reg: &Statistics, exh: &Statistics, writer: &mut W) -> Result<(), std::io::Error> {
    fn sum(s: &[usize]) -> f64 {
        s.iter().sum::<usize>() as f64
    }

    fn avg(s: &[usize]) -> f64 {
        sum(s) / s.len() as f64
    }

    fn percent(hundred: f64, ours: f64) -> String {
        if hundred == 0.0 {
            format!("")
        } else {
            format!("{:.2}\\%", 100.0 * ours / hundred)
        }
    }

    let header_align = std::iter::repeat("r l").take(MAX_DIM).collect::<Vec<&str>>().join(" ");
    let header_top = (0..MAX_DIM).map(|d| format!("\\multicolumn{{2}}{{c}}{{$d={}$}} ", d)).collect::<Vec<_>>().join(" & ");
    write!(
        writer,
        r#"\begin{{tabular}}{{r r {}}}
  Name & Total & {} \\
  \toprule
"#, header_align, header_top
    )?;

    let total = reg.num_of_dimens.iter().sum::<usize>();
    let mut line = format!("Simplices & {}", total);
    for &n in reg.num_of_dimens.iter() {
        line.push_str(&format!("& \\hspace{{0.5cm}}{} & {}", thousand_split(n), percent(total as f64, n as f64)));
    }
    writeln!(writer, "{}\\\\", line)?;

    let total = reg.skip_d.iter().sum::<usize>();
    let mut line = format!("Skipped & {}", total);
    for &n in reg.skip_d.iter() {
        line.push_str(&format!("& {} & {}", thousand_split(n), percent(total as f64, n as f64)));
    }
    writeln!(writer, "{}\\\\", line)?;

    let total = reg.zeroed_d.iter().sum::<usize>();
    let mut line = format!("Zeroed & {}", total);
    for &n in reg.zeroed_d.iter() {
        line.push_str(&format!("& {} & {}", thousand_split(n), percent(total as f64, n as f64)));
    }
    writeln!(writer, "{}\\\\", line)?;

    let total = reg.placed_full_d.iter().sum::<usize>();
    let mut line = format!("Stored & {}", total);
    for &n in reg.placed_full_d.iter() {
        line.push_str(&format!("& {} & {}", thousand_split(n), percent(total as f64, n as f64)));
    }
    writeln!(writer, "{}\\\\", line)?;

    writeln!(writer, "\\end{{tabular}}\n")?;
    writeln!(writer, "\\vspace{{1cm}}")?;

    let mut v: Vec<(String, String, String, String)> = Vec::new();

    v.push((
        format!("Column additions"),
        thousand_split(reg.col_adds),
        thousand_split(exh.col_adds),
        percent(reg.col_adds as f64, exh.col_adds as f64),
    ));

    v.push((
        format!("Estimate \\texttt{{column\\_add}} cost"),
        thousand_split(sum(&reg.add_size_sum) as usize),
        thousand_split(sum(&exh.add_size_sum) as usize),
        percent(sum(&reg.add_size_sum), sum(&exh.add_size_sum)),
    ));

    v.push((
        format!("\\texttt{{simplex\\_with\\_low}} queries"),
        thousand_split(reg.simplex_with_low_queries),
        thousand_split(exh.simplex_with_low_queries),
        percent(
            reg.simplex_with_low_queries as f64,
            exh.simplex_with_low_queries as f64,
        ),
    ));

    v.push((
        format!("Time spent (ns)"),
        thousand_split(reg.time as usize),
        thousand_split(exh.time as usize),
        percent(reg.time as f64, exh.time as f64),
    ));

    let as0 = |s: &[(usize, usize)]| s.iter().map(|(a, _)| a).cloned().collect::<Vec<_>>();
    v.push((
        format!("Size of first operand"),
        format!("{:.4}", avg(&as0(&reg.add_size))),
        format!("{:.4}", avg(&as0(&exh.add_size))),
        percent(avg(&as0(&reg.add_size)), avg(&as0(&exh.add_size))),
    ));

    let as1 = |s: &[(usize, usize)]| s.iter().map(|(_, b)| b).cloned().collect::<Vec<_>>();
    v.push((
        format!("Size of second operand"),
        format!("{:.4}", avg(&as1(&reg.add_size))),
        format!("{:.4}", avg(&as1(&exh.add_size))),
        percent(avg(&as1(&reg.add_size)), avg(&as1(&exh.add_size))),
    ));

    v.push((
        format!("Iters for reducing a column"),
        format!("{:.4}", avg(&reg.num_iters)),
        format!("{:.4}", avg(&exh.num_iters)),
        percent(avg(&reg.num_iters), avg(&exh.num_iters)),
    ));

    v.push((
        format!("Cost of one column addition"),
        format!("{:.4}", avg(&reg.add_size_sum)),
        format!("{:.4}", avg(&exh.add_size_sum)),
        percent(avg(&reg.add_size_sum), avg(&exh.add_size_sum)),
    ));

    v.push((
        format!("Iters to exh.\\ reduce a col"),
        format!(""),
        format!("{:.4}", avg(&exh.ex_reductions)),
        format!(""),
    ));

    v.push((
        format!("Searches to find \\texttt{{j=low(k)}}"),
        format!(""),
        format!("{:.4}", avg(&exh.ex_searches)),
        format!(""),
    ));

    // TODO: 4 could be bigger!
    let mut reg_counts = [0; 4];
    for &d in &reg.col_add_dimens {
        reg_counts[(d + 1) as usize] += 1;
    }
    let mut exh_counts = [0; 4];
    for &d in &exh.col_add_dimens {
        exh_counts[(d + 1) as usize] += 1;
    }

    for (d1, (&r, &e)) in reg_counts.iter().zip(exh_counts.iter()).enumerate().skip(1) {
        v.push((
            format!("Column adds with $d={}$", d1 as isize - 1),
            thousand_split(r),
            thousand_split(e),
            percent(r as f64, e as f64)
        ));
    }

    write!(
        writer,
        r#"\begin{{tabular}}{{r r r r}}
  Counter & Regular & Exhausitve & Difference \\
  \toprule
"#
    )?;
    for (l, a, b, p) in v {
        write!(writer, "  {} & {} & {} & {}\\\\\n", l, a, b, p)?;
    }
    writeln!(writer, "\\end{{tabular}}")?;

    Ok(())
}

/// Write out a dotfile with the directed edges given by `simplex_with_low`
fn write_graphviz<W: Write>(p: &Persistence, std: &Reduction, exh: &Reduction, w: &mut W) -> Result<(), std::io::Error> {
    const SCALE: f64 = 5.0;
    writeln!(w, "digraph {{")?;
    writeln!(w, "  rankdir=LR;")?;
    // for simplex in &p.simplices {
    //     if simplex.dim() == 0 {
    //         let pt = p.points[simplex.j - 1];
    //         writeln!(w, "  {} [label=\"v{} [{}]\",pos=\"{},{}\"];", simplex.j,
    //                  simplex.j,
    //                  std.j_stored_at[simplex.j],
    //                  pt[0] * SCALE, -pt[1] * SCALE)?;
    //     } else if simplex.dim() == 1 {
    //         let x = Point(p.points[simplex.faces[0] - 1]);
    //         let y = Point(p.points[simplex.faces[1] - 1]);
    //         let center = (x + y) / 2.0;
    //         writeln!(w, "  {} [label=\"e{}\",pos=\"{},{}\"];",
    //                  simplex.j,
    //                  simplex.j,
    //                  center.0[0] * SCALE, -center.0[1] * SCALE)?;
    //     } else if simplex.dim() == 2 {
    //         let (x_i, y_i, z_i) = get_triangle_xyz(simplex, &p.simplices);
    //         let x = Point(p.points[x_i - 1]);
    //         let y = Point(p.points[y_i - 1]);
    //         let z = Point(p.points[z_i - 1]);
    //         let center = (x + y + z) / 3.0;
    //         writeln!(w, "  {} [label=\"f{}\",pos=\"{},{}\"];",
    //                  simplex.j,
    //                  simplex.j,
    //                  center.0[0] * SCALE, -center.0[1] * SCALE)?;
    //     }
    // }

    // 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();
            assert_eq!(std.j_stored_at[i], exh.j_stored_at[i]);
            if s.len() == 2 {
                writeln!(w, "  {} -> {} [color=\"red\"]", s[1], s[0])?;
                writeln!(w, "  {} -> {} [color=\"green\"]", e[1], e[0])?;
            }
        }
    }

    for &(from, to) in &exh.additions {
        if to > 80 { continue; }
        let d = p.simplices[from].dim();
        if  d == 1 {
            writeln!(w, "  {} -> {} [color=\"green\"]", to, from)?;
        }
        if  d == 2 {
            writeln!(w, "  {} -> {} [color=\"green\"]", to, from)?;
        }
    }

    for &(from, to) in &std.additions {
        if to > 80 { continue; }
        let d = p.simplices[from].dim();
        if  d == 1 {
            writeln!(w, "  {} -> {} [color=\"red\"]", to, from)?;
        }
        if  d == 2 {
            writeln!(w, "  {} -> {} [color=\"red\"]", to, from)?;
        }
    }

    writeln!(w, "}}")?;
    Ok(())
}

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

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


@@ 1455,39 733,39 @@ fn main() {
    if false {
        let mut out = std::io::stderr();
        standard_reduction.write_table(&mut out, &persistence.simplices).unwrap();
        write!(out, "\n\n");
        write!(out, "\n\n").unwrap();
        exhaustive_reduction.write_table(&mut out, &persistence.simplices).unwrap();
    }

    let mut stats_file = File::create(out_path.join("perf-stats")).unwrap();

    write!(stats_file, "# Standard\n").unwrap();
    r_stats.write_stats(&mut stats_file).unwrap();
    output::stats(&r_stats, &mut stats_file).unwrap();
    write!(stats_file, "\n\n").unwrap();
    write!(stats_file, "# Exhaustive\n").unwrap();
    e_stats.write_stats(&mut stats_file).unwrap();
    output::stats(&e_stats, &mut stats_file).unwrap();

    let mut tex_file = File::create(out_path.join("perf-stats.tex")).unwrap();
    tex_table(&r_stats, &e_stats, &mut tex_file);
    output::tex_table(&r_stats, &e_stats, &mut tex_file).unwrap();

    output_histogram(
    output::histogram(
        &e_stats.ex_searches,
        "Iters for finding k=low(i) for any i",
        &out_path.join("ex_searches.pdf"),
    );
    output_histogram(
    output::histogram(
        &e_stats.ex_reductions,
        "Iters for exhaustively reducing a column",
        &out_path.join("ex_reductions.pdf"),
    );

    output_2histogram(
    output::histogram2(
        &r_stats.add_size_sum,
        &e_stats.add_size_sum,
        "Cost estimate for column additions",
        &out_path.join("add_size_sum.pdf"),
    );
    output_2histogram(
    output::histogram2(
        &r_stats.num_iters,
        &e_stats.num_iters,
        "Number of column additions before fully reducing one column",


@@ 1495,7 773,7 @@ fn main() {
    );

    let mut f = File::create("test.dot").unwrap();
    write_graphviz(&persistence, &standard_reduction, &exhaustive_reduction, &mut f).unwrap();
    output::graphviz(&persistence, &standard_reduction, &exhaustive_reduction, &mut f).unwrap();

    let birth_death_pairs = standard_reduction
        .pairings


@@ 1507,9 785,9 @@ fn main() {
            )
        })
        .collect::<Vec<(f64, f64)>>();
    output_persistence(&birth_death_pairs, &out_path.join("persistence.pdf"));
    output::persistence(&birth_death_pairs, &out_path.join("persistence.pdf"));

    output_svg(&persistence, &standard_reduction, &exhaustive_reduction, &out_path.join("rendering.svg"));
    output::svg(&persistence, &standard_reduction, &exhaustive_reduction, &out_path.join("rendering.svg"));
}

#[cfg(test)]

A cra/src/output.rs => cra/src/output.rs +743 -0
@@ 0,0 1,743 @@
use std::io::{Write};
use std::path::Path;
use std::iter;
use std::fs::File;
use std::process::{Command, Stdio};
use std::fmt::Write as FmtWrite;

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


 
/// Format the given `usize` to a `String` with underscores in between each third digit.
///
/// ## Example
///
/// ```
/// assert_eq!(thousand_split(1234567), "1_234_567");
/// ```
fn thousand_split(mut num: usize) -> String {
    let mut out = String::new();
    let mut nums = vec![];
    while num > 0 {
        nums.push(num % 1_000);
        num /= 1_000;
    }
    for i in (1..nums.len()).rev() {
        if i == nums.len() - 1 {
            write!(out, "{},", nums[i]).unwrap();
        } else {
            write!(out, "{:03},", nums[i]).unwrap();
        }
    }
    if nums.len() == 0 {
        write!(out, "0").unwrap();
    } else if nums.len() == 1 {
        write!(out, "{}", nums[0]).unwrap();
    } else {
        write!(out, "{:03}", nums[0]).unwrap();
    }
    out
}

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

    const PADDING: f64 = 50f64;

    let mut min_x = persistence.points[0][0];
    let mut max_x = persistence.points[0][0];
    let mut min_y = persistence.points[0][1];
    let mut max_y = persistence.points[0][1];
    for &point in &persistence.points {
        min_x = min_x.min(point[0]);
        min_y = min_y.min(point[1]);
        max_x = max_x.max(point[0]);
        max_y = max_y.max(point[1]);
    }
    let xr = max_x - min_x;
    let yr = max_y - min_y;

    write!(
        f,
        r#"<svg xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink"
                height="{}" width="{}">\n"#,
        yr + PADDING * 2.0,
        xr + PADDING * 2.0
    )
    .unwrap();

    let xo = PADDING - min_x;
    let yo = PADDING - min_y;

    for simplex in persistence.simplices.iter() {
        if simplex.dim() == 2 {
            let xy_i = simplex.faces[0];
            let xy = &persistence.simplices[xy_i];
            let yz_i = simplex.faces[1];
            let yz = &persistence.simplices[yz_i];
            let zx_i = simplex.faces[2];
            let zx = &persistence.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 p0 = Point(persistence.points[x_i - 1]);
            let p1 = Point(persistence.points[y_i - 1]);
            let p2 = Point(persistence.points[z_i - 1]);

            if simplex.r_when_born >= R_CUTOFF {
                if xy.r_when_born < 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();
                    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();
                    }
                }

                if zx.r_when_born < 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();
                    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();
                    }
                }

                if yz.r_when_born < 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();
                    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();
                    }
                }
            } else {
                writeln!(f,
                    r#"<polygon points="{},{} {},{} {},{}" style="fill:#e0a0f0;stroke:black;stroke-width:1;" />"#,
                    p0.0[0] + xo,
                    p0.0[1] + yo,
                    p1.0[0] + xo,
                    p1.0[1] + yo,
                    p2.0[0] + xo,
                    p2.0[1] + yo
                ).unwrap();

                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();
                    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();
                    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();
                }
            }
        }
    }

    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 s.len() == 2 {
                std_jumps[i] = s[0];
                exh_jumps[i] = e[0];
            }
        }
    }

    for (i, v) in persistence.points.iter().enumerate() {
        write!(
            f,
            r#"<circle cx="{}" cy="{}" r="2" fill="red"/>\n"#,
            v[0] + PADDING - min_x,
            v[1] + PADDING - min_y
        )
        .unwrap();
        if DRAW_TEXT {
            write!(f,
                r#"<text style="fill:blue;font-size:10px" x="{}" y="{}" >{}[{}/{}]</text>\n"#,
                v[0] + PADDING - min_x,
                v[1] + PADDING - min_y,
                i + 1,
                std_jumps[i+1],
                exh_jumps[i+1],
            )
            .unwrap();
        }
    }
    write!(f, "</svg>\n").unwrap();
}

fn make_histogram(data: &[usize]) -> Vec<usize> {
    if data.len() == 0 {
        panic!("Trying to create histogram of empty data set");
    }
    let max = data.iter().max().unwrap();
    let mut freq = iter::repeat(0).take(max + 1).collect::<Vec<_>>();
    for &d in data.iter() {
        freq[d] += 1;
    }
    freq
}

pub fn persistence(pairs: &[(f64, f64)], out_file: &Path) {
    let in_file_path = "temp-persistence";
    let mut file = File::create(in_file_path).unwrap();

    for (a, b) in pairs {
        write!(file, "{} {}\n", a, b).unwrap();
    }

    let plot_script = format!(
        r#"
set terminal pdf size 12cm,10cm
set xlabel "Birth"
set ylabel "Death"
set key outside center top horizontal
plot '{}' using 1:2 pt 7 ps 0.25 title "pairs", x"#,
        in_file_path
    );

    run_gnuplot(&plot_script, out_file);
    std::fs::remove_file(in_file_path).unwrap();
}

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;
    for (i, f) in frequency.iter_mut().enumerate() {
        if *f == 1 {
            *f = 0;
        }
        if *f > 1 {
            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
set grid
set terminal pdf size 18cm,10cm
set logscale y
set xtics out
set ytics out
set xrange [STATS_min_x - 0.5:]
set xlabel "{}"
set ylabel "Frequency"
set style fill solid 1.0 noborder
set boxwidth 0.5
plot '{file_path}' using ($1):($2) with boxes lc rgb"gray40" notitle
"#,
        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() {
        for _ in 0..(b_hist.len() - a_hist.len()) {
            a_hist.push(0);
        }
    } else if a_hist.len() > b_hist.len() {
        for _ in 0..(a_hist.len() - b_hist.len()) {
            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())
        .skip_while(|(&a, &b)| a == 0 && b == 0)
    {
        write!(file, "{} {}\n", a, b).unwrap();
    }

    let plot_script = format!(
        r#"
stats '{file_path}' nooutput
set terminal pdf size 18cm,10cm
set logscale y
set xtics out
set ytics out
set xrange [STATS_min_x - 0.5:]
set xlabel "{}"
set ylabel "Frequency"
set boxwidth 0.5
set style fill solid 1.0 noborder
set key outside center top horizontal
a = 0
cum_a(x, i) = (a=a + x * i, a)
b = 0
cum_b(x, i) = (b=b + x * i, b)
set grid
plot '{file_path}' using ($0):(cum_a($1, $0)) with lines   lc rgb"gray20" title "Regular",\
     '{file_path}' using ($0):(cum_b($2, $0)) with lines   lc rgb"0x88bb88" title "Exhaustive",\
     '{file_path}' using ($0 - 0.25):($1+1) with boxes lc rgb"gray20" title "Regular",\
     '{file_path}' using ($0 + 0.25):($2+1) with boxes lc rgb"0x88bb88" title "Exhaustive"
"#,
        label,
        file_path = file_path.display()
    );

    run_gnuplot(&plot_script, out_file);
    std::fs::remove_file(file_path).unwrap();
}

pub fn tex_table<W: Write>(reg: &Statistics, exh: &Statistics, writer: &mut W) -> Result<(), std::io::Error> {
    fn sum(s: &[usize]) -> f64 {
        s.iter().sum::<usize>() as f64
    }

    fn avg(s: &[usize]) -> f64 {
        sum(s) / s.len() as f64
    }

    fn percent(hundred: f64, ours: f64) -> String {
        if hundred == 0.0 {
            format!("")
        } else {
            format!("{:.2}\\%", 100.0 * ours / hundred)
        }
    }

    let header_align = std::iter::repeat("r l").take(MAX_DIM).collect::<Vec<&str>>().join(" ");
    let header_top = (0..MAX_DIM).map(|d| format!("\\multicolumn{{2}}{{c}}{{$d={}$}} ", d)).collect::<Vec<_>>().join(" & ");
    write!(
        writer,
        r#"\begin{{tabular}}{{r r {}}}
  Name & Total & {} \\
  \toprule
"#, header_align, header_top
    )?;

    let total = reg.num_of_dimens.iter().sum::<usize>();
    let mut line = format!("Simplices & {}", total);
    for &n in reg.num_of_dimens.iter() {
        line.push_str(&format!("& \\hspace{{0.5cm}}{} & {}", thousand_split(n), percent(total as f64, n as f64)));
    }
    writeln!(writer, "{}\\\\", line)?;

    let total = reg.skip_d.iter().sum::<usize>();
    let mut line = format!("Skipped & {}", total);
    for &n in reg.skip_d.iter() {
        line.push_str(&format!("& {} & {}", thousand_split(n), percent(total as f64, n as f64)));
    }
    writeln!(writer, "{}\\\\", line)?;

    let total = reg.zeroed_d.iter().sum::<usize>();
    let mut line = format!("Zeroed & {}", total);
    for &n in reg.zeroed_d.iter() {
        line.push_str(&format!("& {} & {}", thousand_split(n), percent(total as f64, n as f64)));
    }
    writeln!(writer, "{}\\\\", line)?;

    let total = reg.placed_full_d.iter().sum::<usize>();
    let mut line = format!("Stored & {}", total);
    for &n in reg.placed_full_d.iter() {
        line.push_str(&format!("& {} & {}", thousand_split(n), percent(total as f64, n as f64)));
    }
    writeln!(writer, "{}\\\\", line)?;

    writeln!(writer, "\\end{{tabular}}\n")?;
    writeln!(writer, "\\vspace{{1cm}}")?;

    let mut v: Vec<(String, String, String, String)> = Vec::new();

    v.push((
        format!("Column additions"),
        thousand_split(reg.col_adds),
        thousand_split(exh.col_adds),
        percent(reg.col_adds as f64, exh.col_adds as f64),
    ));

    v.push((
        format!("Estimate \\texttt{{column\\_add}} cost"),
        thousand_split(sum(&reg.add_size_sum) as usize),
        thousand_split(sum(&exh.add_size_sum) as usize),
        percent(sum(&reg.add_size_sum), sum(&exh.add_size_sum)),
    ));

    v.push((
        format!("\\texttt{{simplex\\_with\\_low}} queries"),
        thousand_split(reg.simplex_with_low_queries),
        thousand_split(exh.simplex_with_low_queries),
        percent(
            reg.simplex_with_low_queries as f64,
            exh.simplex_with_low_queries as f64,
        ),
    ));

    v.push((
        format!("Time spent (ns)"),
        thousand_split(reg.time as usize),
        thousand_split(exh.time as usize),
        percent(reg.time as f64, exh.time as f64),
    ));

    let as0 = |s: &[(usize, usize)]| s.iter().map(|(a, _)| a).cloned().collect::<Vec<_>>();
    v.push((
        format!("Size of first operand"),
        format!("{:.4}", avg(&as0(&reg.add_size))),
        format!("{:.4}", avg(&as0(&exh.add_size))),
        percent(avg(&as0(&reg.add_size)), avg(&as0(&exh.add_size))),
    ));

    let as1 = |s: &[(usize, usize)]| s.iter().map(|(_, b)| b).cloned().collect::<Vec<_>>();
    v.push((
        format!("Size of second operand"),
        format!("{:.4}", avg(&as1(&reg.add_size))),
        format!("{:.4}", avg(&as1(&exh.add_size))),
        percent(avg(&as1(&reg.add_size)), avg(&as1(&exh.add_size))),
    ));

    v.push((
        format!("Iters for reducing a column"),
        format!("{:.4}", avg(&reg.num_iters)),
        format!("{:.4}", avg(&exh.num_iters)),
        percent(avg(&reg.num_iters), avg(&exh.num_iters)),
    ));

    v.push((
        format!("Cost of one column addition"),
        format!("{:.4}", avg(&reg.add_size_sum)),
        format!("{:.4}", avg(&exh.add_size_sum)),
        percent(avg(&reg.add_size_sum), avg(&exh.add_size_sum)),
    ));

    v.push((
        format!("Iters to exh.\\ reduce a col"),
        format!(""),
        format!("{:.4}", avg(&exh.ex_reductions)),
        format!(""),
    ));

    v.push((
        format!("Searches to find \\texttt{{j=low(k)}}"),
        format!(""),
        format!("{:.4}", avg(&exh.ex_searches)),
        format!(""),
    ));

    // TODO: 4 could be bigger!
    let mut reg_counts = [0; 4];
    for &d in &reg.col_add_dimens {
        reg_counts[(d + 1) as usize] += 1;
    }
    let mut exh_counts = [0; 4];
    for &d in &exh.col_add_dimens {
        exh_counts[(d + 1) as usize] += 1;
    }

    for (d1, (&r, &e)) in reg_counts.iter().zip(exh_counts.iter()).enumerate().skip(1) {
        v.push((
            format!("Column adds with $d={}$", d1 as isize - 1),
            thousand_split(r),
            thousand_split(e),
            percent(r as f64, e as f64)
        ));
    }

    write!(
        writer,
        r#"\begin{{tabular}}{{r r r r}}
  Counter & Regular & Exhausitve & Difference \\
  \toprule
"#
    )?;
    for (l, a, b, p) in v {
        write!(writer, "  {} & {} & {} & {}\\\\\n", l, a, b, p)?;
    }
    writeln!(writer, "\\end{{tabular}}")?;

    Ok(())
}

/// Write out a dotfile with the directed edges given by `simplex_with_low`
pub fn graphviz<W: Write>(p: &Persistence, std: &Reduction, 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;")?;
    if false {
        for simplex in &p.simplices {
            if simplex.dim() == 0 {
                let pt = p.points[simplex.j - 1];
                writeln!(w, "  {} [label=\"v{} [{}]\",pos=\"{},{}\"];", simplex.j,
                         simplex.j,
                         std.j_stored_at[simplex.j],
                         pt[0] * SCALE, -pt[1] * SCALE)?;
            } else if simplex.dim() == 1 {
                let x = Point(p.points[simplex.faces[0] - 1]);
                let y = Point(p.points[simplex.faces[1] - 1]);
                let center = (x + y) / 2.0;
                writeln!(w, "  {} [label=\"e{}\",pos=\"{},{}\"];",
                         simplex.j,
                         simplex.j,
                         center.0[0] * SCALE, -center.0[1] * SCALE)?;
            } else if simplex.dim() == 2 {
                let (x_i, y_i, z_i) = get_triangle_xyz(simplex, &p.simplices);
                let x = Point(p.points[x_i - 1]);
                let y = Point(p.points[y_i - 1]);
                let z = Point(p.points[z_i - 1]);
                let center = (x + y + z) / 3.0;
                writeln!(w, "  {} [label=\"f{}\",pos=\"{},{}\"];",
                         simplex.j,
                         simplex.j,
                         center.0[0] * SCALE, -center.0[1] * SCALE)?;
            }
        }
    }

    // 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();
            assert_eq!(std.j_stored_at[i], exh.j_stored_at[i]);
            if s.len() == 2 {
                writeln!(w, "  {} -> {} [color=\"red\"]", s[1], s[0])?;
                writeln!(w, "  {} -> {} [color=\"green\"]", e[1], e[0])?;
            }
        }
    }

    for &(from, to) in &exh.additions {
        if to > 80 { continue; }
        let d = p.simplices[from].dim();
        if  d == 1 {
            writeln!(w, "  {} -> {} [color=\"green\"]", to, from)?;
        }
        if  d == 2 {
            writeln!(w, "  {} -> {} [color=\"green\"]", to, from)?;
        }
    }

    for &(from, to) in &std.additions {
        if to > 80 { continue; }
        let d = p.simplices[from].dim();
        if  d == 1 {
            writeln!(w, "  {} -> {} [color=\"red\"]", to, from)?;
        }
        if  d == 2 {
            writeln!(w, "  {} -> {} [color=\"red\"]", to, from)?;
        }
    }

    writeln!(w, "}}")?;
    Ok(())
}

// Run `gnuplot` with the given script and out file for output.
pub fn run_gnuplot(script: &str, out_file: &Path) {
    let mut gnuplot = Command::new("gnuplot")
        .stdin(Stdio::piped())
        .stdout(Stdio::piped())
        .spawn()
        .expect("Failed to run `gnuplot`");
    {
        let stdin = gnuplot
            .stdin
            .as_mut()
            .expect("Failed to open stdin for `gnuplot`");
        stdin
            .write_all(script.as_bytes())
            .expect("Failed to write to stdin");
    }
    let output = gnuplot
        .wait_with_output()
        .expect("`wait_for_output` failed!");
    match output.status.code() {
        Some(0) => {}
        other => panic!("`gnuplot` process returned error: {:?}", other),
    }

    let mut f = File::create(out_file).unwrap();
    f.write_all(&output.stdout).unwrap();
}

pub fn stats<W: Write>(st: &Statistics, writer: &mut W) -> Result<(), std::io::Error> {
        let mut print_pairs = Vec::new();

        print_pairs.push(("Column additions", thousand_split(st.col_adds)));

        let adds = st.add_size_sum.iter().sum::<usize>() as usize;
        print_pairs.push(("Estimate total 'column_add' cost", thousand_split(adds)));

        print_pairs.push((
            "`simplex_with_low` queries",
            thousand_split(st.simplex_with_low_queries),
        ));

        print_pairs.push(("Time spent (ns)", thousand_split(st.time as usize)));

        // ~10 cycles per L2 hit seems reasonable?
        let instr_cost_estimate = st.simplex_with_low_queries * 10 + adds;
        print_pairs.push((
            "Instruction cost estimate",
            thousand_split(instr_cost_estimate),
        ));

        let add_size0 =
            st.add_size.iter().map(|(a, _)| a).sum::<usize>() as f64 / st.add_size.len() as f64;
        print_pairs.push(("Average size of first operand", format!("{}", add_size0)));

        let add_size1 =
            st.add_size.iter().map(|(_, b)| b).sum::<usize>() as f64 / st.add_size.len() as f64;
        print_pairs.push(("Average size of second operand", format!("{}", add_size1)));

        let num_iters = st.num_iters.iter().sum::<usize>() as f64 / st.num_iters.len() as f64;
        print_pairs.push((
            "Average #iterations before reducing a column",
            format!("{}", num_iters),
        ));

        let column_add_cost =
            st.add_size_sum.iter().sum::<usize>() as f64 / st.add_size_sum.len() as f64;
        print_pairs.push((
            "Average cost of one column addition",
            format!("{}", column_add_cost),
        ));

        if st.ex_reductions.len() > 0 {
            let ex_reductions =
                st.ex_reductions.iter().sum::<usize>() as f64 / st.ex_reductions.len() as f64;
            print_pairs.push((
                "Average iterations to exh. reduce a col",
                format!("{}", ex_reductions),
            ));
        }

        if st.ex_searches.len() > 0 {
            let ex_searches =
                st.ex_searches.iter().sum::<usize>() as f64 / st.ex_searches.len() as f64;
            print_pairs.push((
                "Average searches to find `j=low(k)`",
                format!("{}", ex_searches),
            ));
        }

        print_pairs.push(("Number of zeroed columns", thousand_split(st.zeroed)));
        print_pairs.push((
            "Number of single columns",
            thousand_split(st.placed_single),
        ));
        print_pairs.push((
            "Number of remaining columns",
            thousand_split(st.placed_full),
        ));

        let mut last_pos = 0;
        for (i, &n) in st.num_of_dimens.iter().enumerate() {
            if n > 0 {
                last_pos = i;
            }
        }
        let strs = (0..last_pos + 1)
            .map(|d| format!("Simplices with d={}", d))
            .collect::<Vec<_>>();
        for d in 0..last_pos + 1 {
            print_pairs.push((&strs[d], thousand_split(st.num_of_dimens[d])));
        }

        let max_str = print_pairs.iter().map(|(s, _)| s.len()).max().unwrap();
        for (left, right) in print_pairs {
            write!(writer, "{:>2$}: {}\n", left, right, max_str)?;
        }
        Ok(())
}