fast_pl/
rpls.rs

1#![warn(
2     clippy::all,
3     clippy::pedantic,
4     clippy::nursery,
5     clippy::cargo,
6 )]
7
8use std::cmp;
9use float_ord::FloatOrd;
10
11use crate::birthdeath::BirthDeath;
12use crate::persistencelandscape;
13use crate::barcode;
14
15/// # Errors
16///
17/// Will return 'Err' if failed to compute persistencelandscape from `bd_pairs`
18pub fn pairs_to_landscape(bd_pairs: Vec<BirthDeath>, k:usize, debug:bool, disable_filter: bool) -> Result<Vec<Vec<(f64,f64)>>, &'static str>{
19    let bd_pairs: Vec<BirthDeath> = bd_pairs
20        .into_iter()
21        .filter(|bd| (bd.birth - bd.death).abs() > f64::EPSILON)
22        .collect();
23    if bd_pairs.is_empty() {
24        return Err("No BirthDeath pairs found in file");
25    }
26
27    if debug {
28        println!("{bd_pairs:?}");
29    }
30    let filtered_pairs = if disable_filter{
31        bd_pairs
32    }
33    else{
34        let filtered_pairs = barcode::filter(bd_pairs, k);
35        if debug {
36            println!("{filtered_pairs:?}");
37        }
38        filtered_pairs
39    };
40    let landscape = persistencelandscape::generate(filtered_pairs, k, debug);
41    if debug {
42        println!("{landscape:?}");
43    }
44    Ok(landscape)
45}
46
47fn area_under_line_segment(a: (f64,f64), b: (f64,f64)) ->f64 {
48    let height = (a.1 - b.1).abs();
49    let base = a.0 - b.0;
50    let triangle = (height * base) / 2.0;
51    assert!(triangle > 0.0);
52
53    let s1 = base;
54    let s2 = cmp::min(FloatOrd(a.1), FloatOrd(b.1)).0;
55    let rectangle = s1 * s2;
56    assert!(rectangle >= 0.0);
57
58    triangle + rectangle
59}
60
61fn landscape_norm(landscape: &[(f64,f64)]) -> f64 {
62    landscape
63        .iter()
64        .zip(landscape.iter().skip(1))
65        .map(|(a, b)| area_under_line_segment(*a, *b))
66        .sum::<f64>()
67}
68
69fn is_sorted<T>(data: &[T]) -> bool
70where
71    T: Ord,
72{
73    data.windows(2).all(|w| w[0] <= w[1])
74}
75
76/// # Panics
77///
78/// Will panic if areas are not strictly decreasing or equal
79#[must_use]
80pub fn l2_norm(landscapes: &[Vec<(f64,f64)>]) -> f64 {
81    let areas = landscapes
82        .iter()
83        .map(|l| FloatOrd(landscape_norm(l)))
84        .collect::<Vec<FloatOrd<f64>>>();
85    assert!(is_sorted(& areas));
86
87        areas.iter().map(|x| x.0).sum()
88}
89
90/// # Errors
91///
92/// Will return 'Err' if failed to compute persistencelandscape from `bd_pairs`
93pub fn pairs_to_l2_norm(bd_paris: Vec<BirthDeath>, k:usize, debug:bool, disable_filter: bool) -> Result<f64, &'static str>{
94    Ok(l2_norm(pairs_to_landscape(bd_paris, k, debug, disable_filter)?.as_slice()))
95}