~mht/cra

ref: 71377f32979f2152b84f18f8d6c87a86093634e4 cra/src/persistence.rs -rw-r--r-- 12.7 KiB
71377f32 — Martin Hafskjold Thoresen Fix bug in `reduce` 1 year, 5 months 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
use crate::{log, Point};
use std::io::Write;

/// The three variants of the algorithm.
#[derive(PartialEq, Debug, Clone, Copy)]
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,
}
use self::Variant::*;

/// One simplex
#[derive(Debug, Clone)]
pub struct Simplex {
    /// Original index. We need this for bookkeeping.
    pub j: usize,
    /// Indices of the faces of this simplex
    pub faces: Vec<usize>,
    // The ball radius when this simplex is born.
    // pub r_when_born: f64,
    pub dim: isize,
}

impl Simplex {
    /// Get the dimension of this simplex.
    pub fn dim(&self) -> isize {
        self.dim
    }
    /// Make the empty simplex
    pub fn empty() -> Self {
        Simplex {
            j: 0,
            faces: Vec::new(),
            dim: -1,
        }
    }
}

/// Data we use for stuff.
pub struct Persistence {
    /// All simplices in the data set
    pub simplices: Vec<Simplex>,
    /// Point location data, for rendering.
    pub points: Vec<[f64; 2]>,
    /// Alpha values for all simplices.
    /// TODO: reconsider having this here; what if we don't have any alpha?
    pub alphas: Vec<f64>,
}

/// The result from reducing.
#[derive(Debug)]
pub struct Reduction {
    /// When we reduce a simplex we might finding that it's `low` is unique (that `simplex_with_low[low] == None`).
    /// We then store the reduced simplex in that slot. This corresponds to pairing `low` with `j`.
    pub pairings: Vec<(usize, usize)>,
    /// This maps an index to the reduced simplex which maximum boundary element is that index.
    pub simplex_with_low: Vec<Option<Simplex>>,
    /// Here we track which elements are added to which columns.
    /// This was probably used for the output?
    pub additions: Vec<(usize, usize)>,

    pub log: log::Reduction,
}

impl Reduction {
    /// Write debug output for this reduction to `w`. This includes how each simplex was reduced.
    pub fn debug_output<W: Write>(
        &self,
        w: &mut W,
        simplices: &[Simplex],
    ) -> Result<(), std::io::Error> {
        use std::collections::HashMap;

        let mut added_to = HashMap::<usize, Vec<usize>>::new();
        for &(left, right) in &self.additions {
            added_to.entry(left).or_insert(Vec::new()).push(right);
        }

        let mut killed = HashMap::new();
        for &(victim, killer) in &self.pairings {
            killed.insert(killer, victim);
        }

        writeln!(w, "j;faces;adds;reduced")?;
        for s in simplices {
            let empty_vec = Vec::new();
            let j = s.j;
            let initial_faces = &s.faces;
            let reduced_faces = if let Some(&place) = killed.get(&j) {
                let v = self.simplex_with_low[place]
                    .as_ref()
                    .expect("Simplex was killing, but wasn't stored?");
                &v.faces
            } else {
                &empty_vec
            };
            let adds = added_to.get(&s.j).unwrap_or(&empty_vec);
            writeln!(
                w,
                "{};{:?};{:?};{:?}",
                j, initial_faces, adds, reduced_faces
            )?;
        }

        Ok(())
    }
}

/// 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) -> Reduction {
    let mut pairings = Vec::new();
    let mut additions = Vec::new();
    let mut reduced_with_low: Vec<Option<Simplex>> = vec![None; p.simplices.len()];

    let mut reduction_log = log::Reduction::new();

    // We're reducing each simplex.
    for simplex in &p.simplices {
        // we'll store an index here, if we don't fully reduce the simplex
        let mut place_index = None;
        // With no boundary we're done.
        if simplex.faces.last().is_none() {
            continue;
        }
        let mut current_reduced = simplex.clone();
        current_reduced.faces.sort();

        let mut simplex_log = log::SimplexReduction::new(&simplex);

        while current_reduced.faces.len() > 0 {
            let low = *current_reduced.faces.last().unwrap();
            if let Some(ref other) = reduced_with_low[low] {
                additions.push((current_reduced.j, other.j));
                simplex_log.add(log::Add::regular(&current_reduced, other));
                simplex_reduce(&mut current_reduced, other);
            } else {
                place_index = Some(low);
                pairings.push((low, simplex.j));
                simplex_log.place_index = place_index;
                break;
            }
        }

        if variant == Exhaustive {
            // Reduce the rest of the simplex. Go backwards. We can skip the last index, since
            // that's us. Then we count the numbers of boundary simplices we know we can't reduce
            // any further. Since we're going bottom up and looking at `reduced_with_low`, these
            // cannot change with later redutions.
            let mut safe_simplices = 1;
            let mut n = current_reduced.faces.len();
            while safe_simplices < n {
                let cand_i = current_reduced.faces[n - 1 - safe_simplices];
                if let Some(ref cand) = reduced_with_low[cand_i] {
                    additions.push((current_reduced.j, cand.j));
                    simplex_log.add(log::Add::exhaustive(&current_reduced, cand));
                    simplex_reduce(&mut current_reduced, cand);
                    n = current_reduced.faces.len();
                } else {
                    safe_simplices += 1;
                }
            }
        }

        if current_reduced.faces.len() != 0 {
            let i = place_index.expect(
                "Reduced simplex isn't empty, yet we haven't set the index in which to store it",
            );
            reduced_with_low[i] = Some(current_reduced);
        }

        reduction_log.add(simplex_log);
    }

    Reduction {
        pairings,
        simplex_with_low: reduced_with_low,
        additions,
        log: reduction_log,
    }
}

