app-proto/examples/kaleidocycle.rs

@@ -1,81 +0,0 @@
-use nalgebra::DMatrix;
-use std::{array, f64::consts::PI};
-
-use dyna3::engine::{Q, point, realize_gram, PartialMatrix};
-
-fn main() {
-    // set up a kaleidocycle, made of points with fixed distances between them,
-    // and find its tangent space
-    const N_POINTS: usize = 12;
-    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 = {
-        const N_HINGES: usize = 6;
-        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, &frozen,
-        1.0e-12, 0.5, 0.9, 1.1, 200, 110
-    );
-    print!("Completed Gram matrix:{}", config.tr_mul(&*Q) * &config);
-    print!("Configuration:{}", config);
-    if success {
-        println!("Target accuracy achieved!");
-    } else {
-        println!("Failed to reach target accuracy");
-    }
-    println!("Steps: {}", history.scaled_loss.len() - 1);
-    println!("Loss: {}\n", history.scaled_loss.last().unwrap());
-    
-    // find the kaleidocycle's twist motion
-    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)
-            ]
-        }
-    ).sum();
-    let normalization = 5.0 / twist_motion[(2, 0)];
-    print!("Twist motion:{}", normalization * twist_motion);
-}

app-proto/run-examples

@@ -9,4 +9,3 @@
 cargo run --example irisawa-hexlet
 cargo run --example three-spheres
 cargo run --example point-on-sphere
-cargo run --example kaleidocycle

app-proto/src/assembly.rs

@@ -5,7 +5,7 @@ use std::{collections::BTreeSet, sync::atomic::{AtomicU64, Ordering}};
 use sycamore::prelude::*;
 use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */
 
-use crate::engine::{realize_gram, local_unif_to_std, ConfigSubspace, PartialMatrix, Q};
+use crate::engine::{realize_gram, ConfigSubspace, PartialMatrix, Q};
 
 // the types of the keys we use to access an assembly's elements and constraints
 pub type ElementKey = usize;
@@ -120,7 +120,6 @@ pub struct Constraint {
     pub active: Signal<bool>
 }
 
-// the velocity is expressed in uniform coordinates
 pub struct ElementMotion<'a> {
     pub key: ElementKey,
     pub velocity: DVectorView<'a, f64>
@@ -360,14 +359,7 @@ impl Assembly {
                 // this element didn't have a column index when we started, so
                 // by invariant (2), it's unconstrained
                 let mut target_column = motion_proj.column_mut(column_index);
-                let unif_to_std = self.elements.with_untracked(
+                target_column += elt_motion.velocity;
-                    |elts| {
-                        elts[elt_motion.key].representation.with_untracked(
-                            |rep| local_unif_to_std(rep.as_view())
-                        )
-                    }
-                );
-                target_column += unif_to_std * elt_motion.velocity;
             }
         }
         

app-proto/src/display.rs

