working! want to clean up a bit, though

fractal_dimension/circle_counting_new/data/cube.txt

@@ -1,30 +1,9 @@
-{{{1, 0, 0, 0, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0, 0, 0}, {2, 1, 0, 0, 0, 3, -1,
+{{1, -1, -3, -1, -1, -3, -5, -3}, {-1, 
-0}, {3, 0, 0, 0, 0, 3, -1, 0}, {1, -1, 0, 0, 0, 1, 0, 0}, {0, 0, 0, 0, 0, 1, 0,
+  1, -1, -3, -3, -1, -3, -5}, {-3, -1, 
-0}, {2, 0, 0, 0, 0, 4, -1, 0}, {3, -1, 0, 0, 0, 4, -1, 0}}, {{1, 0, 0, 0, 0, 0,
+  1, -1, -5, -3, -1, -3}, {-1, -3, -1, 
-0, 0}, {0, 1, 0, 0, 0, 0, 0, 0}, {0, 0, 1, 0, 0, 0, 0, 0}, {1, -1, 1, 0, 0, 0,
+  1, -3, -5, -3, -1}, {-1, -3, -5, -3, 
-0, 0}, {3, 1, 2, 0, 0, -1, 0, 0}, {2, 2, 2, 0, 0, -1, 0, 0}, {2, 1, 3, 0, 0,
+  1, -1, -3, -1}, {-3, -1, -3, -5, -1, 
--1, 0, 0}, {3, 0, 3, 0, 0, -1, 0, 0}}, {{1, 0, 0, 0, 0, 0, 0, 0}, {3, 0, 0, 0,
+  1, -1, -3}, {-5, -3, -1, -3, -3, -1, 
-0, 0, -1, 3}, {2, 0, 0, 1, 0, 0, -1, 3}, {0, 0, 0, 1, 0, 0, 0, 0}, {1, 0, 0,
+  1, -1}, {-3, -5, -3, -1, -1, -3, -1, 1}}
--1, 0, 0, 0, 1}, {3, 0, 0, -1, 0, 0, -1, 4}, {2, 0, 0, 0, 0, 0, -1, 4}, {0, 0,
-0, 0, 0, 0, 0, 1}}, {{-1, 4, 0, 0, 0, 0, 2, 0}, {0, 1, 0, 0, 0, 0, 0, 0}, {0,
-1, 0, 0, 0, -1, 1, 0}, {-1, 4, 0, 0, 0, -1, 3, 0}, {-1, 3, 0, 0, 0, 1, 2, 0},
-{0, 0, 0, 0, 0, 1, 0, 0}, {0, 0, 0, 0, 0, 0, 1, 0}, {-1, 3, 0, 0, 0, 0, 3, 0}},
-{{0, 0, 0, 3, 0, -1, 3, 0}, {0, 0, 0, 3, 0, -1, 4, -1}, {0, 0, 0, 1, 0, 0, 1,
--1}, {0, 0, 0, 1, 0, 0, 0, 0}, {0, 0, 0, 2, 0, -1, 3, 1}, {0, 0, 0, 2, 0, -1,
-4, 0}, {0, 0, 0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, 0, 0, 0, 1}}, {{0, 0, 0, -1, 0,
-3, -1, 4}, {0, 0, 0, -1, 0, 3, 0, 3}, {0, 0, 0, -1, 0, 2, 1, 3}, {0, 0, 0, -1,
-0, 2, 0, 4}, {0, 0, 0, 0, 0, 1, -1, 1}, {0, 0, 0, 0, 0, 1, 0, 0}, {0, 0, 0, 0,
-0, 0, 1, 0}, {0, 0, 0, 0, 0, 0, 0, 1}}}
-
-{{0., 8.485281374238571, -1., 0.}, {0., -2.8284271247461903, 
-  1., 0.}, {1.4142135623730951, 18.38477631085024, -5., 
-  1.4142135623730951}, {1.4142135623730951, 
-  7.0710678118654755, -3., 
-  1.4142135623730951}, {1.4142135623730951, 
-  18.38477631085024, -5., -1.4142135623730951}, 
-{1.4142135623730951, 
-  7.0710678118654755, -3., -1.4142135623730951}, 
-{2.8284271247461903, 16.970562748477143, -7., 
-  0.}, {2.8284271247461903, 28.284271247461902, -9., 0.}}
 
 {{0, 1, 5, 4}, {0, 3, 2, 1}, {0, 4, 7, 3}, {1, 2, 6, 5}, {2, 3, 7, 6}, {4, 5, 6, 7}}

