working?

fractal_dimension/circle_counting_new/data/tetrahedron.txt

@@ -1,30 +1,15 @@
-{
+{{{1.33333, 1.16068, -2.1547, -1.24402}, {0.095729, 
-    {
+   1.33333, -0.618802, -0.357266}, {0.309401, 
-        {-1, 2, 2, 2},
+   1.07735, -1., -1.1547}, {0.178633, 0.622008, -1.1547, 
-        {0, 1, 0, 0},
+   0.333333}}, {{1.33333, 1.16068, 2.1547, -1.24402}, {0.095729, 
-        {0, 0, 1, 0},
+   1.33333, 0.618802, -0.357266}, {-0.309401, -1.07735, -1., 
-        {0, 0, 0, 1}
+   1.1547}, {0.178633, 0.622008, 1.1547, 0.333333}}, {{0., 0.25, 0., 
-    },
+   0.}, {4., 0., 0., 0.}, {0., 0., 1., 0.}, {0., 0., 0., 
-    {
+   1.}}, {{1.33333, 1.16068, 0., 2.48803}, {0.095729, 1.33333, 0., 
-        {1, 0, 0, 0},
+   0.714531}, {0., 0., 1., 0.}, {-0.357266, -1.24402, 0., -1.66667}}}
-        {2, -1, 2, 2},
-        {0, 0, 1, 0},
-        {0, 0, 0, 1}
-    },
-    {
-        {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, 2, 2, -1}
-    }
-}
 
-{-6, 11, 14, 15}
+{{-1.86603, 0.535898, 0., 0.}, {0.288675, 1.1547, 0., 
+  1.1547}, {0.288675, 1.1547, 1., -0.57735}, {0.288675, 
+  1.1547, -1., -0.57735}}
 
 {{1, 2, 3}, {0, 2, 3}, {0, 1, 3}, {0, 1, 2}}

fractal_dimension/circle_counting_new/src/constants.rs

@@ -0,0 +1,2 @@
+pub type I = u64;
+pub type F = f64;

fractal_dimension/circle_counting_new/src/fractal.rs

@@ -0,0 +1,57 @@
+use ansi_term::Color::Yellow;
+use linregress::{FormulaRegressionBuilder, RegressionDataBuilder};
+use nalgebra::{DMatrix, DVector};
+
+use crate::search::Searcher;
+
+pub fn fractal_dimension(
+    generators: Vec<DMatrix<f64>>,
+    root: Vec<DVector<f64>>,
+    upper_bound: f64,
+    n: usize,
+    debug: bool,
+    generations: usize,
+    orthogonal_generators: Vec<Vec<usize>>,
+) -> Result<f64, linregress::Error> {
+    let xs: Vec<f64> = (1..=n)
+        .map(|x| (x as f64 * upper_bound / (2.0 * n as f64) + upper_bound / 2.0))
+        .collect();
+
+    let mut searcher = Searcher::new(&xs, &generators);
+    for circle in root {
+        searcher.search(&circle, std::usize::MAX, 0, generations);
+    }
+
+    if debug {
+        println!("{}", Yellow.paint("Done counting circles!"));
+        println!("sample points:\n{:?}", xs);
+        println!(
+            "\nnumber of circles fewer than each of those sample points:\n{:?}",
+            searcher.counts
+        );
+    }
+
+    let xs: Vec<f64> = xs.iter().map(|x| x.ln()).collect();
+    let ys: Vec<f64> = searcher.counts.iter().map(|x| (*x as f64).ln()).collect();
+    let data = vec![("Y", ys), ("X", xs)];
+    let data = RegressionDataBuilder::new().build_from(data)?;
+
+    let formula = "Y ~ X";
+    let model = FormulaRegressionBuilder::new()
+        .data(&data)
+        .formula(formula)
+        .fit()?;
+
+    if debug {
+        println!("\n{}", Yellow.paint("Built regression model!"));
+        println!("Model info:");
+        println!("\tr^2:\t\t{}", model.rsquared);
+        println!("\tp-value:\t{}", model.pvalues.pairs()[0].1);
+        println!("\tintercept:\t{}", model.parameters.intercept_value);
+        println!("\tslope:\t{}", model.parameters.regressor_values[0]);
+    }
+
+    let parameters = model.parameters;
+
+    Ok(parameters.regressor_values[0])
+}

