~mht/cra

1910306bc1ae84cacbd8287130eaacd9531e8d1f — Martin Hafskjold Thoresen 1 year, 3 months ago 92ecc9c
Add a bunch of comments
2 files changed, 281 insertions(+), 177 deletions(-)

M cra/src/main.rs
M cra/src/output.rs
M cra/src/main.rs => cra/src/main.rs +78 -55
@@ 383,25 383,43 @@ pub fn read_input_stdin2() -> Result<Persistence, Box<dyn std::error::Error>> {
/// 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> {
    #[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)?;
            writeln!(
                writer,
                "j={}, faces={:?} => {:?}",
                simplex.j, simplex.faces, il
            )?;
        }
        Ok(())
    }

    fn write_table<W: Write>(&self, w: &mut W, simplices: &[Simplex]) -> Result<(), std::io::Error> {
    /// 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 = {


@@ 441,49 459,65 @@ impl Reduction {
    }
}

/// 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 Variant::*;

/// Given a `Vec<Vec<usize>>` in which the first `Vec` is a list of simplices and the second `Vec`
/// is the indices of its faces, run the "matrix reduction" algorithm.
///
/// We assume the slice is sorted, such that all faces precede their simplices, and that
/// the indices are sorted.
/// 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();

    for simplex in p.simplices.iter() {
	// 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; 
            continue;
        }
        let s_dim = simplex.dim().max(0) as usize; 
        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 = *current_simplex.last().unwrap();
            let low: usize = *current_simplex.last().unwrap();
            let cand = &simplex_with_low[low];
            stats.swl();

            if cand.is_null() {
			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));



@@ 510,7 544,8 @@ 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 = simplex_with_low.as_ptr().offset(low as isize);
                                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()
                            };


@@ 526,14 561,6 @@ pub fn reduce(p: &Persistence, variant: Variant, stats: &mut Statistics) -> Redu
                        break;
                    }
                }

                break;
            }

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


@@ 572,6 599,8 @@ 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)]
enum IndList {
    Null,


@@ 592,7 621,7 @@ impl IndList {
        }
    }

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


@@ 617,7 646,12 @@ impl IndList {

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


@@ 732,9 766,13 @@ fn main() {

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

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


@@ 773,7 811,13 @@ fn main() {
    );

    let mut f = File::create("test.dot").unwrap();
    output::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


@@ 787,31 831,10 @@ fn main() {
        .collect::<Vec<(f64, f64)>>();
    output::persistence(&birth_death_pairs, &out_path.join("persistence.pdf"));

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

#[cfg(test)]
mod test {
    use super::*;
    #[test]
    fn indexlist_xor() {
        let s = &mut Statistics::new();

        let mut il = IndexList(vec![0, 1, 2, 3]);

        il.xor(&IndexList(vec![2, 3, 4, 5]), s);
        assert_eq!(il.0, vec![0, 1, 4, 5]);

        il.xor(&IndexList(vec![4, 5, 6]), s);
        assert_eq!(il.0, vec![0, 1, 6]);

        il.xor(&IndexList(vec![2]), s);
        assert_eq!(il.0, vec![0, 1, 2, 6]);

        il.xor(&IndexList(vec![0, 2, 6]), s);
        assert_eq!(il.0, vec![1]);

        il.xor(&IndexList(vec![1]), s);
        assert_eq!(il.0, vec![]);
    }
    output::svg(
        &persistence,
        &standard_reduction,
        &exhaustive_reduction,
        &out_path.join("rendering.svg"),
    );
}

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

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

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

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


@@ 41,9 39,11 @@ 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;
    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.
    const PADDING: f64 = 50f64;

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


@@ 113,8 113,14 @@ pub fn svg(persistence: &Persistence, std: &Reduction, exh: &Reduction, out_path
                    .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();
                        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();
                    }
                }



@@ 131,8 137,14 @@ pub fn svg(persistence: &Persistence, std: &Reduction, exh: &Reduction, out_path
                    .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();
                        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();
                    }
                }



@@ 149,8 161,14 @@ pub fn svg(persistence: &Persistence, std: &Reduction, exh: &Reduction, out_path
                    .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();
                        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 {


@@ 166,14 184,32 @@ pub fn svg(persistence: &Persistence, std: &Reduction, exh: &Reduction, out_path

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


@@ 202,13 238,14 @@ pub fn svg(persistence: &Persistence, std: &Reduction, exh: &Reduction, out_path
        )
        .unwrap();
        if DRAW_TEXT {
            write!(f,
            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],
                std_jumps[i + 1],
                exh_jumps[i + 1],
            )
            .unwrap();
        }


@@ 355,7 392,11 @@ plot '{file_path}' using ($0):(cum_a($1, $0)) with lines   lc rgb"gray20" title 
    std::fs::remove_file(file_path).unwrap();
}

pub fn tex_table<W: Write>(reg: &Statistics, exh: &Statistics, writer: &mut W) -> Result<(), std::io::Error> {
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
    }


@@ 372,41 413,64 @@ pub fn tex_table<W: Write>(reg: &Statistics, exh: &Statistics, writer: &mut W) -
        }
    }

    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(" & ");
    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