fractal_dimension/circle_counting_new/data/octahedron.txt

@@ -1,28 +1,5 @@
-{{{1., 0., 0., 0., 0., 0.}, {0., 1., 0., 0., 0., 0.}, {4., 3., 0., 
+{{1, -1, -1, -1, -1, -3}, {-1, 1, -1, -1, -3, -1}, {-1, -1, 
-   3., -1., 0.}, {0., 0., 0., 1., 0., 0.}, {4., 2., 0., 4., -1., 
+  1, -3, -1, -1}, {-1, -1, -3, 1, -1, -1}, {-1, -3, -1, -1, 
-   0.}, {3., 3., 0., 4., -1., 0.}}, {{1., 0., 0., 0., 0., 0.}, {0., 
+  1, -1}, {-3, -1, -1, -1, -1, 1}}
-   1., 0., 0., 0., 0.}, {0., 0., 1., 0., 0., 0.}, {3., 4., 3., 0., 
-   0., -1.}, {3., 3., 4., 0., 0., -1.}, {2., 4., 4., 0., 
-   0., -1.}}, {{1., 0., 0., 0., 0., 0.}, {3., 0., 0., 4., 
-   3., -1.}, {3., 0., 0., 3., 4., -1.}, {0., 0., 0., 1., 0., 0.}, {0.,
-    0., 0., 0., 1., 0.}, {2., 0., 0., 4., 4., -1.}}, {{1., 0., 0., 0.,
-    0., 0.}, {3., 0., 4., 0., 3., -1.}, {0., 0., 1., 0., 0., 0.}, {3.,
-    0., 3., 0., 4., -1.}, {0., 0., 0., 0., 1., 0.}, {2., 0., 4., 0., 
-   4., -1.}}, {{-1., 4., 4., 0., 0., 2.}, {0., 1., 0., 0., 0., 
-   0.}, {0., 0., 1., 0., 0., 0.}, {-1., 4., 3., 0., 0., 3.}, {-1., 3.,
-    4., 0., 0., 3.}, {0., 0., 0., 0., 0., 1.}}, {{0., 3., 0., 4., -1.,
-    3.}, {0., 1., 0., 0., 0., 0.}, {0., 3., 0., 3., -1., 4.}, {0., 0.,
-    0., 1., 0., 0.}, {0., 2., 0., 4., -1., 4.}, {0., 0., 0., 0., 0., 
-   1.}}, {{-1., 0., 4., 0., 4., 2.}, {-1., 0., 4., 0., 3., 3.}, {0., 
-   0., 1., 0., 0., 0.}, {-1., 0., 3., 0., 4., 3.}, {0., 0., 0., 0., 
-   1., 0.}, {0., 0., 0., 0., 0., 1.}}, {{0., 0., -1., 3., 4., 
-   3.}, {0., 0., -1., 3., 3., 4.}, {0., 0., -1., 2., 4., 4.}, {0., 0.,
-    0., 1., 0., 0.}, {0., 0., 0., 0., 1., 0.}, {0., 0., 0., 0., 0., 
-   1.}}}
-
-{{10., 8., -5.656854249492381, 7.}, {24., 
-  24., -16.970562748477143, 17.}, {8., 10., -5.656854249492381, 
-  7.}, {12., 10., -8.485281374238571, 7.}, {-4., -4., 
-  2.8284271247461903, -3.}, {10., 12., -8.485281374238571, 7.}}
 
 {{0, 1, 3}, {0, 2, 1}, {0, 3, 4}, {0, 4, 2}, {1, 2, 5}, {1, 5, 3}, {2, 4, 5}, {3, 5, 4}}

fractal_dimension/circle_counting_new/data/tetrahedron.txt

