extern crate clap;
extern crate lazy_static;
extern crate serde;
extern crate serde_json;
// std
use std::fs::File;
use std::io::{BufRead, BufReader, Read, Write};
use std::path::Path;
// mods
mod log;
mod output;
mod persistence;
// mods use
use persistence::Variant::*;
use persistence::{compute_alpha_values_2d, reduce, Persistence, Simplex};
// consts
const R_CUTOFF: f64 = 250.0;
#[derive(Copy, Clone, Debug)]
pub struct Point(pub [f64; 2]);
impl std::ops::Add<Point> for Point {
type Output = Point;
fn add(self, other: Point) -> Point {
Point([self.0[0] + other.0[0], self.0[1] + other.0[1]])
}
}
impl std::ops::Sub<Point> for Point {
type Output = Point;
fn sub(self, other: Point) -> Point {
Point([self.0[0] - other.0[0], self.0[1] - other.0[1]])
}
}
impl std::ops::Div<f64> for Point {
type Output = Point;
fn div(self, f: f64) -> Point {
Point([self.0[0] / f, self.0[1] / f])
}
}
impl Point {
fn len(self) -> f64 {
(self.0[0].powi(2) + self.0[1].powi(2)).sqrt()
}
fn dot(self, other: Self) -> f64 {
self.0[0] * other.0[0] + self.0[1] * other.0[1]
}
fn angle(self, other: Point) -> f64 {
let n = self.dot(other) / (self.len() * other.len());
n.acos()
}
}
/// Read a boundary file. The first line should be the number of vertices,
/// and consequent lines should be whitespace separated indices for the boundary
/// of that simplex, using 1-indexing.
/// Everything after a '#' in a line is ignored, and empty lines are ignored.
///
/// For now, we only allow the boundary to be the minimal possible for that
/// dimension. Eg, the dimension of a simplex is the size of its boundary minus 1.
///
/// ## Example
///
/// The following file would define the boundary of a tetrahedron.
///
/// 4 # 4 vertices
/// 1 2 # simplex #5 is the first edge
/// 1 3 # 6
/// 1 4 # 7
/// 2 3 # 8
/// 2 4 # 9
/// 3 4 # 10
/// 5 6 8 # 11 - the first triangle 123
/// 5 7 9 # 12 (124)
/// 6 7 10 # 13 (134)
/// 8 9 10 # 14 (234)
/// 11 12 13 14 # the tetrahedron
///
pub fn read_boundary<R: Read>(r: R) -> Result<Vec<Simplex>, Box<dyn std::error::Error>> {
let mut simplices = vec![Simplex::empty()];
let mut lines = BufReader::new(r).lines();
let num_vertices = lines
.next()
.ok_or("The first line should contain the number of vertices, but was empty")??
.parse::<usize>()?;
for i in 0..num_vertices {
simplices.push(Simplex {
j: i,
faces: vec![0],
dim: 0,
});
}
let mut j = simplices.len();
for line in lines.flat_map(|id| id) {
let line: &str = if let Some(i) = line.find('#') {
&line[..i]
} else {
&line
}
.trim();
if line.is_empty() {
continue;
}
let mut v = Vec::new();
for f in line.split_whitespace() {
v.push(f.parse::<usize>()?);
}
let dim = v.len() as isize - 1;
simplices.push(Simplex { j, faces: v, dim });
j += 1;
}
Ok(simplices)
}
/// Read vertex 2D positions, one for each line.
pub fn read_2d_points<R: Read>(r: R) -> Result<Vec<[f64; 2]>, Box<dyn std::error::Error>> {
let lines = BufReader::new(r).lines();
let mut points = Vec::new();
for line in lines {
let line = line?;
let mut sp = line.split_whitespace();
let x = sp.next().ok_or("Missing x in line")?.parse::<f64>()?;
let y = sp.next().ok_or("Missing x in line")?.parse::<f64>()?;
points.push([x, y]);
}
Ok(points)
}
/// Read a point triangulation from `reader`.
///
/// ### File format
/// The file format used is as follows
///
/// number-of-verts
/// point data of the vertices
/// empty line
/// triangle point indices
///
/// ### Note
/// This only reads a 2D triangulation.
pub fn read_triangulation<R: Read>(
r: &mut R,
) -> Result<(Vec<[f64; 2]>, Vec<Simplex>), Box<dyn std::error::Error>> {
let mut lines = BufReader::new(r).lines();
let first_line = lines
.next()
.ok_or("reader is empty in `read_triangulation`")??;
let number_of_vertices: usize = first_line.parse::<usize>()?;
let points = (&mut lines)
.take(number_of_vertices)
.map(|line| {
// TODO: Errors here aren't caught.
let line = line.expect("Need to catch errors in here!");
let mut buf = [0f64; 2];
for (i, string) in line.split(' ').enumerate() {
assert!(
i < 2,
"the line should not contain more than two numbers! (2D only!)"
);
buf[i] = string
.parse::<f64>()
.expect("Need to catch errors in here!");
}
buf
})
.collect::<Vec<[f64; 2]>>();
assert_eq!(
lines
.next()
.expect("Missing empty line in triangulation format!")?,
""
); // the empty line
// Next we read in the face indices.
let mut faces_to_points = Vec::new();
for line in lines {
let mut buf = [0usize; 3];
for (i, string) in line.unwrap().split(' ').enumerate() {
assert!(
i < 3,
"the line should not contain more than three indices (2D only!)"
);
buf[i] = string.parse::<usize>().unwrap();
}
faces_to_points.push(buf);
}
// TOOD: clean this up: it is bad.
let simplices: Vec<Vec<usize>> = {
let mut simplices: Vec<Vec<usize>> =
Vec::with_capacity(points.len() + faces_to_points.len() + 1);
simplices.push(vec![]); // empty simplex of dim `-1`.
for _ in &points {
simplices.push(vec![0]);
}
// Now we want to make the edges, which are implicit in the face data.
// However, we must also avoid double counting. We make a hashmap to find all
// unique edges, and vecs for edges and faces.
let mut edge_to_index = std::collections::HashMap::new();
// To avoid double counting we store the edges sorted in the hashmap.
fn sorted_a2([x, y]: [usize; 2]) -> [usize; 2] {
if x < y {
[x, y]
} else {
[y, x]
}
}
let mut faces: Vec<Vec<usize>> = Vec::new();
for &face in &faces_to_points {
// Winding order seems to be CW, but not sure if this is guaranteed by CGAL.
// Faces are zero-indexed, but we have the empty simlex, so we must shift the indices
// by one.
let current_edges = [
[face[0] + 1, face[1] + 1],
[face[1] + 1, face[2] + 1],
[face[2] + 1, face[0] + 1],
];
let mut buf = [0usize; 3];
for (i, &edge) in current_edges.iter().enumerate() {
let edge_i = *edge_to_index.entry(sorted_a2(edge)).or_insert_with(|| {
let index = simplices.len();
simplices.push(edge.to_vec());
index
});
buf[i] = edge_i;
}
faces.push(buf.to_vec());
}
simplices.extend_from_slice(&faces);
simplices
};
let simplices = simplices
.into_iter()
.enumerate()
.map(|(j, v)| {
let dim = v.len() as isize - 1;
Simplex { j, faces: v, dim }
})
.collect::<Vec<_>>();
Ok((points, simplices))
}
/// Open the file given, or if the filename is "-", open `stdin`.
fn open_file_or_stdin(
filename: &str,
) -> Result<Box<dyn std::io::Read>, Box<dyn std::error::Error>> {
if filename == "-" {
Ok(Box::new(std::io::stdin()))
} else {
Ok(Box::new(File::open(Path::new(filename))?))
}
}
/// Open the file given, or if the filename is "-", open `stdin`.
fn open_file_or_stdout(
filename: &str,
) -> Result<Box<dyn std::io::Write>, Box<dyn std::error::Error>> {
if filename == "-" {
Ok(Box::new(std::io::stdout()))
} else {
Ok(Box::new(File::create(Path::new(filename))?))
}
}
fn main() {
use clap::*;
let matches = clap_app!(
myapp =>
(name: crate_name!())
(about: crate_description!())
(version: crate_version!())
(author: crate_authors!())
// (@arg triangulation: -t --triangulation[FILE] "Input a 2D triangulation in text form.")
(@arg exhaustive: -e --exhaustive "Run the exhaustive variant of the algorithm.")
// (@arg svg: --svg[FILE] 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[FILE] "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 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.
The first line of the boundary file should be the number of points in the complex.
Subsequent lines each contain the boundary of a simplex as a list of simplex IDs. The ID of a simplex
is it's sequence number in the file. See `data/tet.boundary` for an example of a tetrahedron.
"),
).get_matches();
let mut was_useful = false;
let debug = matches.is_present("debug");
let mut simplices = {
let mut reader = open_file_or_stdin(matches.value_of("boundary").unwrap())
.expect("failed to open given file");
read_boundary(&mut reader).expect("Couldn't read boundary")
};
let mut alphas = None;
let mut points = None;
if matches.is_present("alpha") {
let pts = {
let mut reader = open_file_or_stdin(matches.value_of("points").unwrap())
.expect("failed to open given file");
read_2d_points(&mut reader).expect("Couldn't read 2d points")
};
let (ordered_simplices, alphs) = compute_alpha_values_2d(&pts, &simplices);
alphas = Some(alphs);
points = Some(pts);
simplices = ordered_simplices;
}
let persistence = Persistence {
simplices,
points: points.unwrap_or(Vec::new()),
alphas: alphas.unwrap_or(Vec::new()),
};
let variant = if matches.is_present("exhaustive") {
Exhaustive
} else {
Standard
};
let reduction = reduce(&persistence, variant);
if debug {
let mut out = std::io::stdout();
reduction
.debug_output(&mut out, &persistence.simplices)
.unwrap();
was_useful = true;
}
if let Some(path) = matches.value_of("stats") {
let mut stats_file = open_file_or_stdout(path).unwrap();
let json = serde_json::to_string(&reduction.log.aggregate()).expect("failed to convert stats to .json");
write!(stats_file, "{}", json).expect("failed to write stats to file");
was_useful = true;
}
// if let Some(path) = matches.value_of("textable") {
// let mut tex_file = File::create(Path::new(path)).unwrap();
// TODO: FIXME
// output::tex_table(&stats, &stats, &mut tex_file).unwrap();
// }
if let Some(path) = matches.value_of("graphviz") {
let mut f = File::create(Path::new(path)).unwrap();
output::graphviz(&persistence, &reduction, &reduction, &mut f).unwrap();
was_useful = true;
}
if let Some(path) = matches.value_of("diagram") {
let birth_death_pairs = reduction
.pairings
.iter()
.map(|&(a, b)| (persistence.alphas[a], persistence.alphas[b]))
.collect::<Vec<(f64, f64)>>();
output::persistence(&birth_death_pairs, &Path::new(path));
was_useful = true;
}
if let Some(path) = matches.value_of("svg") {
output::svg(&persistence, &reduction, &reduction, &Path::new(path));
was_useful = true;
}
if !was_useful {
eprintln!("It doesn't look like you actually got anything useful here!");
eprintln!("See `cra --help` if you are confused.");
eprintln!("{}", matches.usage());
}
}