~mht/cra

9cc525935909f3a7f0a705a128c9b46943e3706d — Martin Hafskjold Thoresen 7 months ago deef01a
Begin on improved reduction implementation
1 files changed, 76 insertions(+), 1 deletions(-)

M cra/src/persistence.rs
M cra/src/persistence.rs => cra/src/persistence.rs +76 -1
@@ 68,7 68,7 @@ impl Persistence {
pub struct Reduction {
    /// Two simplices are paired when a simplex has been fully reduced; it is then paired with its highest boundary element.
    pub pairings: Vec<(usize, usize)>,
    ///
    /// The original `j` of the simplex stored in `simplex_with_low`.
    pub j_stored_at: Vec<usize>,
    /// This maps an index to the reduced simplex which maximum boundary element is that index.
    pub simplex_with_low: Vec<Option<Vec<usize>>>,


@@ 153,10 153,13 @@ pub fn reduce(p: &Persistence, variant: Variant, stats: &mut Statistics) -> Redu
    // `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 reduced_with_low: Vec<Option<Simplex>> = vec![None; n];

    let mut pairings = Vec::new();
    let mut additions = Vec::new();

    // We're reducing each simplex.
    'simplex_loop:
    for simplex in &p.simplices {
        let j = simplex.j;
        // With no boundary we're done.


@@ 172,6 175,8 @@ pub fn reduce(p: &Persistence, variant: Variant, stats: &mut Statistics) -> Redu
        // TODO: this probably doens't need to be cloned.
        let mut current_simplex = simplex.faces.clone();
        current_simplex.sort();
        let mut current_reduced = simplex.clone();
        current_reduced.faces.sort();

        let mut r_iter = 0;
        loop {


@@ 180,6 185,28 @@ pub fn reduce(p: &Persistence, variant: Variant, stats: &mut Statistics) -> Redu
            let cand = &simplex_with_low[low];
            stats.swl();

            match &mut reduced_with_low[low] {
                place @ &mut None => {
                    *place = Some(current_reduced);
                    continue 'simplex_loop
                }
                &mut Some(ref other) => {
                    debug_assert!(low <= j);
                    additions.push((other.j, j));
                    simplex_reduce(&mut current_reduced, other);

                    stats.col_add(simplex.dim() as usize);
                    stats.add_sizes(current_reduced.faces.len(), other.faces.len());

                    if current_reduced.faces.len() == 0 {
                        stats.zeroed(simplex.dim() as usize);
                        continue 'simplex_loop;
                    }
                }
            }



            if cand.is_some() {
                additions.push((j_stored_at[low], j));
                column_add(&mut current_simplex, cand.as_ref().unwrap(), s_dim, stats);


@@ 269,6 296,53 @@ pub fn reduce(p: &Persistence, variant: Variant, stats: &mut Statistics) -> Redu
    }
}

/// Reduce `this` by joinint 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;
        }
    }

	if fastpath {
        // fast path
        this.faces.reserve(num_inds);
        this.faces.clear();
        this.faces.copy_from_slice(buf);
    } else {
        std::mem::swap(&mut this.faces, &mut slow_vec);
    }
}

/// 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(


@@ 466,3 540,4 @@ pub fn compute_alpha_values_2d(

    (ordered, alphas)
}