~mht/cra

fc08a9880627545f1f130ca802b2701e64261179 — Martin Hafskjold Thoresen 1 year, 3 months ago 1febdad
Rework alpha_2d a little, and fix warnings
4 files changed, 186 insertions(+), 174 deletions(-)

M cra/Makefile
M cra/src/main.rs
M cra/src/output.rs
M cra/src/persistence.rs
M cra/Makefile => cra/Makefile +4 -2
@@ 15,11 15,13 @@ target/release/cra: src/main.rs src/output.rs src/persistence.rs
# Output from the reduction algorithm
output/.%: data/%.tri target/release/cra
	-mkdir -p output/$*
	env RUST_BACKTRACE=1 target/debug/cra -t="-" < $< \
	env RUST_BACKTRACE=1 target/release/cra -t="-" < $< \
		--graphviz=output/$*/graph.dot \
		--diagram=output/$*/persistence.pdf \
		--svg=output/$*/rendering.svg \
		--tex-table=output/$*/perf-stats.tex
		--tex-table=output/$*/perf-stats.tex \
		--stats=output/$*/stats \
		--points="-"
	convert -size 800x800 output/$*/rendering.svg output/$*/rendering.png
	touch $@


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


@@ 17,16 16,6 @@ use persistence::{compute_alpha_values_2d, reduce, Persistence, Simplex};
const R_CUTOFF: f64 = 250.0;
const MAX_DIM: usize = 3;

fn f64_eq(a: f64, b: f64) -> bool {
    // TOOD: this is stupid and should be removed.
    let u = 1_000_000.0;
    let a = (a * u).round() as usize;
    let b = (b * u).round() as usize;
    a == b
    // limit is ~0.00119
    // (a - b).abs() < 10_000.0 * std::f64::EPSILON
}

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


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

/// Read a 2D triangulation from `stdin` and compute the alpha values for the simplices.
pub fn read_input_stdin2() -> Result<Persistence, Box<dyn std::error::Error>> {
    let mut stdin = std::io::stdin();
    let (points, mut simplices) = read_triangulation(&mut stdin)?;
    let (simplices, alphas) = compute_alpha_values_2d(&points, &simplices);
    Ok(Persistence {
        simplices,
        points,
        alphas,
    })
}

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


