app-proto/src/assembly.rs

@@ -33,9 +33,8 @@ pub struct Element {
     pub serial: u64,
     
     // the configuration matrix column index that was assigned to this element
-    // last time the assembly was realized, or `None` if the element has never
-    // been through a realization
-    column_index: Option<usize>
+    // last time the assembly was realized
+    column_index: usize
 }
 
 impl Element {
@@ -63,7 +62,7 @@ impl Element {
             representation: create_signal(representation),
             constraints: create_signal(BTreeSet::default()),
             serial: serial,
-            column_index: None
+            column_index: 0
         }
     }
     
@@ -120,13 +119,6 @@ pub struct Constraint {
     pub active: Signal<bool>
 }
 
-pub struct ElementMotion<'a> {
-    pub key: ElementKey,
-    pub velocity: DVectorView<'a, f64>
-}
-
-type AssemblyMotion<'a> = Vec<ElementMotion<'a>>;
-
 // a complete, view-independent description of an assembly
 #[derive(Clone)]
 pub struct Assembly {
@@ -134,16 +126,7 @@ pub struct Assembly {
     pub elements: Signal<Slab<Element>>,
     pub constraints: Signal<Slab<Constraint>>,
     
-    // solution variety tangent space. the basis vectors are stored in
-    // configuration matrix format, ordered according to the elements' column
-    // indices. when you realize the assembly, every element that's present
-    // during realization gets a column index and is reflected in the tangent
-    // space. since the methods in this module never assign column indices
-    // without later realizing the assembly, we get the following invariant:
-    //
-    //  (1) if an element has a column index, its tangent motions can be found
-    //      in that column of the tangent space basis matrices
-    //
+    // solution variety tangent space
     pub tangent: Signal<ConfigSubspace>,
     
     // indexing
@@ -155,7 +138,7 @@ impl Assembly {
         Assembly {
             elements: create_signal(Slab::new()),
             constraints: create_signal(Slab::new()),
-            tangent: create_signal(ConfigSubspace::zero(0)),
+            tangent: create_signal(ConfigSubspace::zero()),
             elements_by_id: create_signal(FxHashMap::default())
         }
     }
@@ -169,6 +152,13 @@ impl Assembly {
         let id = elt.id.clone();
         let key = self.elements.update(|elts| elts.insert(elt));
         self.elements_by_id.update(|elts_by_id| elts_by_id.insert(id, key));
+        
+        // realize to update the tangent space
+        /* KLUDGE */
+        // since the newly inserted element is unconstrained, we should be able
+        // to update the tangent space without recomputing the Hessian and its
+        // eigendecomposition
+        self.realize();
     }
     
     pub fn try_insert_element(&self, elt: Element) -> bool {
@@ -219,7 +209,7 @@ impl Assembly {
         // index the elements
         self.elements.update_silent(|elts| {
             for (index, (_, elt)) in elts.into_iter().enumerate() {
-                elt.column_index = Some(index);
+                elt.column_index = index;
             }
         });
         
@@ -231,8 +221,8 @@ impl Assembly {
                 for (_, cst) in csts {
                     if cst.active.get_untracked() && cst.lorentz_prod_valid.get_untracked() {
                         let subjects = cst.subjects;
-                        let row = elts[subjects.0].column_index.unwrap();
-                        let col = elts[subjects.1].column_index.unwrap();
+                        let row = elts[subjects.0].column_index;
+                        let col = elts[subjects.1].column_index;
                         gram_to_be.push_sym(row, col, cst.lorentz_prod.get_untracked());
                     }
                 }
@@ -242,7 +232,7 @@ impl Assembly {
             // Gram matrix
             let mut guess_to_be = DMatrix::<f64>::zeros(5, elts.len());
             for (_, elt) in elts {
-                let index = elt.column_index.unwrap();
+                let index = elt.column_index;
                 gram_to_be.push_sym(index, index, 1.0);
                 guess_to_be.set_column(index, &elt.representation.get_clone_untracked());
             }
@@ -289,7 +279,7 @@ impl Assembly {
             // read out the solution
             for (_, elt) in self.elements.get_clone_untracked() {
                 elt.representation.update(
-                    |rep| rep.set_column(0, &config.column(elt.column_index.unwrap()))
+                    |rep| rep.set_column(0, &config.column(elt.column_index))
                 );
             }
             
@@ -300,67 +290,26 @@ impl Assembly {
     
     // --- deformation ---
     
-    // project the given motion to the tangent space of the solution variety and
-    // move the assembly along it. the implementation is based on invariant (1)
-    // from above and the following additional invariant:
-    //
-    //  (2) if an element is affected by a constraint, it has a column index
-    //
-    // we have this invariant because the assembly gets realized each time you
-    // add a constraint
-    pub fn deform(&self, motion: AssemblyMotion) {
+    pub fn deform(&self, element_motions: Vec<(ElementKey, DVectorView<f64>)>) {
         /* KLUDGE */
-        // when the tangent space is zero, deformation won't do anything, but
-        // the attempt to deform should be registered in the UI. this console
-        // message will do for now
-        if self.tangent.with(|tan| tan.dim() <= 0 && tan.assembly_dim() > 0) {
+        // when the tangent space is zero, we currently need to avoid calling
+        // its `proj` method, because it will panic rather than returning zero.
+        // in the future, we'll want a more intentionally designed system for
+        // handling this case
+        if self.tangent.with(|tan| tan.dim() <= 0) {
             console::log_1(&JsValue::from("The assembly is rigid"));
+            return;
         }
         
-        // give a column index to each moving element that doesn't have one yet.
-        // this temporarily breaks invariant (1), but the invariant will be
-        // restored when we realize the assembly at the end of the deformation.
-        // in the process, we find out how many matrix columns we'll need to
-        // hold the deformation
-        let realized_dim = self.tangent.with(|tan| tan.assembly_dim());
-        let motion_dim = self.elements.update_silent(|elts| {
-            let mut next_column_index = realized_dim;
-            for elt_motion in motion.iter() {
-                let moving_elt = &mut elts[elt_motion.key];
-                if moving_elt.column_index.is_none() {
-                    moving_elt.column_index = Some(next_column_index);
-                    next_column_index += 1;
-                }
-            }
-            next_column_index
-        });
+        const ELEMENT_DIM: usize = 5;
+        let assembly_dim = self.elements.with(|elts| elts.len());
+        let mut motion_proj = DMatrix::zeros(ELEMENT_DIM, assembly_dim);
         
         // project the element motions onto the tangent space of the solution
-        // variety and sum them to get a deformation of the whole assembly. the
-        // matrix `motion_proj` that holds the deformation has extra columns for
-        // any moving elements that aren't reflected in the saved tangent space
-        const ELEMENT_DIM: usize = 5;
-        let mut motion_proj = DMatrix::zeros(ELEMENT_DIM, motion_dim);
-        for elt_motion in motion {
-            // we can unwrap the column index because we know that every moving
-            // element has one at this point
-            let column_index = self.elements.with_untracked(
-                |elts| elts[elt_motion.key].column_index.unwrap()
-            );
-            
-            if column_index < realized_dim {
-                // this element had a column index when we started, so by
-                // invariant (1), it's reflected in the tangent space
-                let mut target_columns = motion_proj.columns_mut(0, realized_dim);
-                target_columns += self.tangent.with(
-                    |tan| tan.proj(&elt_motion.velocity, column_index)
-                );
-            } else {
-                // 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);
-                target_column += elt_motion.velocity;
-            }
+        // variety, and sum them to get a deformation of the whole assembly
+        for (elt_key, elt_motion) in element_motions {
+            let column_index = self.elements.with(|elts| elts[elt_key].column_index);
+            motion_proj += self.tangent.with(|tan| tan.proj(&elt_motion, column_index));
         }
         
         // step each element along the mass shell geodesic that matches its
@@ -370,24 +319,13 @@ impl Assembly {
         // element is on the 1 mass shell
         for (_, elt) in self.elements.get_clone_untracked() {
             elt.representation.update_silent(|rep| {
-                match elt.column_index {
-                    Some(column_index) => {
-                        let rep_next = &*rep + motion_proj.column(column_index);
+                let rep_next = &*rep + motion_proj.column(elt.column_index);
                 let normalizer = rep_next.dot(&(&*Q * &rep_next));
                 rep.set_column(0, &(rep_next / normalizer));
-                    },
-                    None => {
-                        console::log_1(&JsValue::from(
-                            format!("No velocity to unpack for fresh element \"{}\"", elt.id)
-                        ))
-                    }
-                };
             });
         }
         
-        // bring the configuration back onto the solution variety. this also
-        // gets the elements' column indices and the saved tangent space back in
-        // sync
+        // bring the configuration back onto the solution variety
         self.realize();
     }
 }

app-proto/src/display.rs

@@ -14,7 +14,7 @@ use web_sys::{
     wasm_bindgen::{JsCast, JsValue}
 };
 
-use crate::{AppState, assembly::{ElementKey, ElementMotion}};
+use crate::{AppState, assembly::ElementKey};
 
 fn compile_shader(
     context: &WebGl2RenderingContext,
@@ -341,14 +341,7 @@ pub fn Display() -> View {
                         ])
                     };
                     let elt_motion: DVector<f64> = time_step * TRANSLATION_SPEED * vel_field * rep;
-                    assembly_for_raf.deform(
-                        vec![
-                            ElementMotion {
-                                key: sel,
-                                velocity: elt_motion.as_view()
-                            }
-                        ]
-                    );
+                    assembly_for_raf.deform(vec![(sel, elt_motion.as_view())]);
                     scene_changed.set(true);
                 }
             }

app-proto/src/engine.rs

@@ -88,17 +88,11 @@ impl PartialMatrix {
 // --- configuration subspaces ---
 
 #[derive(Clone)]
-pub struct ConfigSubspace {
-    assembly_dim: usize,
-    basis: Vec<DMatrix<f64>>
-}
+pub struct ConfigSubspace(Vec<DMatrix<f64>>);
 
 impl ConfigSubspace {
-    pub fn zero(assembly_dim: usize) -> ConfigSubspace {
-        ConfigSubspace {
-            assembly_dim: assembly_dim,
-            basis: Vec::new()
-        }
+    pub fn zero() -> ConfigSubspace {
+        ConfigSubspace(Vec::new())
     }
     
     // approximate the kernel of a symmetric endomorphism of the configuration
@@ -125,32 +119,25 @@ impl ConfigSubspace {
             format!("Eigenvalues used to find kernel: {}", eig.eigenvalues)
         ));
         
-        ConfigSubspace {
-            assembly_dim: assembly_dim,
-            basis: basis.collect()
-        }
+        ConfigSubspace(basis.collect())
     }
     
     pub fn dim(&self) -> usize {
-        self.basis.len()
-    }
-    
-    pub fn assembly_dim(&self) -> usize {
-        self.assembly_dim
+        let ConfigSubspace(basis) = self;
+        basis.len()
     }
     
     // find the projection onto this subspace of the motion where the element
     // with the given column index has velocity `v`
+    /* TO DO */
+    // for the zero subspace, this method's behavior doesn't match its name: it
+    // panics rather than returning zero
     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(
+        let ConfigSubspace(basis) = self;
+        basis.into_iter().map(
             |b| b.column(column_index).dot(&v) * b
         ).sum()
     }
-    }
 }
 
 // --- descent history ---
@@ -338,14 +325,14 @@ pub fn realize_gram(
                 state = better_state;
                 history.backoff_steps.push(backoff_steps);
             },
-            None => return (state.config, ConfigSubspace::zero(assembly_dim), false, history)
+            None => return (state.config, ConfigSubspace::zero(), false, history)
         };
     }
     let success = state.loss < tol;
     let tangent = if success {
         ConfigSubspace::symmetric_kernel(hess, assembly_dim)
     } else {
-        ConfigSubspace::zero(assembly_dim)
+        ConfigSubspace::zero()
     };
     (state.config, tangent, success, history)
 }