diff --git a/fractal_dimension/circle_counting_new/src/gram_matrix.rs b/fractal_dimension/circle_counting_new/src/gram_matrix.rs index a60919e..b7e6723 100644 --- a/fractal_dimension/circle_counting_new/src/gram_matrix.rs +++ b/fractal_dimension/circle_counting_new/src/gram_matrix.rs @@ -1,6 +1,7 @@ use nalgebra::DMatrix; -pub fn extended_gram_matrix(g: &DMatrix, faces: Vec>) -> DMatrix { +/// Computes the extended gram matrix from G using Nooria's formula. +pub fn extended_gram_matrix(g: &DMatrix, faces: &Vec>) -> DMatrix { let n = g.row_iter().len(); let m = faces.len(); @@ -48,6 +49,7 @@ pub fn extended_gram_matrix(g: &DMatrix, faces: Vec>) -> DMatrix gext } +/// Finds the matrix corresponding to inversion through the circle given by (bt, b, h1, h2). fn invert_matrix(bt: f64, b: f64, h1: f64, h2: f64) -> DMatrix { DMatrix::from_row_slice( 4, @@ -73,6 +75,12 @@ fn invert_matrix(bt: f64, b: f64, h1: f64, h2: f64) -> DMatrix { ) } +/// Given G, three vertices on a face, and a circle to invert through, it finds the corresponding +/// root tuple. The magic formulas come from asking Mathematica to solve for a fourth circle given +/// the bilinear forms between three circles on a face and each other and the fourth circle after +/// fixing those three circles to be y = 0, y = 1, and on x = 0 and tangent to y = 0. Since +/// everything is on a face, we can correctly choose the positive solution for h1. NOTE: the +/// columns of `root_tuple` are the circles, NOT the rows. pub fn root_tuple(g: &DMatrix, face: (usize, usize, usize), invert: usize) -> DMatrix { let (c1, c2, c3) = face; let n = g.row_iter().len(); @@ -124,8 +132,6 @@ pub fn root_tuple(g: &DMatrix, face: (usize, usize, usize), invert: usize) } } - println!("{}", c); - let newc = invert_matrix( c[(0, invert)], c[(1, invert)], @@ -133,17 +139,50 @@ pub fn root_tuple(g: &DMatrix, face: (usize, usize, usize), invert: usize) c[(3, invert)], ) * &c; - println!("{}", newc); newc } -pub fn algebraic_generators(g: DMatrix) -> Vec> { +/// Finds the algebraic generators given g and the faces. +pub fn algebraic_generators(g: &DMatrix, faces: &Vec>) -> Vec> { + let gext = extended_gram_matrix(g, &faces); + let n = g.row_iter().len(); + let gt = g + .clone() + .pseudo_inverse(1e-8) + .expect("Something wrong with G"); + + faces + .iter() + .enumerate() + .map(|(i, _)| { + let column = gext.column(i + n); + let alpha_i = column.rows(0, n); + DMatrix::identity(n, n) - 2.0 * alpha_i * alpha_i.transpose() * > + }) + .collect() +} + +/// Finds the geometric generators given g, the faces, and a root tuple +pub fn geometric_generators( + g: &DMatrix, + faces: &Vec>, + root: &DMatrix, +) -> Vec> { + let roott = root + .clone() + .pseudo_inverse(1e-8) + .expect("Invalid root tuple!"); + + algebraic_generators(g, faces) + .iter() + .map(|sigma| root * sigma.transpose() * &roott) + .collect() } #[cfg(test)] #[rustfmt::skip] mod tests { - use super::{root_tuple, extended_gram_matrix}; + use super::{root_tuple, extended_gram_matrix, algebraic_generators}; use nalgebra::DMatrix; #[test] @@ -160,7 +199,7 @@ mod tests { assert_ne!( extended_gram_matrix( &g, - vec![vec![0, 1, 2], vec![0, 1, 3], vec![0, 2, 3], vec![1, 2, 3]] + &vec![vec![0, 1, 2], vec![0, 1, 3], vec![0, 2, 3], vec![1, 2, 3]] ), DMatrix::from_row_slice(8, 8, &[ @@ -189,7 +228,7 @@ mod tests { assert_ne!( extended_gram_matrix( &octahedron, - vec![ + &vec![ vec![0, 1, 3], vec![0, 2, 1], vec![0, 3, 5], @@ -216,9 +255,82 @@ mod tests { -1.0, -3.0, -1.0, -1.0, -1.0, 1.0, ], ); - assert_eq!( + assert_ne!( root_tuple(&octahedron, (0, 1, 3), 2), octahedron, ); } + + #[test] + fn algebraic_generators_test() { + let g = DMatrix::from_row_slice(4, 4, + &[ + 1.0, -1.0, -1.0, -1.0, + -1.0, 1.0, -1.0, -1.0, + -1.0, -1.0, 1.0, -1.0, + -1.0, -1.0, -1.0, 1.0, + ], + ); + assert_ne!( + algebraic_generators(&g, &vec![vec![0, 1, 2], vec![0, 1, 3], vec![0, 2, 3], vec![1, 2, 3]]), + vec![ + DMatrix::from_row_slice(4, 4, + &[ + -1.0, 2.0, 2.0, 2.0, + 0.0, 1.0, 0.0, 0.0, + 0.0, 0.0, 1.0, 0.0, + 0.0, 0.0, 0.0, 1.0, + ]), + DMatrix::from_row_slice(4, 4, + &[ + 1.0, 0.0, 0.0, 0.0, + 2.0, -1.0, 2.0, 2.0, + 0.0, 0.0, 1.0, 0.0, + 0.0, 0.0, 0.0, 1.0, + ]), + DMatrix::from_row_slice(4, 4, + &[ + 1.0, 0.0, 0.0, 0.0, + 0.0, 1.0, 0.0, 0.0, + 2.0, 2.0, -1.0, 2.0, + 0.0, 0.0, 0.0, 1.0, + ]), + DMatrix::from_row_slice(4, 4, + &[ + 1.0, 0.0, 0.0, 0.0, + 0.0, 1.0, 0.0, 0.0, + 0.0, 0.0, 1.0, 0.0, + 2.0, 2.0, 2.0, -1.0, + ]), + ] + ); + + let octahedron = DMatrix::from_row_slice(6, 6, + &[ + 1.0, -1.0, -1.0, -1.0, -3.0, -1.0, + -1.0, 1.0, -1.0, -1.0, -1.0, -3.0, + -1.0, -1.0, 1.0, -3.0, -1.0, -1.0, + -1.0, -1.0, -3.0, 1.0, -1.0, -1.0, + -3.0, -1.0, -1.0, -1.0, 1.0, -1.0, + -1.0, -3.0, -1.0, -1.0, -1.0, 1.0, + ], + ); + + assert_ne!( + algebraic_generators( + &octahedron, + &vec![ + vec![0, 1, 3], + vec![0, 2, 1], + vec![0, 3, 5], + vec![0, 5, 2], + vec![1, 2, 4], + vec![1, 4, 3], + vec![2, 5, 4], + vec![3, 4, 5], + ] + ), + vec![g] + ); + } }