@@ -130,8 +130,6 @@ pub fn Display() -> View {
     let translate_pos_y = create_signal(0.0);
     let translate_neg_z = create_signal(0.0);
     let translate_pos_z = create_signal(0.0);
-    let shrink_neg = create_signal(0.0);
-    let shrink_pos = create_signal(0.0);
     
     // change listener
     let scene_changed = create_signal(true);
@@ -166,7 +164,6 @@ pub fn Display() -> View {
         
         // manipulation
         const TRANSLATION_SPEED: f64 = 0.15; // in length units per second
-        const SHRINKING_SPEED: f64 = 0.15; // in length units per second
         
         // display parameters
         const OPACITY: f32 = 0.5; /* SCAFFOLDING */
@@ -295,8 +292,6 @@ pub fn Display() -> View {
             let translate_pos_y_val = translate_pos_y.get();
             let translate_neg_z_val = translate_neg_z.get();
             let translate_pos_z_val = translate_pos_z.get();
-            let shrink_neg_val = shrink_neg.get();
-            let shrink_pos_val = shrink_pos.get();
             
             // update the assembly's orientation
             let ang_vel = {
@@ -328,27 +323,24 @@ pub fn Display() -> View {
                 let sel = state.selection.with(
                     |sel| *sel.into_iter().next().unwrap()
                 );
+                let rep = state.assembly.elements.with_untracked(
+                    |elts| elts[sel].representation.get_clone_untracked()
+                );
                 let translate_x = translate_pos_x_val - translate_neg_x_val;
                 let translate_y = translate_pos_y_val - translate_neg_y_val;
                 let translate_z = translate_pos_z_val - translate_neg_z_val;
-                let shrink = shrink_pos_val - shrink_neg_val;
+                if translate_x != 0.0 || translate_y != 0.0 || translate_z != 0.0 {
-                let translating =
+                    let vel_field = {
-                    translate_x != 0.0
+                        let u = Vector3::new(translate_x, translate_y, translate_z).normalize();
-                    || translate_y != 0.0
+                        DMatrix::from_column_slice(5, 5, &[
-                    || translate_z != 0.0;
+                                 0.0,      0.0,      0.0, 0.0, u[0],
-                if translating || shrink != 0.0 {
+                                 0.0,      0.0,      0.0, 0.0, u[1],
-                    let elt_motion = {
+                                 0.0,      0.0,      0.0, 0.0, u[2],
-                        let u = if translating {
+                            2.0*u[0], 2.0*u[1], 2.0*u[2], 0.0,  0.0,
-                            TRANSLATION_SPEED * Vector3::new(
+                                 0.0,      0.0,      0.0, 0.0,  0.0
-                                translate_x, translate_y, translate_z
+                        ])
-                            ).normalize()
-                        } else {
-                            Vector3::zeros()
-                        };
-                        time_step * DVector::from_column_slice(
-                            &[u[0], u[1], u[2], SHRINKING_SPEED * shrink, 0.0]
-                        )
                     };
+                    let elt_motion: DVector<f64> = time_step * TRANSLATION_SPEED * vel_field * rep;
                     assembly_for_raf.deform(
                         vec![
                             ElementMotion {
@@ -509,8 +501,6 @@ pub fn Display() -> View {
             "s" | "S" if shift => translate_pos_z.set(value),
             "w" | "W" => translate_pos_y.set(value),
             "s" | "S" => translate_neg_y.set(value),
-            "]" | "}" => shrink_neg.set(value),
-            "[" | "{" => shrink_pos.set(value),
             _ => manipulating = false
         };
         if manipulating {

app-proto/src/engine.rs

@@ -90,33 +90,32 @@ impl PartialMatrix {
 #[derive(Clone)]
 pub struct ConfigSubspace {
     assembly_dim: usize,
-    basis_std: Vec<DMatrix<f64>>,
+    basis: Vec<DMatrix<f64>>
-    basis_proj: Vec<DMatrix<f64>>
 }
 
 impl ConfigSubspace {
     pub fn zero(assembly_dim: usize) -> ConfigSubspace {
         ConfigSubspace {
             assembly_dim: assembly_dim,
-            basis_proj: Vec::new(),
+            basis: 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>, proj_to_std: DMatrix<f64>, assembly_dim: usize) -> ConfigSubspace {
+    fn symmetric_kernel(a: 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_std = DMatrix::from_columns(
+        let basis = eig_pairs.filter_map(
-            eig_pairs.filter_map(
+            |(λ, v)| (λ.abs() < THRESHOLD).then_some(
-                |(λ, v)| (λ.abs() < THRESHOLD).then_some(v)
+                Into::<DMatrix<f64>>::into(
-            ).collect::<Vec<_>>().as_slice()
+                    v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim))
+                )
+            )
         );
         
         /* DEBUG */
@@ -126,53 +125,30 @@ 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;
-        
-        // 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)
-        ));
-        
         ConfigSubspace {
             assembly_dim: assembly_dim,
-            basis_std: basis_std_orth.column_iter().map(
+            basis: basis.collect()
-                |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_std.len()
+        self.basis.len()
     }
     
     pub fn assembly_dim(&self) -> usize {
         self.assembly_dim
     }
     
-    // find the projection onto this subspace of the motion where the element
+    // find the projection onto this subspace, with respect to the Euclidean
-    // with the given column index has velocity `v`. the velocity is given in
+    // inner product, of the motion where the element with the given column
-    // projection coordinates, and the projection is done with respect to the
+    // index has velocity `v`
-    // 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_proj.iter().zip(self.basis_std.iter()).map(
+            self.basis.iter().map(
-                |(b_proj, b_std)| b_proj.column(column_index).dot(&v) * b_std
+                |b| b.column(column_index).dot(&v) * b
             ).sum()
         }
     }
@@ -239,38 +215,6 @@ fn basis_matrix(index: (usize, usize), nrows: usize, ncols: usize) -> DMatrix<f6
     result
 }
 
-// given a normalized vector `v` representing an element, build a basis for the
-// element's linear configuration space consisting of:
-// - the unit translation motions of the element
-// - the unit shrinking motion of the element, if it's a sphere
-// - one or two vectors whose coefficients vanish on the tangent space of the
-//   normalization variety
-pub fn local_unif_to_std(v: DVectorView<f64>) -> DMatrix<f64> {
-    const ELEMENT_DIM: usize = 5;
-    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, &[
-            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, &[
-                 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,
@@ -400,17 +344,7 @@ pub fn realize_gram(
     }
     let success = state.loss < tol;
     let tangent = if success {
-        // express the uniform basis in the standard basis
+        ConfigSubspace::symmetric_kernel(hess, assembly_dim)
-        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)
     };