Try projecting with the uniform inner product

This projection should be invariant under Euclidean motions. For now, we
assume that all of the elements are spheres.

app-proto/src/engine.rs

@@ -90,32 +90,33 @@ impl PartialMatrix {
 #[derive(Clone)]
 pub struct ConfigSubspace {
     assembly_dim: usize,
-    basis: Vec<DMatrix<f64>>
+    basis_std: Vec<DMatrix<f64>>,
+    basis_proj: Vec<DMatrix<f64>>
 }
 
 impl ConfigSubspace {
     pub fn zero(assembly_dim: usize) -> ConfigSubspace {
         ConfigSubspace {
             assembly_dim: assembly_dim,
-            basis: Vec::new()
+            basis_proj: Vec::new(),
+            basis_std: Vec::new()
         }
     }
     
     // approximate the kernel of a symmetric endomorphism of the configuration
     // space for `assembly_dim` elements. we consider an eigenvector to be part
     // of the kernel if its eigenvalue is smaller than the constant `THRESHOLD`
-    fn symmetric_kernel(a: DMatrix<f64>, assembly_dim: usize) -> ConfigSubspace {
+    fn symmetric_kernel(a: DMatrix<f64>, proj_to_std: DMatrix<f64>, assembly_dim: usize) -> ConfigSubspace {
+        // find a basis for the kernel, expressed in the standard coordinates
         const ELEMENT_DIM: usize = 5;
         const THRESHOLD: f64 = 1.0e-4;
         let eig = SymmetricEigen::new(a);
         let eig_vecs = eig.eigenvectors.column_iter();
         let eig_pairs = eig.eigenvalues.iter().zip(eig_vecs);
-        let basis = eig_pairs.filter_map(
-            |(λ, v)| (λ.abs() < THRESHOLD).then_some(
-                Into::<DMatrix<f64>>::into(
-                    v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim))
-                )
-            )
+        let basis_std = DMatrix::from_columns(
+            eig_pairs.filter_map(
+                |(λ, v)| (λ.abs() < THRESHOLD).then_some(v)
+            ).collect::<Vec<_>>().as_slice()
         );
         
         /* DEBUG */
@@ -125,30 +126,50 @@ impl ConfigSubspace {
             format!("Eigenvalues used to find kernel:{}", eig.eigenvalues)
         ));
         
+        // express the basis in the projection coordinates
+        let basis_proj = proj_to_std.clone().qr().solve(&basis_std).unwrap();
+        
+        // orthonormalize the basis with respect to the projection inner product
+        let basis_proj_orth = basis_proj.qr().q();
+        let basis_std_orth = proj_to_std * &basis_proj_orth;
+        console::log_1(&JsValue::from(
+            format!("Basis in projection coordinates:{}", basis_proj_orth)
+        ));
+        
         ConfigSubspace {
             assembly_dim: assembly_dim,
-            basis: basis.collect()
+            basis_std: basis_std_orth.column_iter().map(
+                |v| Into::<DMatrix<f64>>::into(
+                    v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim))
+                )
+            ).collect(),
+            basis_proj: basis_proj_orth.column_iter().map(
+                |v| Into::<DMatrix<f64>>::into(
+                    v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim))
+                )
+            ).collect()
         }
     }
     
     pub fn dim(&self) -> usize {
-        self.basis.len()
+        self.basis_std.len()
     }
     
     pub fn assembly_dim(&self) -> usize {
         self.assembly_dim
     }
     
-    // find the projection onto this subspace, with respect to the Euclidean
-    // inner product, of the motion where the element with the given column
-    // index has velocity `v`
+    // find the projection onto this subspace of the motion where the element
+    // with the given column index has velocity `v`. the velocity is given in
+    // projection coordinates, and the projection is done with respect to the
+    // projection inner product
     pub fn proj(&self, v: &DVectorView<f64>, column_index: usize) -> DMatrix<f64> {
         if self.dim() == 0 {
             const ELEMENT_DIM: usize = 5;
             DMatrix::zeros(ELEMENT_DIM, self.assembly_dim)
         } else {
-            self.basis.iter().map(
-                |b| b.column(column_index).dot(&v) * b
+            self.basis_proj.iter().zip(self.basis_std.iter()).map(
+                |(b_proj, b_std)| b_proj.column(column_index).dot(&v) * b_std
             ).sum()
         }
     }
@@ -215,6 +236,21 @@ fn basis_matrix(index: (usize, usize), nrows: usize, ncols: usize) -> DMatrix<f6
     result
 }
 
+// given a spacelike unit vector `v`, which represents a sphere, build the basis
+// for the configuration space given by the three unit translation motions of
+// the sphere, the unit shrinking motion of the sphere, and `v`
+pub fn local_unif_to_std(v: DVectorView<f64>) -> DMatrix<f64> {
+    const ELEMENT_DIM: usize = 5;
+    let curv = 2.0*v[3];
+    DMatrix::from_column_slice(ELEMENT_DIM, ELEMENT_DIM, &[
+             curv,       0.0,       0.0,       0.0,            v[0],
+              0.0,      curv,       0.0,       0.0,            v[1],
+              0.0,       0.0,      curv,       0.0,            v[2],
+        curv*v[0], curv*v[1], curv*v[2], curv*v[3], curv*v[4] + 1.0,
+             v[0],      v[1],      v[2],      v[3],            v[4]
+    ])
+}
+
 // use backtracking line search to find a better configuration
 fn seek_better_config(
     gram: &PartialMatrix,
@@ -344,7 +380,17 @@ pub fn realize_gram(
     }
     let success = state.loss < tol;
     let tangent = if success {
-        ConfigSubspace::symmetric_kernel(hess, assembly_dim)
+        // express the uniform basis in the standard basis
+        let mut unif_to_std = DMatrix::<f64>::zeros(total_dim, total_dim);
+        for n in 0..assembly_dim {
+            let block_start = element_dim * n;
+            unif_to_std
+                .view_mut((block_start, block_start), (element_dim, element_dim))
+                .copy_from(&local_unif_to_std(state.config.column(n)));
+        }
+        
+        // find the kernel of the Hessian. give it the uniform inner product
+        ConfigSubspace::symmetric_kernel(hess, unif_to_std, assembly_dim)
     } else {
         ConfigSubspace::zero(assembly_dim)
     };