@@ 317,12 294,21 @@ fn main() {
        (author: crate_authors!())
        (@arg triangulation: -t --triangulation[FILE] "Input a 2D triangulation in text form.")
        (@arg exhaustive: -e --exhaustive "Run the exhaustive variant of the algorithm.")
        (@arg svg: --svg[FILE] "Draw the compex to an .svg.")
        (@arg alpha: -a --alpha "Ignore the ordering of the simplices given and calculate and use the alpha ordering.")
        (@arg svg: --svg[FILE] requires[points] "Draw the compex to an .svg.")
        (@arg alpha: -a --alpha requires[points] "Ignore the ordering of the simplices given and calculate and use the alpha ordering.")
        (@arg stats: -s --stats "Print out statistics.")
        (@arg graphviz: --graphviz[FILE] "Use `graphviz` to graph out which simplices are added together.")
        (@arg diagram: --diagram[FILE] "Use `gnuplot` to print out the persistence diagram.")
        (@arg textable: --("tex-table")[FILE] "Output a `tex` formatted table of statistics.")
        (@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.'",
        ).long_help("\
This file should contain the boundaries of your simplices.
The boundary of the empty complex and the vertices are implied, and so
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") {


@@ 332,7 318,7 @@ fn main() {
        } else {
            Box::new(File::open(path).expect("No such file"))
        };
        let (points, mut simplices) =
        let (points, simplices) =
            read_triangulation(&mut reader).expect("failed to read triangulation");
        let (simplices, alphas) = compute_alpha_values_2d(&points, &simplices);
        Persistence {

M cra/src/output.rs => cra/src/output.rs +118 -119
@@ 1,7 1,6 @@
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};



@@ 60,11 59,11 @@ pub fn svg(persistence: &Persistence, std: &Reduction, exh: &Reduction, out_path
    let max_x = persistence
        .points
        .iter()
        .fold(std::f64::MAX, |m, a| m.max(a[0]));
        .fold(std::f64::MIN, |m, a| m.max(a[0]));
    let max_y = persistence
        .points
        .iter()
        .fold(std::f64::MAX, |m, a| m.max(a[1]));
        .fold(std::f64::MIN, |m, a| m.max(a[1]));

    let xr = max_x - min_x;
    let yr = max_y - min_y;


@@ 265,18 264,6 @@ pub fn svg(persistence: &Persistence, std: &Reduction, exh: &Reduction, out_path
    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
}

/// Use `gnuplot` to make a persistence diagram with the given pairs.
/// This is really just a point plot with labeled axes.
pub fn persistence(pairs: &[(f64, f64)], out_file: &Path) {


@@ 301,110 288,122 @@ plot '{}' using 1:2 pt 7 ps 0.25 title "pairs", x"#,
    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();
}
// 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 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();
// }

/// Output a `tex` formatted table of statistics to the `writer`.
pub fn tex_table<W: Write>(

M cra/src/persistence.rs => cra/src/persistence.rs +52 -27
@@ 1,4 1,4 @@
use crate::{Point, Statistics, R_CUTOFF};
use crate::{Point, Statistics};
use std::io::Write;

/// The three variants of the algorithm.


@@ 369,8 369,7 @@ pub fn compute_alpha_values_2d(
        simplex: &Simplex,
        points: &[[f64; 2]],
        simplices: &[Simplex],
        alphas: &mut [f64],
    ) -> f64 {
    ) -> (f64, Option<usize>) {
        let xy_i = simplex.faces[0];
        let xy = &simplices[xy_i];
        let yz_i = simplex.faces[1];


@@ 411,54 410,80 @@ pub fn compute_alpha_values_2d(
        let circumradius = (xy * yz * xz)
            / ((xy + yz + xz) * (xy + yz - xz) * (xy - yz + xz) * (-xy + yz + xz)).sqrt();

        if am * 180.0 / std::f64::consts::PI >= 90.0 {
		let mut edit = None;
		// TODO: check this out again
        if am * 180.0 / std::f64::consts::PI >= 120.0 {
            // When one angle is > 90.0, the circumcenter of the triangle is outside. Then
            // the balls from the edges that's closest to the circumcenter doesn't meet
            // before they are at the circumcenter. Thus we must adjust the `r_when_born` on
            // this edge.
            if am == ax {
                // debug_assert_eq!(yz / 2.0, simpl(yz_i).r_when_born);
                alphas[yz_i] = circumradius;
                edit = Some(yz_i);
            } else if am == ay {
                // debug_assert_eq!(xz / 2.0, simpl(zx_i).r_when_born);
                alphas[zx_i] = circumradius;
                edit = Some(zx_i);
            } else {
                // debug_assert_eq!(xy / 2.0, simpl(xy_i).r_when_born);
                alphas[xy_i] = circumradius;
                edit = Some(xy_i);
            }
        }
        circumradius
        (circumradius, edit)
    }

    unsafe fn simplex_mut(i: usize, simplices: &[Simplex]) -> &mut Simplex {
        let simplices_ptr = simplices.as_ptr();
        let target: *mut Simplex = simplices_ptr.offset(i as isize) as *mut _;
        &mut *target
	#[derive(PartialEq)]
    struct CmpSimplex {
        dim: usize,
        alpha: f64,
        index: usize,
    }

    let mut tuples = Vec::new();
    let mut alphas = vec![-1.0; simplices.len()];
    impl Eq for CmpSimplex { }

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

    impl Ord for CmpSimplex {
        fn cmp(&self, other: &Self) -> std::cmp::Ordering {
            self.dim.cmp(&other.dim)
            	.then(self.alpha.partial_cmp(&other.alpha).expect("alpha value was NaN or something!"))
        }
    }

	let mut hacked_alphas = std::collections::HashMap::new();
    let mut v = Vec::with_capacity(simplices.len());
    for current_simplex_i in 0..simplices.len() {
        let simplex = &simplices[current_simplex_i];
        let alpha = match simplex.dim() {
            -1 | 0 => 0.0,
            1 => alpha_1d(&simplex, &points),
            2 => alpha_2d(&simplex, &points, simplices, &mut alphas),
            2 => {
                let (a, edit) = alpha_2d(&simplex, &points, simplices);
                if let Some(i) = edit {
                    hacked_alphas.insert(i, a);
                }
                a
            },
            more => panic!("Only support 2d simplices! ({} > 2)", more),
        };
        tuples.push((alpha, current_simplex_i));
        v.push(CmpSimplex {
            dim: simplex.faces.len(),
            alpha,
            index: current_simplex_i
        });
    }

    tuples.sort_by(|(f1, i1), (f2, i2)| f1.partial_cmp(f2).unwrap_or(i1.cmp(i2)));

    let mut ordered_simplices = Vec::new();
    for &(alpha, i) in tuples.iter() {
        ordered_simplices.push(simplices[i].clone());
        if alphas[i] == -1.0 {
            alphas[i] = alpha;
    v.sort();
    let mut ordered: Vec<Simplex> = Vec::with_capacity(simplices.len());
    let mut alphas: Vec<f64> = Vec::with_capacity(simplices.len());
    for cs in v.into_iter() {
        ordered.push(simplices[cs.index].clone());
        if let Some(&a) = hacked_alphas.get(&cs.index) {
            alphas.push(a);
        } else {
            alphas.push(cs.alpha);
        }
    }

    (ordered_simplices, alphas)
    (ordered, alphas)
}