app-proto/examples/kaleidocycle.rs

@@ -1,4 +1,4 @@
-use nalgebra::{DMatrix, DVector};
+use nalgebra::DMatrix;
 use std::{array, f64::consts::PI};
 
 use dyna3::engine::{Q, point, realize_gram, PartialMatrix};
@@ -59,13 +59,22 @@ 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| [
+        |n| {
-            tangent.proj(&up.as_view(), n),
+            let up_field = {
-            tangent.proj(&down.as_view(), n+1)
+                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)
             ]
+        }
     ).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]
+                            &[u[0], u[1], u[2], SHRINKING_SPEED * shrink, 0.0]
                         )
                     };
                     assembly_for_raf.deform(

app-proto/src/engine.rs

@@ -107,14 +107,13 @@ 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. the basis is expressed in the projection
+        // find a basis for the kernel, expressed in the standard coordinates
-        // coordinates, and it's orthonormal with respect to the projection
+        const ELEMENT_DIM: usize = 5;
-        // inner product
+        const THRESHOLD: f64 = 1.0e-4;
-        const THRESHOLD: f64 = 0.1;
+        let eig = SymmetricEigen::new(a);
-        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_proj = DMatrix::from_columns(
+        let basis_std = DMatrix::from_columns(
             eig_pairs.filter_map(
                 |(λ, v)| (λ.abs() < THRESHOLD).then_some(v)
             ).collect::<Vec<_>>().as_slice()
@@ -127,27 +126,29 @@ impl ConfigSubspace {
             format!("Eigenvalues used to find kernel:{}", eig.eigenvalues)
         ));
         
-        // express the basis in the standard coordinates
+        // express the basis in the projection coordinates
-        let basis_std = proj_to_std * &basis_proj;
+        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;
         
         // 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)
+            format!("Basis in projection coordinates:{}", basis_proj_orth)
         ));
         
-        const ELEMENT_DIM: usize = 5;
-        const UNIFORM_DIM: usize = 4;
         ConfigSubspace {
             assembly_dim: assembly_dim,
-            basis_std: basis_std.column_iter().map(
+            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.column_iter().map(
+            basis_proj: basis_proj_orth.column_iter().map(
                 |v| Into::<DMatrix<f64>>::into(
-                    v.reshape_generic(Dyn(UNIFORM_DIM), Dyn(assembly_dim))
+                    v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim))
                 )
             ).collect()
         }
@@ -246,25 +247,26 @@ 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, UNIFORM_DIM, &[
+        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],
+            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, UNIFORM_DIM, &[
+        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
+            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]
         ])
     }
 }
@@ -399,13 +401,11 @@ pub fn realize_gram(
     let success = state.loss < tol;
     let tangent = if success {
         // express the uniform basis in the standard basis
-        const UNIFORM_DIM: usize = 4;
+        let mut unif_to_std = DMatrix::<f64>::zeros(total_dim, total_dim);
-        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, UNIFORM_DIM * n);
+            let block_start = element_dim * n;
             unif_to_std
-                .view_mut(block_start, (element_dim, UNIFORM_DIM))
+                .view_mut((block_start, block_start), (element_dim, element_dim))
                 .copy_from(&local_unif_to_std(state.config.column(n)));
         }