Factor out the kaleidocycle realization

This parallels what we did for the Irisawa hexlet realization. The
kaleidocycle tangent test comes out slightly weaker, because we no
longer confirm that the realized configuration matches the initial
guess. However, we still confirm that the configuration history only has
one entry, which is equivalent as long as the configuration history
starts with the initial guess and is updated after every optimization
step.

app-proto/examples/irisawa-hexlet.rs

@@ -1,4 +1,4 @@
-use dyna3::engine::{Q, irisawa::realize_irisawa_hexlet};
+use dyna3::engine::{Q, examples::realize_irisawa_hexlet};
 
 fn main() {
     const SCALED_TOL: f64 = 1.0e-12;

app-proto/examples/kaleidocycle.rs

@@ -1,53 +1,10 @@
 use nalgebra::{DMatrix, DVector};
-use std::{array, f64::consts::PI};
 
-use dyna3::engine::{Q, point, realize_gram, PartialMatrix};
+use dyna3::engine::{Q, examples::realize_kaleidocycle};
 
 fn main() {
-    // set up a kaleidocycle, made of points with fixed distances between them,
+    const SCALED_TOL: f64 = 1.0e-12;
-    // and find its tangent space
+    let (config, tangent, success, history) = realize_kaleidocycle(SCALED_TOL);
-    const N_POINTS: usize = 12;
-    let gram = {
-        let mut gram_to_be = PartialMatrix::new();
-        for block in (0..N_POINTS).step_by(2) {
-            let block_next = (block + 2) % N_POINTS;
-            for j in 0..2 {
-                // diagonal and hinge edges
-                for k in j..2 {
-                    gram_to_be.push_sym(block + j, block + k, if j == k { 0.0 } else { -0.5 });
-                }
-                
-                // non-hinge edges
-                for k in 0..2 {
-                    gram_to_be.push_sym(block + j, block_next + k, -0.625);
-                }
-            }
-        }
-        gram_to_be
-    };
-    let guess = {
-        const N_HINGES: usize = 6;
-        let guess_elts = (0..N_HINGES).step_by(2).flat_map(
-            |n| {
-                let ang_hor = (n as f64) * PI/3.0;
-                let ang_vert = ((n + 1) as f64) * PI/3.0;
-                let x_vert = ang_vert.cos();
-                let y_vert = ang_vert.sin();
-                [
-                    point(0.0, 0.0, 0.0),
-                    point(ang_hor.cos(), ang_hor.sin(), 0.0),
-                    point(x_vert, y_vert, -0.5),
-                    point(x_vert, y_vert, 0.5)
-                ]
-            }
-        ).collect::<Vec<_>>();
-        DMatrix::from_columns(&guess_elts)
-    };
-    let frozen: [_; N_POINTS] = array::from_fn(|k| (3, k));
-    let (config, tangent, success, history) = realize_gram(
-        &gram, guess, &frozen,
-        1.0e-12, 0.5, 0.9, 1.1, 200, 110
-    );
     print!("Completed Gram matrix:{}", config.tr_mul(&*Q) * &config);
     print!("Configuration:{}", config);
     if success {
@@ -58,7 +15,8 @@ fn main() {
     println!("Steps: {}", history.scaled_loss.len() - 1);
     println!("Loss: {}\n", history.scaled_loss.last().unwrap());
     
-    // find the kaleidocycle's twist motion
+    // find the kaleidocycle's twist motion by projecting onto the tangent space
+    const N_POINTS: usize = 12;
     let up = DVector::from_column_slice(&[0.0, 0.0, 1.0, 0.0]);
     let down = -&up;
     let twist_motion: DMatrix<_> = (0..N_POINTS).step_by(4).flat_map(

app-proto/src/engine.rs

@@ -413,20 +413,20 @@ pub fn realize_gram(
 
 // --- tests ---
 
-// this problem is from a sangaku by Irisawa Shintarō Hiroatsu. the article
-// below includes a nice translation of the problem statement, which was
-// recorded in Uchida Itsumi's book _Kokon sankan_ (_Mathematics, Past and
-// Present_)
-//
-//   "Japan's 'Wasan' Mathematical Tradition", by Abe Haruki
-//   https://www.nippon.com/en/japan-topics/c12801/
-//
 #[cfg(feature = "dev")]
-pub mod irisawa {
+pub mod examples {
     use std::{array, f64::consts::PI};
     
     use super::*;
     
+    // this problem is from a sangaku by Irisawa Shintarō Hiroatsu. the article
+    // below includes a nice translation of the problem statement, which was
+    // recorded in Uchida Itsumi's book _Kokon sankan_ (_Mathematics, Past and
+    // Present_)
+    //
+    //   "Japan's 'Wasan' Mathematical Tradition", by Abe Haruki
+    //   https://www.nippon.com/en/japan-topics/c12801/
+    //
     pub fn realize_irisawa_hexlet(scaled_tol: f64) -> (DMatrix<f64>, ConfigSubspace, bool, DescentHistory) {
         let gram = {
             let mut gram_to_be = PartialMatrix::new();
@@ -480,14 +480,64 @@ pub mod irisawa {
             scaled_tol, 0.5, 0.9, 1.1, 200, 110
         )
     }
+    
+    // set up a kaleidocycle, made of points with fixed distances between them,
+    // and find its tangent space
+    pub fn realize_kaleidocycle(scaled_tol: f64) -> (DMatrix<f64>, ConfigSubspace, bool, DescentHistory) {
+        const N_POINTS: usize = 12;
+        let gram = {
+            let mut gram_to_be = PartialMatrix::new();
+            for block in (0..N_POINTS).step_by(2) {
+                let block_next = (block + 2) % N_POINTS;
+                for j in 0..2 {
+                    // diagonal and hinge edges
+                    for k in j..2 {
+                        gram_to_be.push_sym(block + j, block + k, if j == k { 0.0 } else { -0.5 });
+                    }
+                    
+                    // non-hinge edges
+                    for k in 0..2 {
+                        gram_to_be.push_sym(block + j, block_next + k, -0.625);
+                    }
+                }
+            }
+            gram_to_be
+        };
+        
+        let guess = {
+            const N_HINGES: usize = 6;
+            let guess_elts = (0..N_HINGES).step_by(2).flat_map(
+                |n| {
+                    let ang_hor = (n as f64) * PI/3.0;
+                    let ang_vert = ((n + 1) as f64) * PI/3.0;
+                    let x_vert = ang_vert.cos();
+                    let y_vert = ang_vert.sin();
+                    [
+                        point(0.0, 0.0, 0.0),
+                        point(ang_hor.cos(), ang_hor.sin(), 0.0),
+                        point(x_vert, y_vert, -0.5),
+                        point(x_vert, y_vert, 0.5)
+                    ]
+                }
+            ).collect::<Vec<_>>();
+            DMatrix::from_columns(&guess_elts)
+        };
+        
+        let frozen: [_; N_POINTS] = array::from_fn(|k| (3, k));
+        
+        realize_gram(
+            &gram, guess, &frozen,
+            scaled_tol, 0.5, 0.9, 1.1, 200, 110
+        )
+    }
 }
 
 #[cfg(test)]
 mod tests {
     use nalgebra::Vector3;
-    use std::{array, f64::consts::{FRAC_1_SQRT_2, PI}, iter};
+    use std::{f64::consts::{FRAC_1_SQRT_2, PI}, iter};
     
-    use super::{*, irisawa::realize_irisawa_hexlet};
+    use super::{*, examples::*};
     
     #[test]
     fn sub_proj_test() {
@@ -523,7 +573,7 @@ mod tests {
             entries
         });
         let config = {
-            let a: f64 = 0.75_f64.sqrt();
+            let a = 0.75_f64.sqrt();
             DMatrix::from_columns(&[
                 sphere(1.0, 0.0, 0.0, a),
                 sphere(-0.5, a, 0.0, a),
@@ -534,6 +584,40 @@ mod tests {
         assert!(state.loss.abs() < f64::EPSILON);
     }
     
+    // at the frozen indices, the optimization steps should have exact zeros,
+    // and the realized configuration should match the initial guess
+    #[test]
+    fn frozen_entry_test() {
+        let gram = {
+            let mut gram_to_be = PartialMatrix::new();
+            for j in 0..2 {
+                for k in j..2 {
+                    gram_to_be.push_sym(j, k, if (j, k) == (1, 1) { 1.0 } else { 0.0 });
+                }
+            }
+            gram_to_be
+        };
+        let guess = DMatrix::from_columns(&[
+            point(0.0, 0.0, 2.0),
+            sphere(0.0, 0.0, 0.0, 1.0)
+        ]);
+        let frozen = [(3, 0), (3, 1)];
+        println!();
+        let (config, _, success, history) = realize_gram(
+            &gram, guess.clone(), &frozen,
+            1.0e-12, 0.5, 0.9, 1.1, 200, 110
+        );
+        assert_eq!(success, true);
+        for base_step in history.base_step.into_iter() {
+            for index in frozen {
+                assert_eq!(base_step[index], 0.0);
+            }
+        }
+        for index in frozen {
+            assert_eq!(config[index], guess[index]);
+        }
+    }
+    
     #[test]
     fn irisawa_hexlet_test() {
         // solve Irisawa's problem
@@ -574,12 +658,8 @@ mod tests {
         assert_eq!(success, true);
         assert_eq!(history.scaled_loss.len(), 1);
         
-        // confirm that the tangent space has dimension five or less
+        // list some motions that should form a basis for the tangent space of
-        assert_eq!(tangent.basis_std.len(), 5);
+        // the solution variety
-        
-        // confirm that the tangent space contains all the motions we expect it
-        // to. since we've already bounded the dimension of the tangent space,
-        // this confirms that the tangent space is what we expect it to be
         const UNIFORM_DIM: usize = 4;
         let element_dim = guess.nrows();
         let assembly_dim = guess.ncols();
@@ -605,6 +685,14 @@ mod tests {
                 0.0, 0.0, -1.0,  0.25,  1.0
             ])
         ];
+        
+        // confirm that the dimension of the tangent space is no greater than
+        // expected
+        assert_eq!(tangent.basis_std.len(), tangent_motions_std.len());
+        
+        // confirm that the tangent space contains all the motions we expect it
+        // to. since we've already bounded the dimension of the tangent space,
+        // this confirms that the tangent space is what we expect it to be
         let tol_sq = ((element_dim * assembly_dim) as f64) * SCALED_TOL * SCALED_TOL;
         for (motion_unif, motion_std) in tangent_motions_unif.into_iter().zip(tangent_motions_std) {
             let motion_proj: DMatrix<_> = motion_unif.column_iter().enumerate().map(
@@ -633,59 +721,17 @@ mod tests {
     
     #[test]
     fn tangent_test_kaleidocycle() {
-        // set up a kaleidocycle, made of points with fixed distances between
+        // set up a kaleidocycle and find its tangent space
-        // them, and find its tangent space
-        const N_POINTS: usize = 12;
-        const N_HINGES: usize = 6;
         const SCALED_TOL: f64 = 1.0e-12;
-        let gram = {
+        let (config, tangent, success, history) = realize_kaleidocycle(SCALED_TOL);
-            let mut gram_to_be = PartialMatrix::new();
-            for block in (0..N_POINTS).step_by(2) {
-                let block_next = (block + 2) % N_POINTS;
-                for j in 0..2 {
-                    // diagonal and hinge edges
-                    for k in j..2 {
-                        gram_to_be.push_sym(block + j, block + k, if j == k { 0.0 } else { -0.5 });
-                    }
-                    
-                    // non-hinge edges
-                    for k in 0..2 {
-                        gram_to_be.push_sym(block + j, block_next + k, -0.625);
-                    }
-                }
-            }
-            gram_to_be
-        };
-        let guess = {
-            let guess_elts = (0..N_HINGES).step_by(2).flat_map(
-                |n| {
-                    let ang_hor = (n as f64) * PI/3.0;
-                    let ang_vert = ((n + 1) as f64) * PI/3.0;
-                    let x_vert = ang_vert.cos();
-                    let y_vert = ang_vert.sin();
-                    [
-                        point(0.0, 0.0, 0.0),
-                        point(ang_hor.cos(), ang_hor.sin(), 0.0),
-                        point(x_vert, y_vert, -0.5),
-                        point(x_vert, y_vert, 0.5)
-                    ]
-                }
-            ).collect::<Vec<_>>();
-            DMatrix::from_columns(&guess_elts)
-        };
-        let frozen: [_; N_POINTS] = array::from_fn(|k| (3, k));
-        let (config, tangent, success, history) = realize_gram(
-            &gram, guess.clone(), &frozen,
-            SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
-        );
-        assert_eq!(config, guess);
         assert_eq!(success, true);
         assert_eq!(history.scaled_loss.len(), 1);
         
         // list some motions that should form a basis for the tangent space of
         // the solution variety
-        let element_dim = guess.nrows();
+        const N_HINGES: usize = 6;
-        let assembly_dim = guess.ncols();
+        let element_dim = config.nrows();
+        let assembly_dim = config.ncols();
         let tangent_motions_unif = vec![
             // the translations along the coordinate axes
             translation_motion_unif(&Vector3::new(1.0, 0.0, 0.0), assembly_dim),
@@ -693,9 +739,9 @@ mod tests {
             translation_motion_unif(&Vector3::new(0.0, 0.0, 1.0), assembly_dim),
             
             // the rotations about the coordinate axes
-            rotation_motion_unif(&Vector3::new(1.0, 0.0, 0.0), guess.column_iter().collect()),
+            rotation_motion_unif(&Vector3::new(1.0, 0.0, 0.0), config.column_iter().collect()),
-            rotation_motion_unif(&Vector3::new(0.0, 1.0, 0.0), guess.column_iter().collect()),
+            rotation_motion_unif(&Vector3::new(0.0, 1.0, 0.0), config.column_iter().collect()),
-            rotation_motion_unif(&Vector3::new(0.0, 0.0, 1.0), guess.column_iter().collect()),
+            rotation_motion_unif(&Vector3::new(0.0, 0.0, 1.0), config.column_iter().collect()),
             
             // the twist motion. more precisely: a motion that keeps the center
             // of mass stationary and preserves the distances between the
@@ -720,7 +766,7 @@ mod tests {
         ];
         let tangent_motions_std = tangent_motions_unif.iter().map(
             |motion| DMatrix::from_columns(
-                &guess.column_iter().zip(motion).map(
+                &config.column_iter().zip(motion).map(
                     |(v, elt_motion)| local_unif_to_std(v) * elt_motion
                 ).collect::<Vec<_>>()
             )
@@ -826,38 +872,4 @@ mod tests {
         let tol_sq = ((guess_orig.nrows() * guess_orig.ncols()) as f64) * SCALED_TOL_TFM * SCALED_TOL_TFM;
         assert!((motion_proj_tfm - motion_tfm_proj).norm_squared() < tol_sq);
     }
-    
-    // at the frozen indices, the optimization steps should have exact zeros,
-    // and the realized configuration should match the initial guess
-    #[test]
-    fn frozen_entry_test() {
-        let gram = {
-            let mut gram_to_be = PartialMatrix::new();
-            for j in 0..2 {
-                for k in j..2 {
-                    gram_to_be.push_sym(j, k, if (j, k) == (1, 1) { 1.0 } else { 0.0 });
-                }
-            }
-            gram_to_be
-        };
-        let guess = DMatrix::from_columns(&[
-            point(0.0, 0.0, 2.0),
-            sphere(0.0, 0.0, 0.0, 1.0)
-        ]);
-        let frozen = [(3, 0), (3, 1)];
-        println!();
-        let (config, _, success, history) = realize_gram(
-            &gram, guess.clone(), &frozen,
-            1.0e-12, 0.5, 0.9, 1.1, 200, 110
-        );
-        assert_eq!(success, true);
-        for base_step in history.base_step.into_iter() {
-            for index in frozen {
-                assert_eq!(base_step[index], 0.0);
-            }
-        }
-        for index in frozen {
-            assert_eq!(config[index], guess[index]);
-        }
-    }
 }