From 1e3505dd01540af0005be10bdc6d6b66ab33c3f0 Mon Sep 17 00:00:00 2001
From: Aaron Fenyes <aaron.fenyes@fareycircles.ooo>
Date: Thu, 23 Jan 2025 10:16:04 -0800
Subject: [PATCH] Confirm that projection is Euclidean-equivariant

---
 app-proto/src/engine.rs | 88 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 88 insertions(+)

diff --git a/app-proto/src/engine.rs b/app-proto/src/engine.rs
index 4c90156..0756fac 100644
--- a/app-proto/src/engine.rs
+++ b/app-proto/src/engine.rs
@@ -490,6 +490,9 @@ pub mod irisawa {
 
 #[cfg(test)]
 mod tests {
+    use nalgebra::Vector3;
+    use std::f64::consts::FRAC_1_SQRT_2;
+    
     use super::{*, irisawa::realize_irisawa_hexlet};
     
     #[test]
@@ -617,6 +620,91 @@ mod tests {
         }
     }
     
+    fn translation(u: 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,             u[0],
+                 0.0,      1.0,      0.0, 0.0,             u[1],
+                 0.0,      0.0,      1.0, 0.0,             u[2],
+            2.0*u[0], 2.0*u[1], 2.0*u[2], 1.0, u.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]