@@ -1,9 +1,8 @@
-{{{1., 0., 0., 0.}, {0., 1., 0., 0.}, {2., 2., -1., 2.}, {0., 0., 0., 
+{
-   1.}}, {{1., 0., 0., 0.}, {0., 1., 0., 0.}, {0., 0., 1., 0.}, {2., 
+    {1.0, -1.0, -1.0, -1.0},
-   2., 2., -1.}}, {{1., 0., 0., 0.}, {2., -1., 2., 2.}, {0., 0., 1., 
+    {-1.0, 1.0, -1.0, -1.0},
-   0.}, {0., 0., 0., 1.}}, {{-1., 2., 2., 2.}, {0., 1., 0., 0.}, {0., 
+    {-1.0, -1.0, 1.0, -1.0},
-   0., 1., 0.}, {0., 0., 0., 1.}}}
+    {-1.0, -1.0, -1.0, 1.0},
-
+}
-{{6., 4., -4., 3.}, {4., 4., -4., 1.}, {4., 6., -4., 3.}, {-2., -2., 2., -1.}}
 
 {{1, 2, 3}, {0, 2, 3}, {0, 1, 3}, {0, 1, 2}}

fractal_dimension/circle_counting_new/src/main.rs

@@ -1,6 +1,11 @@
 #![allow(dead_code)]
 
-use crate::{fractal::fractal_dimension, parser::read_file};
+use crate::{
+    fractal::fractal_dimension,
+    gram_matrix::{algebraic_generators, geometric_generators, root_tuple},
+    parser::read_file,
+};
+use nalgebra::DVector;
 use std::process;
 use structopt::StructOpt;
 
@@ -8,9 +13,9 @@ use ansi_term::Color::Yellow;
 
 pub mod constants;
 pub mod fractal;
+pub mod gram_matrix;
 pub mod parser;
 pub mod search;
-pub mod gram_matrix;
 
 /// Compute fractal dimension of crystallographic packings via the circle counting method
 #[derive(StructOpt)]