fractal_dimension/circle_counting_new/src/lib.rs

@@ -1,131 +0,0 @@
-use ansi_term::Color::Yellow;
-use linregress::{FormulaRegressionBuilder, RegressionDataBuilder};
-use nalgebra::{DMatrix, DVector};
-
-pub fn fractal_dimension(
-    generators: Vec<DMatrix<f64>>,
-    root: DVector<f64>,
-    faces: Vec<Vec<usize>>,
-    upper_bound: f64,
-    n: usize,
-    debug: bool,
-    generations: usize,
-    orthogonal_generators: Vec<Vec<usize>>,
-) -> Result<f64, linregress::Error> {
-    let mut totals = vec![root.len(); n];
-    let mut current = vec![(root, std::usize::MAX, false)];
-    let mut next = vec![];
-    let mut nodes: u64 = 1;
-
-    let xs: Vec<f64> = (1..=n)
-        .map(|x| (x as f64 * upper_bound / (2.0 * n as f64) + upper_bound / 2.0))
-        .collect();
-
-    let mut i = 0;
-
-    loop {
-        next.clear();
-        for (tuple, previous_generator, tuple_too_big) in &current {
-            let mut add_children = false;
-            let mut children = vec![];
-            for (i, generator) in generators.iter().enumerate() {
-                let mut skip = false;
-                for orthogonal_pairs in &orthogonal_generators {
-                    if i == orthogonal_pairs[1] && *previous_generator == orthogonal_pairs[0] {
-                        skip = true;
-                        break;
-                    }
-                }
-                if skip {
-                    continue;
-                }
-                if i != *previous_generator {
-                    let new_tuple = generator * tuple;
-                    if new_tuple.iter().sum::<f64>() < tuple.iter().sum() {
-                        continue;
-                    }
-                    let mut too_big = true;
-                    for (j, curvature) in new_tuple.iter().enumerate() {
-                        let mut skip = false;
-                        for vertex in &faces[i] {
-                            if j == *vertex {
-                                skip = true;
-                                break;
-                            }
-                        }
-                        if skip {
-                            continue;
-                        }
-                        for (k, max) in xs.iter().enumerate() {
-                            if curvature <= max {
-                                add_children = true;
-                                too_big = false;
-                                totals[k] += 1;
-                            }
-                        }
-                    }
-                    children.push((new_tuple, i, too_big));
-                }
-            }
-            if add_children || !tuple_too_big {
-                nodes += 1;
-                for child in children {
-                    next.push(child);
-                }
-            }
-        }
-        std::mem::swap(&mut current, &mut next);
-
-        i += 1;
-        if generations != 0 && i > generations {
-            break;
-        }
-        if debug {
-            println!("Generation {}:", i);
-            println!("\tnumber of leaves:\t{}", current.len());
-            if !current.is_empty() {
-                let tuple = &current[current.len() / 2];
-                println!("\trandom tuple:\t\t{}", tuple.0);
-            }
-            println!();
-        }
-
-        if current.is_empty() {
-            break;
-        }
-    }
-
-    if debug {
-        println!("{}", Yellow.paint("Done counting circles!"));
-        println!("sample points:\n{:?}", xs);
-        println!(
-            "\nnumber of circles fewer than each of those sample points:\n{:?}",
-            totals
-            );
-        println!("\nTotal number of nodes:\t{}", nodes);
-    }
-
-    let xs: Vec<f64> = xs.iter().map(|x| x.ln()).collect();
-    let ys: Vec<f64> = totals.iter().map(|x| (*x as f64).ln()).collect();
-    let data = vec![("Y", ys), ("X", xs)];
-    let data = RegressionDataBuilder::new().build_from(data)?;
-
-    let formula = "Y ~ X";
-    let model = FormulaRegressionBuilder::new()
-        .data(&data)
-        .formula(formula)
-        .fit()?;
-
-    if debug {
-        println!("\n{}", Yellow.paint("Built regression model!"));
-        println!("Model info:");
-        println!("\tr^2:\t\t{}", model.rsquared);
-        println!("\tp-value:\t{}", model.pvalues.pairs()[0].1);
-        println!("\tintercept:\t{}", model.parameters.intercept_value);
-        println!("\tslope:\t{}", model.parameters.regressor_values[0]);
-    }
-
-    let parameters = model.parameters;
-
-    Ok(parameters.regressor_values[0])
-}

