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,
-    // and find its tangent space
-    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
-    );
+    const SCALED_TOL: f64 = 1.0e-12;
+    let (config, tangent, success, history) = realize_kaleidocycle(SCALED_TOL);
     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
-        assert_eq!(tangent.basis_std.len(), 5);
-        
-        // 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
+        // list some motions that should form a basis for the tangent space of
+        // the solution variety
         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
-        // them, and find its tangent space
-        const N_POINTS: usize = 12;
-        const N_HINGES: usize = 6;
+        // set up a kaleidocycle and find its tangent space
         const SCALED_TOL: f64 = 1.0e-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 = {
-            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);
+        let (config, tangent, success, history) = realize_kaleidocycle(SCALED_TOL);
         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();
-        let assembly_dim = guess.ncols();
+        const N_HINGES: usize = 6;
+        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(0.0, 1.0, 0.0), guess.column_iter().collect()),
-            rotation_motion_unif(&Vector3::new(0.0, 0.0, 1.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), config.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]);
-        }
-    }
 }