Test the kaleidocycle's first-order motions

The motions we feed into the projection map now need to be expressed in
uniform coordinates. I've verified by hand that `tangent_motions_unif`
and `tangent_motions_std` represent the same motions.

app-proto/src/engine.rs

@@ -490,6 +490,9 @@ pub mod irisawa {
 
 #[cfg(test)]
 mod tests {
+    use nalgebra::Vector3;
+    use std::{array, f64::consts::{FRAC_1_SQRT_2, PI}, iter};
+    
     use super::{*, irisawa::realize_irisawa_hexlet};
     
     #[test]
@@ -552,10 +555,8 @@ mod tests {
     }
     
     #[test]
-    fn tangent_test() {
+    fn tangent_test_three_spheres() {
         const SCALED_TOL: f64 = 1.0e-12;
-        const ELEMENT_DIM: usize = 5;
-        const ASSEMBLY_DIM: usize = 3;
         let gram = {
             let mut gram_to_be = PartialMatrix::new();
             for j in 0..3 {
@@ -580,32 +581,258 @@ mod tests {
         assert_eq!(history.scaled_loss.len(), 1);
         
         // confirm that the tangent space has dimension five or less
-        let ConfigSubspace(ref tangent_basis) = tangent;
+        assert_eq!(tangent.basis_std.len(), 5);
-        assert_eq!(tangent_basis.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
-        let tangent_motions = vec![
+        const UNIFORM_DIM: usize = 4;
-            basis_matrix((0, 1), ELEMENT_DIM, ASSEMBLY_DIM),
+        let element_dim = guess.nrows();
-            basis_matrix((1, 1), ELEMENT_DIM, ASSEMBLY_DIM),
+        let assembly_dim = guess.ncols();
-            basis_matrix((0, 2), ELEMENT_DIM, ASSEMBLY_DIM),
+        let tangent_motions_unif = vec![
-            basis_matrix((1, 2), ELEMENT_DIM, ASSEMBLY_DIM),
+            basis_matrix((0, 1), UNIFORM_DIM, assembly_dim),
-            DMatrix::<f64>::from_column_slice(ELEMENT_DIM, 3, &[
+            basis_matrix((1, 1), UNIFORM_DIM, assembly_dim),
-                0.0, 0.0, 0.0, 0.0, 0.0,
+            basis_matrix((0, 2), UNIFORM_DIM, assembly_dim),
-                0.0, 0.0, -1.0, -0.25, -1.0,
+            basis_matrix((1, 2), UNIFORM_DIM, assembly_dim),
-                0.0, 0.0, -1.0, 0.25, 1.0
+            DMatrix::<f64>::from_column_slice(UNIFORM_DIM, assembly_dim, &[
+                0.0, 0.0,  0.0,  0.0,
+                0.0, 0.0, -0.5, -0.5,
+                0.0, 0.0, -0.5,  0.5
             ])
         ];
-        let tol_sq = ((ELEMENT_DIM * ASSEMBLY_DIM) as f64) * SCALED_TOL * SCALED_TOL;
+        let tangent_motions_std = vec![
-        for motion in tangent_motions {
+            basis_matrix((0, 1), element_dim, assembly_dim),
-            let motion_proj: DMatrix<_> = motion.column_iter().enumerate().map(
+            basis_matrix((1, 1), element_dim, assembly_dim),
+            basis_matrix((0, 2), element_dim, assembly_dim),
+            basis_matrix((1, 2), element_dim, assembly_dim),
+            DMatrix::<f64>::from_column_slice(element_dim, assembly_dim, &[
+                0.0, 0.0,  0.0,  0.00,  0.0,
+                0.0, 0.0, -1.0, -0.25, -1.0,
+                0.0, 0.0, -1.0,  0.25,  1.0
+            ])
+        ];
+        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(
                 |(k, v)| tangent.proj(&v, k)
             ).sum();
-            assert!((motion - motion_proj).norm_squared() < tol_sq);
+            assert!((motion_std - motion_proj).norm_squared() < tol_sq);
         }
     }
     
+    fn translation_motion_unif(vel: &Vector3<f64>, assembly_dim: usize) -> Vec<DVector<f64>> {
+        let mut elt_motion = DVector::zeros(4);
+        elt_motion.fixed_rows_mut::<3>(0).copy_from(vel);
+        iter::repeat(elt_motion).take(assembly_dim).collect()
+    }
+    
+    fn rotation_motion_unif(ang_vel: &Vector3<f64>, points: Vec<DVectorView<f64>>) -> Vec<DVector<f64>> {
+        points.into_iter().map(
+            |pt| {
+                let vel = ang_vel.cross(&pt.fixed_rows::<3>(0));
+                let mut elt_motion = DVector::zeros(4);
+                elt_motion.fixed_rows_mut::<3>(0).copy_from(&vel);
+                elt_motion
+            }
+        ).collect()
+    }
+    
+    #[test]
+    fn tangent_test_kaleidocycle() {
+        const SCALED_TOL: f64 = 1.0e-12;
+        // 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;
+        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);
+        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();
+        let tangent_motions_unif = vec![
+            // the translations along the coordinate axes
+            translation_motion_unif(&Vector3::new(1.0, 0.0, 0.0), assembly_dim),
+            translation_motion_unif(&Vector3::new(0.0, 1.0, 0.0), assembly_dim),
+            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()),
+            
+            // the twist motion. more precisely: a motion that keeps the center
+            // of mass stationary and preserves the distances between the
+            // vertices to first order. this has to be the twist as long as:
+            //   - twisting is the kaleidocycle's only internal degree of
+            //     freedom
+            //   - every first-order motion of the kaleidocycle comes from an
+            //     actual motion
+            (0..N_HINGES).step_by(2).flat_map(
+                |n| {
+                    let ang_vert = ((n + 1) as f64) * PI/3.0;
+                    let vel_vert_x = 4.0 * ang_vert.cos();
+                    let vel_vert_y = 4.0 * ang_vert.sin();
+                    [
+                        DVector::from_column_slice(&[0.0, 0.0, 5.0, 0.0]),
+                        DVector::from_column_slice(&[0.0, 0.0, 1.0, 0.0]),
+                        DVector::from_column_slice(&[-vel_vert_x, -vel_vert_y, -3.0, 0.0]),
+                        DVector::from_column_slice(&[vel_vert_x, vel_vert_y, -3.0, 0.0])
+                    ]
+                }
+            ).collect::<Vec<_>>()
+        ];
+        let tangent_motions_std = tangent_motions_unif.iter().map(
+            |motion| DMatrix::from_columns(
+                &guess.column_iter().zip(motion).map(
+                    |(v, elt_motion)| local_unif_to_std(v) * elt_motion
+                ).collect::<Vec<_>>()
+            )
+        ).collect::<Vec<_>>();
+        
+        // confirm that the dimension of the tangent space is no greater than
+        // expected
+        assert_eq!(tangent.basis_std.len(), tangent_motions_unif.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.into_iter().enumerate().map(
+                |(k, v)| tangent.proj(&v.as_view(), k)
+            ).sum();
+            assert!((motion_std - motion_proj).norm_squared() < tol_sq);
+        }
+    }
+    
+    fn translation(dis: Vector3<f64>) -> DMatrix<f64> {
+        const ELEMENT_DIM: usize = 5;
+        DMatrix::from_column_slice(ELEMENT_DIM, ELEMENT_DIM, &[
+                   1.0,        0.0,        0.0, 0.0,             dis[0],
+                   0.0,        1.0,        0.0, 0.0,             dis[1],
+                   0.0,        0.0,        1.0, 0.0,             dis[2],
+            2.0*dis[0], 2.0*dis[1], 2.0*dis[2], 1.0, dis.norm_squared(),
+                   0.0,        0.0,        0.0, 0.0,                1.0
+        ])
+    }
+    
+    // confirm that projection onto a configuration subspace is equivariant with
+    // respect to Euclidean motions
+    #[test]
+    fn proj_equivar_test() {
+        // find a pair of spheres that meet at 120°
+        const SCALED_TOL: f64 = 1.0e-12;
+        let gram = {
+            let mut gram_to_be = PartialMatrix::new();
+            gram_to_be.push_sym(0, 0, 1.0);
+            gram_to_be.push_sym(1, 1, 1.0);
+            gram_to_be.push_sym(0, 1, 0.5);
+            gram_to_be
+        };
+        let guess_orig = DMatrix::from_columns(&[
+            sphere(0.0, 0.0, 0.5, 1.0),
+            sphere(0.0, 0.0, -0.5, 1.0)
+        ]);
+        let (config_orig, tangent_orig, success_orig, history_orig) = realize_gram(
+            &gram, guess_orig.clone(), &[],
+            SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
+        );
+        assert_eq!(config_orig, guess_orig);
+        assert_eq!(success_orig, true);
+        assert_eq!(history_orig.scaled_loss.len(), 1);
+        
+        // find another pair of spheres that meet at 120°. we'll think of this
+        // solution as a transformed version of the original one
+        let guess_tfm = {
+            let a = 0.5 * FRAC_1_SQRT_2;
+            DMatrix::from_columns(&[
+                sphere(a, 0.0, 7.0 + a, 1.0),
+                sphere(-a, 0.0, 7.0 - a, 1.0)
+            ])
+        };
+        let (config_tfm, tangent_tfm, success_tfm, history_tfm) = realize_gram(
+            &gram, guess_tfm.clone(), &[],
+            SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
+        );
+        assert_eq!(config_tfm, guess_tfm);
+        assert_eq!(success_tfm, true);
+        assert_eq!(history_tfm.scaled_loss.len(), 1);
+        
+        // project a nudge to the tangent space of the solution variety at the
+        // original solution
+        let motion_orig = DVector::from_column_slice(&[0.0, 0.0, 1.0, 0.0]);
+        let motion_orig_proj = tangent_orig.proj(&motion_orig.as_view(), 0);
+        
+        // project the equivalent nudge to the tangent space of the solution
+        // variety at the transformed solution
+        let motion_tfm = DVector::from_column_slice(&[FRAC_1_SQRT_2, 0.0, FRAC_1_SQRT_2, 0.0]);
+        let motion_tfm_proj = tangent_tfm.proj(&motion_tfm.as_view(), 0);
+        
+        // take the transformation that sends the original solution to the
+        // transformed solution and apply it to the motion that the original
+        // solution makes in response to the nudge
+        const ELEMENT_DIM: usize = 5;
+        let rot = DMatrix::from_column_slice(ELEMENT_DIM, ELEMENT_DIM, &[
+            FRAC_1_SQRT_2, 0.0, -FRAC_1_SQRT_2, 0.0, 0.0,
+                      0.0, 1.0,            0.0, 0.0, 0.0,
+            FRAC_1_SQRT_2, 0.0,  FRAC_1_SQRT_2, 0.0, 0.0,
+                      0.0, 0.0,            0.0, 1.0, 0.0,
+                      0.0, 0.0,            0.0, 0.0, 1.0
+        ]);
+        let transl = translation(Vector3::new(0.0, 0.0, 7.0));
+        let motion_proj_tfm = transl * rot * motion_orig_proj;
+        
+        // confirm that the projection of the nudge is equivariant. we loosen
+        // the comparison tolerance because the transformation seems to
+        // introduce some numerical error
+        const SCALED_TOL_TFM: f64 = 1.0e-9;
+        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]