~mht/cra

ref: f53bbcd49cbee65e0810cdee1eddbd785c5ced4c cra/src/main.rs -rw-r--r-- 12.6 KiB
f53bbcd4 — Martin Hafskjold Thoresen Move code to top level 1 year, 3 days ago
                                                                                
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
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());
    }
}