"#,
        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)));
        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)));
        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)));
        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)));
        line.push_str(&format!(
            "& {} & {}",
            thousand_split(n),
            percent(total as f64, n as f64)
        ));
    }
    writeln!(writer, "{}\\\\", line)?;



@@ 505,7 569,7 @@ pub fn tex_table<W: Write>(reg: &Statistics, exh: &Statistics, writer: &mut W) -
            format!("Column adds with $d={}$", d1 as isize - 1),
            thousand_split(r),
            thousand_split(e),
            percent(r as f64, e as f64)
            percent(r as f64, e as f64),
        ));
    }



@@ 525,7 589,12 @@ pub fn tex_table<W: Write>(reg: &Statistics, exh: &Statistics, writer: &mut W) -
}

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


@@ 539,7 608,6 @@ pub fn graphviz<W: Write>(p: &Persistence, std: &Reduction, exh: &Reduction, w: 
            // 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] {


@@ 558,28 626,41 @@ pub fn graphviz<W: Write>(p: &Persistence, std: &Reduction, exh: &Reduction, w: 
        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)?;
                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)?;
                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)?;
                writeln!(
                    w,
                    "  {} [label=\"f{}\",pos=\"{},{}\"];",
                    simplex.j,
                    simplex.j,
                    center.0[0] * SCALE,
                    -center.0[1] * SCALE
                )?;
            }
        }
    }


@@ 598,23 679,27 @@ pub fn graphviz<W: Write>(p: &Persistence, std: &Reduction, exh: &Reduction, w: 
    }

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

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


@@ 652,92 737,88 @@ pub fn run_gnuplot(script: &str, out_file: &Path) {
}

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)));
    let mut print_pairs = Vec::new();

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

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

        // ~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),
        ));
    print_pairs.push((
        "`simplex_with_low` queries",
        thousand_split(st.simplex_with_low_queries),
    ));

        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)));
    print_pairs.push(("Time spent (ns)", thousand_split(st.time as usize)));

        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)));
    // ~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 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 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 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),
        ));
    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)));

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

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

        print_pairs.push(("Number of zeroed columns", thousand_split(st.zeroed)));
    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((
            "Number of single columns",
            thousand_split(st.placed_single),
            "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((
            "Number of remaining columns",
            thousand_split(st.placed_full),
            "Average searches to find `j=low(k)`",
            format!("{}", ex_searches),
        ));
    }

        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])));
        }
    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 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)?;
    let mut last_pos = 0;
    for (i, &n) in st.num_of_dimens.iter().enumerate() {
        if n > 0 {
            last_pos = i;
        }
        Ok(())
    }
    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(())
}