fractal_dimension/circle_counting_new/src/main.rs

@@ -1,13 +1,15 @@
 #![allow(dead_code)]
 
-use crate::{lib::fractal_dimension, parser::read_file};
+use crate::{fractal::fractal_dimension, parser::read_file};
 use std::process;
 use structopt::StructOpt;
 
 use ansi_term::Color::Yellow;
 
-mod lib;
+pub mod fractal;
-mod parser;
+pub mod parser;
+pub mod constants;
+pub mod search;
 
 /// Compute fractal dimension of crystallographic packings via the circle counting method
 #[derive(StructOpt)]
@@ -89,7 +91,7 @@ fn main() {
             Yellow.paint("Root Tuple"),
             opt.data_file
         );
-        println!("{}", root);
+        println!("{:?}", root);
         println!(
             "{} (parsed from file {}):",
             Yellow.paint("Faces"),
@@ -107,7 +109,6 @@ fn main() {
     let delta = fractal_dimension(
         generators,
         root,
-        faces,
         opt.max,
         opt.n,
         debug,

fractal_dimension/circle_counting_new/src/parser.rs

@@ -214,7 +214,7 @@ fn matrix_to_rust_value_flat(value: &Value) -> Result<Vec<f64>, String> {
 }
 
 pub type Generator = DMatrix<f64>;
-pub type Root = DVector<f64>;
+pub type Root = Vec<DVector<f64>>;
 pub type FaceList = Vec<Vec<usize>>;
 pub type OrthogonalGenerators = Vec<Vec<usize>>;
 pub type Data = (Vec<Generator>, Root, FaceList, OrthogonalGenerators);
@@ -240,7 +240,7 @@ pub fn read_file(filename: &str) -> Result<Data, String> {
         ));
     }
 
-    let root = vector_to_rust_value(&results[1])?;
+    let root = matrix_to_rust_value(&results[1])?;
     let faces = matrix_to_rust_value(&results[2])?;
     let faces = faces
         .iter()
@@ -263,7 +263,8 @@ pub fn read_file(filename: &str) -> Result<Data, String> {
         }
     }
 
-    Ok((generators, DVector::from_column_slice(&root), faces, orthogonal_generators))
+    let new_root = root.iter().map(|circ| DVector::from_column_slice(circ)).collect();
+    Ok((generators, new_root, faces, orthogonal_generators))
 }
 
 #[cfg(test)]

fractal_dimension/circle_counting_new/src/search.rs

@@ -0,0 +1,47 @@
+use nalgebra::{DMatrix, DVector};
+
+pub struct Searcher<'a> {
+    pub counts: Vec<u64>,
+    maxes: &'a Vec<f64>,
+    generators: &'a Vec<DMatrix<f64>>,
+}
+
+impl<'a> Searcher<'a> {
+    pub fn new(maxes: &'a Vec<f64>, generators: &'a Vec<DMatrix<f64>>) -> Self {
+        Self {
+            counts: vec![0; maxes.len()],
+            maxes,
+            generators,
+        }
+    }
+
+    pub fn search(
+        &mut self,
+        circle: &DVector<f64>,
+        previous_generator: usize,
+        generations: usize,
+        max_generations: usize,
+    ) {
+        if generations > max_generations {
+            return;
+        }
+        for (i, generator) in self.generators.iter().enumerate() {
+            if i != previous_generator {
+                let new_circle = generator * circle;
+                {
+                    let mut seen = false;
+                    for (j, max) in self.maxes.iter().enumerate() {
+                        if seen || new_circle[1] > circle[1] && circle[1] < *max {
+                            seen = true;
+                            self.counts[j] += 1;
+                        }
+                    }
+                    if !seen {
+                        continue;
+                    }
+                } // ensure seen is dropped before recursive call to minimize stack frame size
+                self.search(&new_circle, i, generations + 1, max_generations);
+            }
+        }
+    }
+}