~mht/cra

9b5383c95502841a37b5e526dc770a7957aeee5a — Martin Hafskjold Thoresen 11 months ago 9cc5259
Greatly simplify the `reduce` method.
2 files changed, 54 insertions(+), 149 deletions(-)

M cra/src/output.rs
M cra/src/persistence.rs
M cra/src/output.rs => cra/src/output.rs +6 -6
@@ 217,9 217,9 @@ pub fn svg(persistence: &Persistence, std: &Reduction, exh: &Reduction, out_path
        if std.simplex_with_low[i].is_some() {
            let s = std.simplex_with_low[i].as_ref().unwrap();
            let e = exh.simplex_with_low[i].as_ref().unwrap();
            if s.len() == 2 {
                std_jumps[i] = s[0];
                exh_jumps[i] = e[0];
            if s.faces.len() == 2 {
                std_jumps[i] = s.faces[0];
                exh_jumps[i] = e.faces[0];
            }
        }
    }


@@ 669,9 669,9 @@ pub fn graphviz<W: Write>(
            let s = std.simplex_with_low[i].as_ref().unwrap();
            let e = exh.simplex_with_low[i].as_ref().unwrap();
            assert_eq!(std.j_stored_at[i], exh.j_stored_at[i]);
            if s.len() == 2 {
                writeln!(w, "  {} -> {} [color=\"red\"]", s[1], s[0])?;
                writeln!(w, "  {} -> {} [color=\"green\"]", e[1], e[0])?;
            if s.faces.len() == 2 {
                writeln!(w, "  {} -> {} [color=\"red\"]", s.faces[1], s.faces[0])?;
                writeln!(w, "  {} -> {} [color=\"green\"]", e.faces[1], e.faces[0])?;
            }
        }
    }

M cra/src/persistence.rs => cra/src/persistence.rs +48 -143
@@ 64,14 64,14 @@ impl Persistence {
}

/// The result from reducing.
#[derive(Debug, PartialEq)]
#[derive(Debug)]
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>>>,
    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)>,


@@ 145,154 145,60 @@ impl Reduction {

/// 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, stats: &mut Statistics) -> Reduction {
    let t0 = time::precise_time_ns();
    let n = p.simplices.len();
    // `simplex_with_low` maps the highest index `m` to the simplex containing it.
    let mut simplex_with_low: Vec<Option<Vec<usize>>> = vec![None; n];
    // `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();

pub fn reduce(p: &Persistence, variant: Variant, _stats: &mut Statistics) -> Reduction {
    let mut reduced_with_low: Vec<Option<Simplex>> = vec![None; p.simplices.len()];
    // We're reducing each simplex.
    'simplex_loop:
    for simplex in &p.simplices {
        let j = simplex.j;
        let mut place_index = None;
        // With no boundary we're done.
        if simplex.faces.last().is_none() {
            continue;
        }
        let s_dim = simplex.dim().max(0) as usize;
        //        // If this simplex is born too late we also skip it.
        //        if simplex.r_when_born >= R_CUTOFF {
        //            stats.skip(s_dim);
        //            continue;
        //        }
        // 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 {
            r_iter += 1;
            let low: usize = *current_simplex.last().unwrap();
            let cand = &simplex_with_low[low];
            stats.swl();

            match &mut reduced_with_low[low] {
                place @ &mut None => {
                    *place = Some(current_reduced);
                    continue 'simplex_loop
        while current_reduced.faces.len() > 0 {
            let low: usize = *current_reduced.faces.last().unwrap();
            match reduced_with_low[low] {
                None => {
                    place_index = Some(low);
                    break;
                }
                &mut Some(ref other) => {
                    debug_assert!(low <= j);
                    additions.push((other.j, j));
                Some(ref other) => {
                    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);
                if current_simplex.is_empty() {
                    stats.zeroed(s_dim);
                    break;
                }
            } else {
                debug_assert!(low <= j);
                pairings.push((low, j));

                stats.placed_full(s_dim);
                simplex_with_low[low] = Some(current_simplex);
                j_stored_at[low] = j;

                if variant == Exhaustive || (variant == ExhaustiveWhenD1 && simplex.dim() == 1) {
                    // We want to reduce the column `j` with `low(j) = low`. Try to remove all indices in this
                    // column.
                    let mut iter = 0;
                    'search: loop {
                        let this = simplex_with_low[low].as_ref().unwrap();
                        let list_len = this.len();
                        for (iter_i, face_i) in (0..list_len - 1).rev().enumerate() {
                            let this_index = this[face_i];
                            let other = &simplex_with_low[this_index];
                            stats.swl();
                            if other.is_none() {
                                // nothing to see here
                                continue;
                            }
                            let this_mut: &mut Vec<usize> = unsafe {
                                // We are indexing in twice into the same array, once mutably, but this is fine,
                                // since the indices are not the same.
                                debug_assert!(this_index != low);
                                let ptr: *const Option<Vec<usize>> =
                                    simplex_with_low.as_ptr().offset(low as isize);
                                let mptr = ptr as *mut Option<Vec<usize>>;
                                (*mptr).as_mut().unwrap()
                            };

                            additions.push((j_stored_at[this_index], j));
                            column_add(this_mut, other.as_ref().unwrap(), s_dim, stats);
                            iter += 1;
                            stats.ex_search(iter_i + 1);
                            continue 'search;
                        }
                        stats.ex_search(list_len);
                        stats.ex_reductions(iter);
                        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 faces = current_reduced.faces.len();
            while faces > safe_simplices {
                let cand_i = current_reduced.faces[faces - 1 - safe_simplices];
                if let Some(ref cand) = reduced_with_low[cand_i] {
                    simplex_reduce(&mut current_reduced, cand);
                    faces = current_reduced.faces.len();
                } else {
                    safe_simplices += 1;
                }
                break;
            }
        }
        stats.reductions(r_iter);

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

    // Find the essential cycles: these are entries that did not give death, and are `null` in the
    // `simplex_with_low` (in matrix terms, this means that row `i` has no `low`s in it).

    // for &cyc in zeroed.difference(&killed) {
    //     let s = p.simplices.iter().find(|&s| s.j == cyc).unwrap();
    //     eprintln!("[reduce] Essential cycle: {:?}", s);
    // }

    // for (i, ptr) in simplex_with_low.iter().enumerate() {
    //     eprint!("{:2}: ", i);
    //     if ptr.is_null() {
    //         eprintln!("null");
    //     } else if ptr.is_implicit() {
    //         eprintln!("{}", i);
    //     } else {
    //         for i in ptr.get().0.iter() {
    //             eprint!("{} ", i);
    //         }
    //         eprintln!("");
    //     }
    // }

    let t1 = time::precise_time_ns();
    stats.time = t1 - t0;
    Reduction {
        pairings,
        j_stored_at,
        simplex_with_low,
        additions,
        pairings: Vec::new(),
        j_stored_at: Vec::new(),
        simplex_with_low: reduced_with_low,
        additions: Vec::new(),
    }
}



@@ 300,17 206,17 @@ pub fn reduce(p: &Persistence, variant: Variant, stats: &mut Statistics) -> Redu
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 _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;


@@ 333,7 239,7 @@ fn simplex_reduce(this: &mut Simplex, other: &Simplex) {
        }
    }

	if fastpath {
    if fastpath {
        // fast path
        this.faces.reserve(num_inds);
        this.faces.clear();


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

    (ordered, alphas)
}