/// Reduce `this` by joining it with `other`.
fn simplex_reduce(this: &mut Simplex, other: &Simplex) {
    // TODO: This is terrible!
    let mut slow_vec = Vec::new();
    let mut _buffer = [0usize; 32];
    let num_inds = this.faces.len() + other.faces.len();
    let fastpath = num_inds <= 32;

    let mut buf: &mut [usize] = &mut _buffer;
    if !fastpath {
        // We might need more space. Allocate.
        // TODO: this is awkward
        slow_vec = vec![0usize; num_inds];
        buf = &mut slow_vec;
    }

    let mut our_i = 0;
    let mut their_i = 0;
    let mut buffer_i = 0;
    let our_last = this.faces.len();
    let their_last = other.faces.len();

    while our_i < our_last && their_i < their_last {
        if this.faces[our_i] < other.faces[their_i] {
            buf[buffer_i] = this.faces[our_i];
            our_i += 1;
            buffer_i += 1;
        } else if this.faces[our_i] == other.faces[their_i] {
            our_i += 1;
            their_i += 1;
        } else {
            buf[buffer_i] = other.faces[their_i];
            their_i += 1;
            buffer_i += 1;
        }
    }

    while our_i < our_last {
        buf[buffer_i] = this.faces[our_i];
        our_i += 1;
        buffer_i += 1;
    }

    while their_i < their_last {
        buf[buffer_i] = this.faces[their_i];
        their_i += 1;
        buffer_i += 1;
    }

    if fastpath {
        this.faces.reserve(num_inds);
        this.faces.clear();
        for &n in &buf[..buffer_i] {
            this.faces.push(n);
        }
    } else {
        std::mem::swap(&mut this.faces, &mut slow_vec);
    }
}

/// Compute the alpha values for the simplices, given the point location.
/// As the name suggests, this is only for dimension 2 and less.
pub fn compute_alpha_values_2d(
    points: &[[f64; 2]],
    simplices: &[Simplex],
) -> (Vec<Simplex>, Vec<f64>) {
    fn alpha_1d(simplex: &Simplex, points: &[[f64; 2]]) -> f64 {
        let p0 = Point(points[simplex.faces[0]]);
        let p1 = Point(points[simplex.faces[1]]);
        let diff = p1 - p0;
        0.5 * diff.len()
    }

    fn alpha_2d(
        simplex: &Simplex,
        points: &[[f64; 2]],
        simplices: &[Simplex],
    ) -> (f64, Option<usize>) {
        let xy_i = simplex.faces[0];
        let xy = &simplices[xy_i];
        let yz_i = simplex.faces[1];
        let _yz = &simplices[yz_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.
        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 x = Point(points[x_i - 1]);
        let y = Point(points[y_i - 1]);
        let z = Point(points[z_i - 1]);

        let xy = (y - x).len();
        let xz = (z - x).len();
        let yz = (z - y).len();

        let ax = (y - x).angle(z - x);
        let ay = (z - y).angle(x - y);
        let az = (x - z).angle(y - z);
        let am = ax.max(ay).max(az);

        // The circumcenter is the point in which the balls from each vertex meets, so
        // at `r=circumradius` the face is filled in.
        let circumradius = (xy * yz * xz)
            / ((xy + yz + xz) * (xy + yz - xz) * (xy - yz + xz) * (-xy + yz + xz)).sqrt();

        let mut edit = None;
        if am * 180.0 / std::f64::consts::PI >= 90.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 {
                edit = Some(yz_i);
            } else if am == ay {
                edit = Some(zx_i);
            } else {
                edit = Some(xy_i);
            }
        }
        (circumradius, edit)
    }

    #[derive(PartialEq)]
    struct CmpSimplex {
        dim: usize,
        alpha: f64,
        index: usize,
    }

    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 => {
                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),
        };
        v.push(CmpSimplex {
            dim: simplex.faces.len(),
            alpha,
            index: current_simplex_i,
        });
    }

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

    // Now we have simplices in a new order. Next we need to rewrite the faces of all simplices,
    // since the indices are in the old ordering. To do this we build an index mapping from old to new
    // index.
    let mut old_to_new = vec![0; simplices.len()];
    for (i, s) in ordered.iter().enumerate() {
        old_to_new[s.j] = i;
    }
    // Now we rewrite all faces.
    for (i, s) in ordered.iter_mut().enumerate() {
        s.j = i;
        for j in s.faces.iter_mut() {
            *j = old_to_new[*j];
        }
    }

    (ordered, alphas)
}