Do symmetric_kernel in projection coordinates

Instead of finding the kernel in the standard coordinates and then
expressing it in the projection coordinates, work in the projection
coordinates from the beginning by applying a change of basis to the
input matrix.

app-proto/examples/kaleidocycle.rs

@@ -1,4 +1,4 @@
-use nalgebra::DMatrix;
+use nalgebra::{DMatrix, DVector};
 use std::{array, f64::consts::PI};
 
 use dyna3::engine::{Q, point, realize_gram, PartialMatrix};
@@ -59,22 +59,13 @@ fn main() {
     println!("Loss: {}\n", history.scaled_loss.last().unwrap());
     
     // find the kaleidocycle's twist motion
+    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(
-        |n| {
-            let up_field = {
-                DMatrix::from_column_slice(5, 5, &[
-                    0.0, 0.0, 0.0, 0.0, 0.0,
-                    0.0, 0.0, 0.0, 0.0, 0.0,
-                    0.0, 0.0, 0.0, 0.0, 1.0,
-                    0.0, 0.0, 2.0, 0.0, 0.0,
-                    0.0, 0.0, 0.0, 0.0, 0.0
-                ])
-            };
-            [
-                tangent.proj(&(&up_field * config.column(n)).as_view(), n),
-                tangent.proj(&(-&up_field * config.column(n+1)).as_view(), n+1)
-            ]
-        }
+        |n| [
+            tangent.proj(&up.as_view(), n),
+            tangent.proj(&down.as_view(), n+1)
+        ]
     ).sum();
     let normalization = 5.0 / twist_motion[(2, 0)];
     print!("Twist motion:{}", normalization * twist_motion);

app-proto/src/display.rs

@@ -346,7 +346,7 @@ pub fn Display() -> View {
                             Vector3::zeros()
                         };
                         time_step * DVector::from_column_slice(
-                            &[u[0], u[1], u[2], SHRINKING_SPEED * shrink, 0.0]
+                            &[u[0], u[1], u[2], SHRINKING_SPEED * shrink]
                         )
                     };
                     assembly_for_raf.deform(

app-proto/src/engine.rs

@@ -107,13 +107,14 @@ impl ConfigSubspace {
     // 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>, 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);
+        // find a basis for the kernel. the basis is expressed in the projection
+        // coordinates, and it's orthonormal with respect to the projection
+        // inner product
+        const THRESHOLD: f64 = 0.1;
+        let eig = SymmetricEigen::new(proj_to_std.tr_mul(&a) * &proj_to_std);
         let eig_vecs = eig.eigenvectors.column_iter();
         let eig_pairs = eig.eigenvalues.iter().zip(eig_vecs);
-        let basis_std = DMatrix::from_columns(
+        let basis_proj = DMatrix::from_columns(
             eig_pairs.filter_map(
                 |(λ, v)| (λ.abs() < THRESHOLD).then_some(v)
             ).collect::<Vec<_>>().as_slice()
@@ -126,29 +127,27 @@ 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;
+        // express the basis in the standard coordinates
+        let basis_std = proj_to_std * &basis_proj;
         
         // print the projection basis in projection coordinates
         #[cfg(all(target_family = "wasm", target_os = "unknown"))]
         console::log_1(&JsValue::from(
-            format!("Basis in projection coordinates:{}", basis_proj_orth)
+            format!("Basis in projection coordinates:{}", basis_proj)
         ));
         
+        const ELEMENT_DIM: usize = 5;
+        const UNIFORM_DIM: usize = 4;
         ConfigSubspace {
             assembly_dim: assembly_dim,
-            basis_std: basis_std_orth.column_iter().map(
+            basis_std: basis_std.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(
+            basis_proj: basis_proj.column_iter().map(
                 |v| Into::<DMatrix<f64>>::into(
-                    v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim))
+                    v.reshape_generic(Dyn(UNIFORM_DIM), Dyn(assembly_dim))
                 )
             ).collect()
         }
@@ -247,26 +246,25 @@ fn basis_matrix(index: (usize, usize), nrows: usize, ncols: usize) -> DMatrix<f6
 //   normalization variety
 pub fn local_unif_to_std(v: DVectorView<f64>) -> DMatrix<f64> {
     const ELEMENT_DIM: usize = 5;
+    const UNIFORM_DIM: usize = 4;
     let curv = 2.0*v[3];
     if v.dot(&(&*Q * v)) < 0.5 {
         // `v` represents a point. the normalization condition says that the
         // curvature component of `v` is 1/2
-        DMatrix::from_column_slice(ELEMENT_DIM, ELEMENT_DIM, &[
+        DMatrix::from_column_slice(ELEMENT_DIM, UNIFORM_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],
-            v[0], v[1], v[2], v[3], v[4],
              0.0,  0.0,  0.0,  0.0,  1.0
         ])
     } else {
         // `v` represents a sphere. the normalization condition says that the
         // Lorentz product of `v` with itself is 1
-        DMatrix::from_column_slice(ELEMENT_DIM, ELEMENT_DIM, &[
+        DMatrix::from_column_slice(ELEMENT_DIM, UNIFORM_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]
+            curv*v[0], curv*v[1], curv*v[2], curv*v[3], curv*v[4] + 1.0
         ])
     }
 }
@@ -401,11 +399,13 @@ pub fn realize_gram(
     let success = state.loss < tol;
     let tangent = if success {
         // express the uniform basis in the standard basis
-        let mut unif_to_std = DMatrix::<f64>::zeros(total_dim, total_dim);
+        const UNIFORM_DIM: usize = 4;
+        let total_dim_unif = UNIFORM_DIM * assembly_dim;
+        let mut unif_to_std = DMatrix::<f64>::zeros(total_dim, total_dim_unif);
         for n in 0..assembly_dim {
-            let block_start = element_dim * n;
+            let block_start = (element_dim * n, UNIFORM_DIM * n);
             unif_to_std
-                .view_mut((block_start, block_start), (element_dim, element_dim))
+                .view_mut(block_start, (element_dim, UNIFORM_DIM))
                 .copy_from(&local_unif_to_std(state.config.column(n)));
         }