@@ -53,8 +58,7 @@ fn main() {
 
     let beginning = std::time::Instant::now();
 
-    let (mut generators, root, faces, orthogonal_generators) = read_file(&opt.data_file)
+    let (gram_matrix, faces) = read_file(&opt.data_file).unwrap_or_else(|err| {
-        .unwrap_or_else(|err| {
         eprintln!("{}", err);
         process::exit(-1);
     });
@@ -82,52 +86,109 @@ fn main() {
         );
     }
 
-    if !opt.geometric {
-        let c = nalgebra::DMatrix::from_columns(&root);
-        let ct = c
-            .clone()
-            .pseudo_inverse(1e-5)
-            .expect(format!("Invalid root matrix {}", c).as_str());
-        println!("{}", c);
-        println!("{}", ct);
-
-        generators = generators
-            .iter()
-            .map(|sigma| &c * sigma.transpose() * &ct)
-            .collect();
-    }
-
     if debug {
-        println!(
-            "{} (parsed from file {})",
-            Yellow.paint("Generators"),
-            opt.data_file
-        );
-        for generator in &generators {
-            println!("{}", generator);
-        }
         println!(
             "{} (parsed from file {}):",
-            Yellow.paint("Root Tuple"),
+            Yellow.paint("Gram Matrix"),
             opt.data_file
         );
-        for circle in &root {
+        println!("{}\n", gram_matrix);
-            println!("{}", circle);
-        }
         println!(
             "{} (parsed from file {}):",
             Yellow.paint("Faces"),
             opt.data_file
         );
         println!("{:?}\n", faces);
-        println!(
+        /* println!(
             "{} (parsed from file {}):",
             Yellow.paint("Orthogonal Generators"),
             opt.data_file
         );
-        println!("{:?}\n", orthogonal_generators);
+        println!("{:?}\n", orthogonal_generators); */
     }
 
+    // finding a suitable root tuple
+    let mut root = None;
+    for face in faces.iter().skip(1) {
+        for vertex in face {
+            let mut valid = true;
+            for vertex2 in &faces[0] {
+                if *vertex == *vertex2 {
+                    valid = false;
+                    break;
+                }
+            }
+            if valid {
+                let temp = root_tuple(
+                    &gram_matrix,
+                    (faces[0][0], faces[0][1], faces[0][2]),
+                    *vertex,
+                );
+
+                let mut valid = true;
+
+                let generators = geometric_generators(&gram_matrix, &faces, &temp);
+                for gen in generators {
+                    if gen[(1, 0)].abs() <= 1e-8 {
+                        valid = false;
+                        break;
+                    }
+                }
+
+                if valid {
+                    root = Some(temp);
+                    println!("{}\t{:?}", vertex, faces[0]);
+                    break;
+                }
+            }
+        }
+        if root.is_some() {
+            break;
+        }
+    }
+
+    if root.is_none() {
+        panic!("Invalid face scheme!");
+    }
+
+    let root = root.unwrap();
+
+    let generators = geometric_generators(&gram_matrix, &faces, &root);
+
+    if debug {
+        println!(
+            "{} (computed using gram_matrix):",
+            Yellow.paint("Root Matrix")
+        );
+        println!("{}", root);
+
+        println!(
+            "{} (computed using gram_matrix and root):",
+            Yellow.paint("Geometric generators")
+        );
+        for generator in &generators {
+            println!("{}", generator);
+        }
+
+        let generators = algebraic_generators(&gram_matrix, &faces);
+        println!(
+            "{} (computed using gram_matrix and root):",
+            Yellow.paint("Algebraic generators")
+        );
+        for generator in &generators {
+            println!("{}", generator);
+        }
+    }
+
+    let mut temp = vec![];
+    for circle in root.column_iter() {
+        temp.push(DVector::from_iterator(
+            4,
+            circle.iter().map(|val| val.clone()),
+        ));
+    }
+    let root = temp;
+
     let delta = fractal_dimension(
         generators,
         root,
@@ -136,7 +197,7 @@ fn main() {
         debug,
         depth,
         faces,
-        orthogonal_generators,
+        vec![],
     )
     .unwrap();
     let after_computing = std::time::Instant::now();

fractal_dimension/circle_counting_new/src/parser.rs

@@ -1,6 +1,6 @@
 use std::fs;
 
-use nalgebra::{DMatrix, DVector};
+use nalgebra::DMatrix;
 use ansi_term::Color::Red;
 
 #[derive(Debug, PartialEq)]
@@ -213,11 +213,9 @@ fn matrix_to_rust_value_flat(value: &Value) -> Result<Vec<f64>, String> {
     Ok(result)
 }
 
-pub type Generator = DMatrix<f64>;
+pub type Matrix = DMatrix<f64>;
-pub type Root = Vec<DVector<f64>>;
 pub type FaceList = Vec<Vec<usize>>;
-pub type OrthogonalGenerators = Vec<Vec<usize>>;
+pub type Data = (Matrix, FaceList);
-pub type Data = (Vec<Generator>, Root, FaceList, OrthogonalGenerators);
 pub fn read_file(filename: &str) -> Result<Data, String> {
     let contents = match fs::read_to_string(filename) {
         Ok(contents) => contents,
@@ -226,28 +224,16 @@ pub fn read_file(filename: &str) -> Result<Data, String> {
 
     let results = eval_full(contents)?;
 
-    let mut generators = vec![];
+    let n = matrix_to_rust_value(&results[0])?.len();
-    if let Value::List(gens) = &results[0] {
+    let gram_matrix = DMatrix::from_row_slice(n, n, &matrix_to_rust_value_flat(&results[0])?);
-        for val in gens {
-            let mat = matrix_to_rust_value(val)?;
-            let mat2 = matrix_to_rust_value_flat(val)?;
-            generators.push(DMatrix::from_row_slice(mat.len(), mat[0].len(), &mat2));
-        }
-    } else {
-        return Err(format!(
-            "Expected first value in file to be the generators; got {:?}",
-            results[0]
-        ));
-    }
 
-    let root = matrix_to_rust_value(&results[1])?;
+    let faces = matrix_to_rust_value(&results[1])?;
-    let faces = matrix_to_rust_value(&results[2])?;
     let faces = faces
         .iter()
         .map(|v| v.iter().map(|x| x.round() as usize).collect())
         .collect();
 
-    let orthogonal_generators: Vec<Vec<usize>> = if results.len() < 4 {
+    let orthogonal_generators: Vec<Vec<usize>> = if results.len() < 3 {
         vec![]
     } else {
         matrix_to_rust_value(&results[3])?
@@ -263,8 +249,7 @@ pub fn read_file(filename: &str) -> Result<Data, String> {
         }
     }
 
-    let new_root = root.iter().map(|circ| DVector::from_column_slice(circ)).collect();
+    Ok((gram_matrix, faces))
-    Ok((generators, new_root, faces, orthogonal_generators))
 }
 
 #[cfg(test)]