~mht/cra

1febdad4078b6ad8c87ebe9d43fad8bf0ed943da — Martin Hafskjold Thoresen 1 year, 3 months ago 63085b4
Move alpha values to `Persistence`

This makes more sense, but it's not optimal still.
3 files changed, 64 insertions(+), 59 deletions(-)

M cra/src/main.rs
M cra/src/output.rs
M cra/src/persistence.rs
M cra/src/main.rs => cra/src/main.rs +14 -28
@@ 289,11 289,7 @@ pub fn read_triangulation<R: Read>(
    let simplices = simplices
        .into_iter()
        .enumerate()
        .map(|(j, v)| Simplex {
            j,
            faces: v,
            r_when_born: 0.0,
        })
        .map(|(j, v)| Simplex { j, faces: v })
        .collect::<Vec<_>>();

    Ok((points, simplices))


@@ 303,8 299,12 @@ pub fn read_triangulation<R: Read>(
pub fn read_input_stdin2() -> Result<Persistence, Box<dyn std::error::Error>> {
    let mut stdin = std::io::stdin();
    let (points, mut simplices) = read_triangulation(&mut stdin)?;
    compute_alpha_values_2d(&points, &mut simplices);
    Ok(Persistence { simplices, points })
    let (simplices, alphas) = compute_alpha_values_2d(&points, &simplices);
    Ok(Persistence {
        simplices,
        points,
        alphas,
    })
}

fn main() {


@@ 334,8 334,12 @@ fn main() {
        };
        let (points, mut simplices) =
            read_triangulation(&mut reader).expect("failed to read triangulation");
        compute_alpha_values_2d(&points, &mut simplices);
        Persistence { simplices, points }
        let (simplices, alphas) = compute_alpha_values_2d(&points, &simplices);
        Persistence {
            simplices,
            points,
            alphas,
        }
    } else {
        // Need to take some other format, then!
        unimplemented!();


@@ 343,19 347,6 @@ fn main() {

    // TODO: this should be done somewhere else!
    {
        // Fix the ordering of the simplices!
        persistence.simplices.sort_by(|a, b| {
            if b.faces.contains(&a.j) {
                Less
            } else if a.faces.contains(&b.j) {
                Greater
            } else if f64_eq(a.r_when_born, b.r_when_born) {
                a.dim().cmp(&b.dim())
            } else {
                a.r_when_born.partial_cmp(&b.r_when_born).unwrap()
            }
        });

        // Map j to sorted index.
        let sorted_index_of_j = {
            let mut v = (0..persistence.simplices.len()).collect::<Vec<_>>();


@@ 466,12 457,7 @@ fn main() {
        let birth_death_pairs = standard_reduction
            .pairings
            .iter()
            .map(|&(a, b)| {
                (
                    persistence.simplices[a].r_when_born,
                    persistence.simplices[b].r_when_born,
                )
            })
            .map(|&(a, b)| (persistence.alphas[a], persistence.alphas[b]))
            .collect::<Vec<(f64, f64)>>();
        output::persistence(&birth_death_pairs, &Path::new(path));
    }

M cra/src/output.rs => cra/src/output.rs +28 -17
@@ 39,24 39,33 @@ fn thousand_split(mut num: usize) -> String {
    out
}

/// Write an `svg` file to the file `out_path` showing the simplicial complex.
///
pub fn svg(persistence: &Persistence, std: &Reduction, exh: &Reduction, out_path: &Path) {
    // Whether to draw text on all simplices outputted.
    const DRAW_TEXT: bool = false;
    let mut f = std::io::BufWriter::new(std::fs::File::create(out_path).unwrap());

    // The padding we use for the .svg, so that we know text will be in the image.
    const PADDING: f64 = 50f64;

    let mut min_x = persistence.points[0][0];
    let mut max_x = persistence.points[0][0];
    let mut min_y = persistence.points[0][1];
    let mut max_y = persistence.points[0][1];
    for &point in &persistence.points {
        min_x = min_x.min(point[0]);
        min_y = min_y.min(point[1]);
        max_x = max_x.max(point[0]);
        max_y = max_y.max(point[1]);
    }
    let mut f = std::io::BufWriter::new(std::fs::File::create(out_path).unwrap());

    let min_x = persistence
        .points
        .iter()
        .fold(std::f64::MAX, |m, a| m.min(a[0]));
    let min_y = persistence
        .points
        .iter()
        .fold(std::f64::MAX, |m, a| m.min(a[1]));
    let max_x = persistence
        .points
        .iter()
        .fold(std::f64::MAX, |m, a| m.max(a[0]));
    let max_y = persistence
        .points
        .iter()
        .fold(std::f64::MAX, |m, a| m.max(a[1]));

    let xr = max_x - min_x;
    let yr = max_y - min_y;



@@ 72,7 81,9 @@ pub fn svg(persistence: &Persistence, std: &Reduction, exh: &Reduction, out_path
    let xo = PADDING - min_x;
    let yo = PADDING - min_y;

    for simplex in persistence.simplices.iter() {
    let alpha = &persistence.alphas;

    for (simplex_i, simplex) in persistence.simplices.iter().enumerate() {
        if simplex.dim() == 2 {
            let xy_i = simplex.faces[0];
            let xy = &persistence.simplices[xy_i];


@@ 100,8 111,8 @@ pub fn svg(persistence: &Persistence, std: &Reduction, exh: &Reduction, out_path
            let p1 = Point(persistence.points[y_i - 1]);
            let p2 = Point(persistence.points[z_i - 1]);

            if simplex.r_when_born >= R_CUTOFF {
                if xy.r_when_born < R_CUTOFF {
            if alpha[simplex_i] >= R_CUTOFF {
                if alpha[xy_i] < R_CUTOFF {
                    write!(
                        f,
                        r#"<line x1="{}" y1="{}" x2="{}" y2="{}"


@@ 125,7 136,7 @@ pub fn svg(persistence: &Persistence, std: &Reduction, exh: &Reduction, out_path
                    }
                }

                if zx.r_when_born < R_CUTOFF {
                if alpha[zx_i] < R_CUTOFF {
                    write!(
                        f,
                        r#"<line x1="{}" y1="{}" x2="{}" y2="{}"


@@ 149,7 160,7 @@ pub fn svg(persistence: &Persistence, std: &Reduction, exh: &Reduction, out_path
                    }
                }

                if yz.r_when_born < R_CUTOFF {
                if alpha[yz_i] < R_CUTOFF {
                    write!(
                        f,
                        r#"<line x1="{}" y1="{}" x2="{}" y2="{}"

M cra/src/persistence.rs => cra/src/persistence.rs +22 -14
@@ 39,6 39,9 @@ pub struct Persistence {
    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>,
}

impl Persistence {


@@ 362,15 365,18 @@ pub fn compute_alpha_values_2d(
        0.5 * diff.len()
    }

    fn alpha_2d(simplex: &Simplex, points: &[[f64; 2]], simplices: &[Simplex]) -> f64 {
        let simpl = |i| unsafe { simplex_mut(i, simplices) };

    fn alpha_2d(
        simplex: &Simplex,
        points: &[[f64; 2]],
        simplices: &[Simplex],
        alphas: &mut [f64],
    ) -> f64 {
        let xy_i = simplex.faces[0];
        let xy = simpl(xy_i);
        let xy = &simplices[xy_i];
        let yz_i = simplex.faces[1];
        let _yz = simpl(yz_i);
        let _yz = &simplices[yz_i];
        let zx_i = simplex.faces[2];
        let zx = simpl(zx_i);
        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] {


@@ 412,13 418,13 @@ pub fn compute_alpha_values_2d(
            // this edge.
            if am == ax {
                // debug_assert_eq!(yz / 2.0, simpl(yz_i).r_when_born);
                simpl(yz_i).r_when_born = circumradius;
                alphas[yz_i] = circumradius;
            } else if am == ay {
                // debug_assert_eq!(xz / 2.0, simpl(zx_i).r_when_born);
                simpl(zx_i).r_when_born = circumradius;
                alphas[zx_i] = circumradius;
            } else {
                // debug_assert_eq!(xy / 2.0, simpl(xy_i).r_when_born);
                simpl(xy_i).r_when_born = circumradius;
                alphas[xy_i] = circumradius;
            }
        }
        circumradius


@@ 431,14 437,14 @@ pub fn compute_alpha_values_2d(
    }

    let mut tuples = Vec::new();
    let mut alphas = vec![-1.0; simplices.len()];

    for current_simplex_i in 0..simplices.len() {
        let (former, latter) = simplices.split_at_mut(current_simplex_i);
        let simplex = &mut latter[0];
        let simplex = &simplices[current_simplex_i];
        let alpha = match simplex.dim() {
            -1 | 0 => 0.0,
            1 => alpha_1d(&simplex, &points),
            2 => alpha_2d(&simplex, &points, former),
            2 => alpha_2d(&simplex, &points, simplices, &mut alphas),
            more => panic!("Only support 2d simplices! ({} > 2)", more),
        };
        tuples.push((alpha, current_simplex_i));


@@ 447,10 453,12 @@ pub fn compute_alpha_values_2d(
    tuples.sort_by(|(f1, i1), (f2, i2)| f1.partial_cmp(f2).unwrap_or(i1.cmp(i2)));

    let mut ordered_simplices = Vec::new();
    for &(_, i) in tuples.iter() {
    for &(alpha, i) in tuples.iter() {
        ordered_simplices.push(simplices[i].clone());
        if alphas[i] == -1.0 {
            alphas[i] = alpha;
        }
    }
    let alphas = tuples.into_iter().map(|t| t.0).collect::<Vec<_>>();

    (ordered_simplices, alphas)
}