From 7c5880309c33df60503735f796186442cbda4955 Mon Sep 17 00:00:00 2001 From: Taylor Henderson Date: Tue, 28 Jan 2025 19:46:41 -0500 Subject: [PATCH] little speed ups --- src/barcode.rs | 4 +- src/birthdeath.rs | 4 +- src/main.rs | 7 +- src/persistencelandscape.rs | 263 ++++++++++++++++++++---------------- src/plot.rs | 2 +- src/rpls.rs | 14 +- 6 files changed, 162 insertions(+), 132 deletions(-) diff --git a/src/barcode.rs b/src/barcode.rs index 14f5193..190cf14 100644 --- a/src/barcode.rs +++ b/src/barcode.rs @@ -64,10 +64,10 @@ impl Eq for Node {} #[derive(Debug, Clone)] struct Event { event_type: EventType, - value: f64, + value: f32, } -const fn create_event(birth: f64, death: f64, i: usize) -> Node { +const fn create_event(birth: f32, death: f32, i: usize) -> Node { Node { birth_event: Event { event_type: EventType::Birth, diff --git a/src/birthdeath.rs b/src/birthdeath.rs index e61b3bb..3ad9068 100644 --- a/src/birthdeath.rs +++ b/src/birthdeath.rs @@ -9,8 +9,8 @@ use std::str::FromStr; #[derive(Debug)] pub struct BirthDeath { - pub birth: f64, - pub death: f64, + pub birth: f32, + pub death: f32, } impl FromStr for BirthDeath { diff --git a/src/main.rs b/src/main.rs index ca06106..958c117 100644 --- a/src/main.rs +++ b/src/main.rs @@ -9,6 +9,7 @@ use clap::Parser; use csv::Writer; use std::error::Error; use std::fs; +use std::time::Instant; /// Generates the PL for a set of birth death pairs #[derive(Parser, Debug)] @@ -40,6 +41,7 @@ struct Args { fn main() -> Result<(), Box> { let args = Args::parse(); + let now = Instant::now(); let bd_paris: Vec = fs::read_to_string(args.name)? .lines() .filter(|s| !s.contains("inf") && !s.is_empty()) @@ -49,6 +51,9 @@ fn main() -> Result<(), Box> { let landscapes = fast_pl::rpls::pairs_to_landscape(bd_paris, args.k, args.debug)?; + let elapsed = now.elapsed(); + println!("Elapsed: {elapsed:.?}"); + if !args.csv.is_empty() { let mut wtr = Writer::from_path(args.csv)?; for landscape in &landscapes { @@ -71,7 +76,7 @@ fn main() -> Result<(), Box> { #[cfg(test)] mod tests { - fn test_runner(k: usize, bd_pairs_vec: Vec<(f64, f64)>, answer_vec: Vec>) { + fn test_runner(k: usize, bd_pairs_vec: Vec<(f32, f32)>, answer_vec: Vec>) { let bd_pairs = bd_pairs_vec .into_iter() .map(|(x, y)| fast_pl::birthdeath::BirthDeath { birth: x, death: y }) diff --git a/src/persistencelandscape.rs b/src/persistencelandscape.rs index 4d59251..50988c9 100644 --- a/src/persistencelandscape.rs +++ b/src/persistencelandscape.rs @@ -13,7 +13,7 @@ use geo::{ use std::cmp::min; use std::collections::{BinaryHeap, VecDeque}; -#[derive(Debug, Clone, Copy)] +#[derive(Debug)]//, Clone, Copy)] struct PersistenceMountain { position: Option, slope_rising: bool, @@ -25,8 +25,8 @@ struct PersistenceMountain { #[derive(Debug, Clone, Copy, PartialEq, Eq)] pub struct PointOrd { - pub x: FloatOrd, - pub y: FloatOrd, + pub x: FloatOrd, + pub y: FloatOrd, } impl Ord for PointOrd { @@ -70,13 +70,13 @@ impl Ord for Event { if self.value != other.value{ return self.value.cmp(&other.value).reverse() } - if self.event_type != other.event_type{ - return self.event_type.cmp(&other.event_type).reverse(); - } - if self.parent_mountain_id != other.parent_mountain_id { - return self.parent_mountain_id.cmp(&other.parent_mountain_id).reverse(); - } - self.parent_mountain2_id.cmp(&other.parent_mountain2_id).reverse() + // if self.event_type != other.event_type{ + self.event_type.cmp(&other.event_type).reverse() + // } + // if self.parent_mountain_id != other.parent_mountain_id { + // return self.parent_mountain_id.cmp(&other.parent_mountain_id).reverse(); + // } + // self.parent_mountain2_id.cmp(&other.parent_mountain2_id).reverse() } } @@ -95,7 +95,7 @@ impl PartialEq for Event { impl Eq for Event {} -fn create_mountain(birth: f64, death: f64, index: usize) -> PersistenceMountain { +fn create_mountain(birth: f32, death: f32, index: usize) -> PersistenceMountain { let half_dist = (death - birth) / 2.0; PersistenceMountain { @@ -123,12 +123,12 @@ fn generate_mountains(bd_pairs: Vec) -> Vec { .filter(|BirthDeath { birth, death }| death.is_finite() && birth.is_finite()) .enumerate() .map(|(i, BirthDeath { birth, death })| create_mountain(birth, death, i)) - .collect::>() + .collect::>() } -fn generate_initial_events(mountains: &Vec) -> Vec { +fn generate_initial_events(mountains: &mut Vec<&mut PersistenceMountain>) -> Vec { mountains - .into_iter() + .iter() .flat_map( |PersistenceMountain { birth, @@ -162,7 +162,7 @@ fn generate_initial_events(mountains: &Vec) -> Vec { .collect() } -const fn current_segment_start(mountain: PersistenceMountain) -> (f64, f64) { +const fn current_segment_start(mountain: &PersistenceMountain) -> (f32, f32) { if mountain.slope_rising { (mountain.birth.x.0, mountain.birth.y.0) } else { @@ -170,7 +170,7 @@ const fn current_segment_start(mountain: PersistenceMountain) -> (f64, f64) { } } -const fn current_segment_end(mountain: PersistenceMountain) -> (f64, f64) { +const fn current_segment_end(mountain: &PersistenceMountain) -> (f32, f32) { if mountain.slope_rising { (mountain.middle.x.0, mountain.middle.y.0) } else { @@ -178,18 +178,18 @@ const fn current_segment_end(mountain: PersistenceMountain) -> (f64, f64) { } } -fn create_line_segment(mountain: PersistenceMountain) -> Line { +fn create_line_segment(mountain: &PersistenceMountain) -> Line { Line { - start: current_segment_start(mountain).into(), - end: current_segment_end(mountain).into(), + start: current_segment_start(&mountain).into(), + end: current_segment_end(&mountain).into(), } } -fn intersects_with_neighbor(m1: PersistenceMountain, m2: PersistenceMountain) -> Option { +fn intersects_with_neighbor(m1: &PersistenceMountain, m2: &PersistenceMountain) -> Option { if m1.slope_rising == m2.slope_rising { return None; } - let inter = line_intersection(create_line_segment(m1), create_line_segment(m2)); + let inter = line_intersection(create_line_segment(&m1), create_line_segment(&m2)); match inter { Some(LineIntersection::SinglePoint { intersection: Coord { x, y }, @@ -201,16 +201,21 @@ fn intersects_with_neighbor(m1: PersistenceMountain, m2: PersistenceMountain) -> // Ignore all colinnear, not proper and no intersection results these will be resolved on // slope change or do not matter Some(i) => match i { - LineIntersection::SinglePoint { intersection: _, is_proper: _ } => None, + LineIntersection::SinglePoint { intersection: _, is_proper: _ } | LineIntersection::Collinear { intersection: _ } => None }, - None => None + _ => None } } +fn float_point_check(p1: PointOrd, p2:PointOrd)-> bool{ + (p1.x.0 - p2.x.0).abs() < f32::EPSILON && + (p1.y.0 - p2.y.0).abs() < f32::EPSILON +} + fn log_to_landscape( - mountain: PersistenceMountain, + mountain: &PersistenceMountain, value: PointOrd, landscapes: &mut [Vec], k: usize, @@ -223,49 +228,56 @@ fn log_to_landscape( // assert_ne!(*landscapes[position].last().unwrap(), value); // } // Ensure points are increasing x (except if points are exactly the same) - if ! landscapes[position].is_empty() && - *landscapes[position].last().unwrap() != value{ - assert!(landscapes[position].last().unwrap().x.0 < value.x.0, - "Last x in landscape {} is {} but new point to be added has an x of {}", - position, - landscapes[position].last().unwrap().x.0, - value.x.0); - } - // Ensure birth/death is in bottom most landscape (exception if the nearest is a tie, they - // are just dieing out of order and the other must die right after) - let below = position + 1; - if value.y.0 == 0.0 && - below < landscapes.len() && - ! landscapes[below].is_empty(){ - if landscapes[below].last().unwrap() == landscapes[position].last().unwrap(){ - // This is fine, ignore. See above comment - } - else{ - // println!("{:?}", landscapes[below].last().unwrap()); - // println!("{:?}", landscapes[position].last().unwrap()); - // println!("{:?}", mountain); - assert!(landscapes[below].last().unwrap().y.0 == 0.0, - "Attempting to add a birth/death ({},{}) to higher landscape {} when {} is non zero ({},{})", - value.x.0, - value.y.0, - position, - below, - landscapes[below].last().unwrap().x.0, - landscapes[below].last().unwrap().y.0, - ); - } - } + // if ! landscapes[position].is_empty(){ + // if float_point_check(*landscapes[position].last().unwrap(), value) { + // // Ignore, this is fine. They are the same + // } + // else{ + // assert!(landscapes[position].last().unwrap().x.0 < value.x.0, + // "Last x in landscape {} is ({},{}) but new point to be added has an x of ({},{})", + // position, + // landscapes[position].last().unwrap().x.0, + // landscapes[position].last().unwrap().y.0, + // value.x.0, + // value.y.0 + // ); + // } + // } + // // Ensure birth/death is in bottom most landscape (exception if the nearest is a tie, they + // // are just dieing out of order and the other must die right after) + // let below = position + 1; + // if value.y.0 == 0.0 && + // below < landscapes.len() && + // ! landscapes[below].is_empty(){ + // if float_point_check(*landscapes[below].last().unwrap(), *landscapes[position].last().unwrap()){ + // // This is fine, ignore. See above comment + // } + // else{ + // // println!("{:?}", landscapes[below].last().unwrap()); + // // println!("{:?}", landscapes[position].last().unwrap()); + // // println!("{:?}", mountain); + // assert!(landscapes[below].last().unwrap().y.0 == 0.0, + // "Attempting to add a birth/death ({},{}) to higher landscape {} when {} is non zero ({},{})", + // value.x.0, + // value.y.0, + // position, + // below, + // landscapes[below].last().unwrap().x.0, + // landscapes[below].last().unwrap().y.0, + // ); + // } + // } landscapes[position].push(value); } } fn find_intersection( status: &VecDeque, - m1: PersistenceMountain, - mountains: &[PersistenceMountain], + parent_mountain_id: usize, + mountains: &mut Vec<&mut PersistenceMountain>, direction_to_check: Direction, ) -> Option { - let position = m1.position.expect("Intersection check for dead mountain"); + let position = mountains[parent_mountain_id].position.expect("Intersection check for dead mountain"); // Stop underflow of unsigned number if position == 0 && direction_to_check == Direction::Above { return None; @@ -276,15 +288,15 @@ fn find_intersection( }; if let Some(neighbor) = status.get(neighbor_index) { - if let Some(intersection) = intersects_with_neighbor(m1, mountains[*neighbor]) { - let intersection = Event { + if let Some(intersection) = intersects_with_neighbor(&mountains[parent_mountain_id], &mountains[*neighbor]) { + return Some( Event { value: intersection, event_type: EventType::Intersection, - parent_mountain_id: m1.id, + parent_mountain_id, parent_mountain2_id: Some(*neighbor), - }; + }); // println!("{intersection:?}"); - return Some(intersection); + // return Some(intersection); } } None @@ -309,69 +321,80 @@ fn handle_up(state: &mut State, event: &Event){ state.mountains[event.parent_mountain_id].position = Some(position); // Add to output if needed log_to_landscape( - state.mountains[event.parent_mountain_id], + &state.mountains[event.parent_mountain_id], event.value, &mut state.landscapes, state.k, ); // Check and handle all intersections let new_event = find_intersection( - & state.status, - state.mountains[event.parent_mountain_id], - &state.mountains, + &state.status, + event.parent_mountain_id, + state.mountains, Direction::Above, ); if let Some(intersection) = new_event{ - handle_intersection(state, &intersection); + handle_intersection(state, intersection); } } -fn handle_intersection(state: &mut State, event: &Event) -> Option{ - let parent_mountain2_id = event - .parent_mountain2_id - .expect("Intersection event with no second mountain"); - // Add to ouput if needed - log_to_landscape( - state.mountains[event.parent_mountain_id], - event.value, - &mut state.landscapes, - state.k, +fn handle_intersection(state: &mut State, event: Event){ + let weird_q = &mut VecDeque::new(); + weird_q.push_back(event); + while ! weird_q.is_empty(){ + let event = weird_q.pop_front().unwrap(); + let parent_mountain2_id = event + .parent_mountain2_id + .expect("Intersection event with no second mountain"); + // Add to ouput if needed + log_to_landscape( + state.mountains[event.parent_mountain_id], + event.value, + &mut state.landscapes, + state.k, ); - log_to_landscape( - state.mountains[parent_mountain2_id], - event.value, - &mut state.landscapes, - state.k + log_to_landscape( + state.mountains[parent_mountain2_id], + event.value, + &mut state.landscapes, + state.k ); - let (lower, upper) = if state.mountains[event.parent_mountain_id].slope_rising {( - state.mountains[event.parent_mountain_id], - state.mountains[parent_mountain2_id]) - } else{( - state.mountains[parent_mountain2_id], - state.mountains[event.parent_mountain_id], - )}; - // Swap - state.status.swap( - upper.position.expect("Dead mountain in intersection event"), - lower.position.expect("Dead mountain in intersection event"), + let lower_id = if state.mountains[event.parent_mountain_id].slope_rising { + event.parent_mountain_id + } else{ + parent_mountain2_id + // ) + }; + let upper_id = if state.mountains[event.parent_mountain_id].slope_rising { + parent_mountain2_id + } else{ + event.parent_mountain_id + }; + // Swap + state.status.swap( + state.mountains[upper_id].position.expect("Dead mountain in intersection event"), + state.mountains[lower_id].position.expect("Dead mountain in intersection event"), ); - assert!(upper.position != lower.position); - (state.mountains[lower.id].position, state.mountains[upper.id].position) = - (upper.position, lower.position); - // Check for intersections - // Must check both ways because of no sorting, intersections can be discovered in both - // directions - if let Some(new_event) = - find_intersection(&state.status, state.mountains[lower.id], &state.mountains, Direction::Above) + assert!(state.mountains[upper_id].position != state.mountains[lower_id].position); + let tmp = state.mountains[lower_id].position; + state.mountains[lower_id].position = state.mountains[upper_id].position; + state.mountains[upper_id].position = tmp; + // Check for intersections + // Must check both ways because of no sorting, intersections can be discovered in both + // directions + if let Some(new_event) = + find_intersection(&state.status, lower_id, state.mountains, Direction::Above) { - handle_intersection(state, &new_event); + // handle_intersection(state, &new_event); + weird_q.push_back(new_event); } - if let Some(new_event) = - find_intersection(&state.status, state.mountains[upper.id], &state.mountains, Direction::Below) + if let Some(new_event) = + find_intersection(&state.status, upper_id, state.mountains, Direction::Below) { - handle_intersection(state, &new_event); + // handle_intersection(state, new_event); + weird_q.push_back(new_event); } - None + } } @@ -389,7 +412,7 @@ fn handle_death(state: &mut State, event: &Event){ // } // Add to ouput if needed log_to_landscape( - state.mountains[event.parent_mountain_id], + & state.mountains[event.parent_mountain_id], event.value, &mut state.landscapes, state.k, @@ -416,7 +439,7 @@ fn handle_down(state: &mut State, event: &Event){ state.mountains[event.parent_mountain_id].slope_rising = false; // Add to ouput if needed log_to_landscape( - state.mountains[event.parent_mountain_id], + & state.mountains[event.parent_mountain_id], event.value, &mut state.landscapes, state.k, @@ -424,21 +447,21 @@ fn handle_down(state: &mut State, event: &Event){ // Check for intersections let new_event = find_intersection( &state.status, - state.mountains[event.parent_mountain_id], - &state.mountains, + event.parent_mountain_id, + state.mountains, Direction::Below, ); if let Some(intersection) = new_event{ - handle_intersection(state, &intersection); + handle_intersection(state, intersection); } } #[derive(Debug)] -struct State{ +struct State<'a>{ status: VecDeque, - mountains: Vec, + mountains: &'a mut Vec<&'a mut PersistenceMountain>, landscapes: Vec>, events: BinaryHeap, k: usize, @@ -449,12 +472,14 @@ struct State{ /// Will panic if invalid state is discovered during generation #[must_use] pub fn generate(bd_pairs: Vec, k: usize, debug: bool) -> Vec> { - let mountains = generate_mountains(bd_pairs); + let mut binding = generate_mountains(bd_pairs); + let mut mountains: Vec<&mut PersistenceMountain> + = binding.iter_mut().collect(); let mut state = State{ - events: BinaryHeap::from(generate_initial_events(&mountains)), + events: BinaryHeap::from(generate_initial_events(&mut mountains)), status: VecDeque::new(), - mountains, + mountains: &mut mountains, landscapes: empty_landscape(k), k, }; @@ -474,7 +499,7 @@ pub fn generate(bd_pairs: Vec, k: usize, debug: bool) -> Vec { handle_death(&mut state, &event); } - _ => unreachable!("Event type should not be here") + EventType::Intersection => unreachable!("Event type should not be here") } } diff --git a/src/plot.rs b/src/plot.rs index 09267be..d9d6281 100644 --- a/src/plot.rs +++ b/src/plot.rs @@ -24,7 +24,7 @@ pub fn landscape( width: u32, ) -> Result<(), Box> { // Set up the data - let to_plot: Vec> = landscape + let to_plot: Vec> = landscape .into_iter() .map(|s| s.into_iter().map(|PointOrd { x, y }| (x.0, y.0)).collect()) .collect(); diff --git a/src/rpls.rs b/src/rpls.rs index b662ef7..4401815 100644 --- a/src/rpls.rs +++ b/src/rpls.rs @@ -18,7 +18,7 @@ use crate::barcode; pub fn pairs_to_landscape(bd_pairs: Vec, k:usize, debug:bool) -> Result>, &'static str>{ let bd_pairs: Vec = bd_pairs .into_iter() - .filter(|bd| (bd.birth - bd.death).abs() > f64::EPSILON) + .filter(|bd| (bd.birth - bd.death).abs() > f32::EPSILON) .collect(); if bd_pairs.is_empty() { if debug { @@ -41,7 +41,7 @@ pub fn pairs_to_landscape(bd_pairs: Vec, k:usize, debug:bool) -> Res Ok(landscape) } -fn area_under_line_segment(a: persistencelandscape::PointOrd, b: persistencelandscape::PointOrd) ->f64 { +fn area_under_line_segment(a: persistencelandscape::PointOrd, b: persistencelandscape::PointOrd) ->f32 { let height = (a.y.0 - b.y.0).abs(); let base = a.x.0 - b.x.0; let triangle = (height * base) / 2.0; @@ -55,12 +55,12 @@ fn area_under_line_segment(a: persistencelandscape::PointOrd, b: persistenceland triangle + rectangle } -fn landscape_norm(landscape: &[persistencelandscape::PointOrd]) -> f64 { +fn landscape_norm(landscape: &[persistencelandscape::PointOrd]) -> f32 { landscape .iter() .zip(landscape.iter().skip(1)) .map(|(&a, &b)| area_under_line_segment(a, b)) - .sum::() + .sum::() } fn is_sorted(data: &[T]) -> bool @@ -74,11 +74,11 @@ where /// /// Will panic if areas are not strictly decreasing or equal #[must_use] -pub fn l2_norm(landscapes: &[Vec]) -> f64 { +pub fn l2_norm(landscapes: &[Vec]) -> f32 { let areas = landscapes .iter() .map(|l| FloatOrd(landscape_norm(l))) - .collect::>>(); + .collect::>>(); assert!(is_sorted(& areas)); areas.iter().map(|x| x.0).sum() @@ -87,6 +87,6 @@ pub fn l2_norm(landscapes: &[Vec]) -> f64 { /// # Errors /// /// Will return 'Err' if failed to compute persistencelandscape from `bd_pairs` -pub fn pairs_to_l2_norm(bd_paris: Vec, k:usize, debug:bool) -> Result{ +pub fn pairs_to_l2_norm(bd_paris: Vec, k:usize, debug:bool) -> Result{ Ok(l2_norm(pairs_to_landscape(bd_paris, k, debug)?.as_slice())) }