Remove trailing space in console log

Also, correct the check for whether an element had a column index when
we started. The previous revision would've gotten the wrong answer for
an element without a column index that appeared more than once in the
motion.

app-proto/examples/irisawa-hexlet.rs

@@ -2,7 +2,7 @@ use dyna3::engine::{Q, irisawa::realize_irisawa_hexlet};
 
 fn main() {
     const SCALED_TOL: f64 = 1.0e-12;
-    let (config, success, history) = realize_irisawa_hexlet(SCALED_TOL);
+    let (config, _, success, history) = realize_irisawa_hexlet(SCALED_TOL);
     print!("\nCompleted Gram matrix:{}", config.tr_mul(&*Q) * &config);
     if success {
         println!("Target accuracy achieved!");

app-proto/examples/point-on-sphere.rs

@@ -18,7 +18,7 @@ fn main() {
     ]);
     let frozen = [(3, 0)];
     println!();
-    let (config, success, history) = realize_gram(
+    let (config, _, success, history) = realize_gram(
         &gram, guess, &frozen,
         1.0e-12, 0.5, 0.9, 1.1, 200, 110
     );

app-proto/examples/three-spheres.rs

@@ -21,7 +21,7 @@ fn main() {
         ])
     };
     println!();
-    let (config, success, history) = realize_gram(
+    let (config, _, success, history) = realize_gram(
         &gram, guess, &[],
         1.0e-12, 0.5, 0.9, 1.1, 200, 110
     );

app-proto/src/assembly.rs

@@ -1,11 +1,11 @@
-use nalgebra::{DMatrix, DVector, Vector3};
+use nalgebra::{DMatrix, DVector, DVectorView, Vector3};
 use rustc_hash::FxHashMap;
 use slab::Slab;
 use std::{collections::BTreeSet, sync::atomic::{AtomicU64, Ordering}};
 use sycamore::prelude::*;
 use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */
 
-use crate::engine::{realize_gram, PartialMatrix};
+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;
@@ -33,8 +33,9 @@ pub struct Element {
     pub serial: u64,
     
     // the configuration matrix column index that was assigned to this element
-    // last time the assembly was realized
+    // last time the assembly was realized, or `None` if the element has never
-    column_index: usize
+    // been through a realization
+    column_index: Option<usize>
 }
 
 impl Element {
@@ -62,7 +63,7 @@ impl Element {
             representation: create_signal(representation),
             constraints: create_signal(BTreeSet::default()),
             serial: serial,
-            column_index: 0
+            column_index: None
         }
     }
     
@@ -109,7 +110,6 @@ impl Element {
         }
     }
 }
-            
 
 #[derive(Clone)]
 pub struct Constraint {
@@ -120,6 +120,13 @@ 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 {
@@ -127,6 +134,18 @@ 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
+    //
+    pub tangent: Signal<ConfigSubspace>,
+    
     // indexing
     pub elements_by_id: Signal<FxHashMap<String, ElementKey>>
 }
@@ -136,6 +155,7 @@ impl Assembly {
         Assembly {
             elements: create_signal(Slab::new()),
             constraints: create_signal(Slab::new()),
+            tangent: create_signal(ConfigSubspace::zero(0)),
             elements_by_id: create_signal(FxHashMap::default())
         }
     }
@@ -199,7 +219,7 @@ impl Assembly {
         // index the elements
         self.elements.update_silent(|elts| {
             for (index, (_, elt)) in elts.into_iter().enumerate() {
-                elt.column_index = index;
+                elt.column_index = Some(index);
             }
         });
         
@@ -211,8 +231,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;
+                        let row = elts[subjects.0].column_index.unwrap();
-                        let col = elts[subjects.1].column_index;
+                        let col = elts[subjects.1].column_index.unwrap();
                         gram_to_be.push_sym(row, col, cst.lorentz_prod.get_untracked());
                     }
                 }
@@ -222,7 +242,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;
+                let index = elt.column_index.unwrap();
                 gram_to_be.push_sym(index, index, 1.0);
                 guess_to_be.set_column(index, &elt.representation.get_clone_untracked());
             }
@@ -247,7 +267,7 @@ impl Assembly {
         }
         
         // look for a configuration with the given Gram matrix
-        let (config, success, history) = realize_gram(
+        let (config, tangent, success, history) = realize_gram(
             &gram, guess, &[],
             1.0e-12, 0.5, 0.9, 1.1, 200, 110
         );
@@ -263,14 +283,111 @@ impl Assembly {
         ));
         console::log_2(&JsValue::from("Steps:"), &JsValue::from(history.scaled_loss.len() - 1));
         console::log_2(&JsValue::from("Loss:"), &JsValue::from(*history.scaled_loss.last().unwrap()));
+        console::log_2(&JsValue::from("Tangent dimension:"), &JsValue::from(tangent.dim()));
         
         if success {
             // 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))
+                    |rep| rep.set_column(0, &config.column(elt.column_index.unwrap()))
                 );
             }
+            
+            // save the tangent space
+            self.tangent.set_silent(tangent);
         }
     }
+    
+    // --- 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) {
+        /* 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) {
+            console::log_1(&JsValue::from("The assembly is rigid"));
+        }
+        
+        // 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
+        });
+        
+        // 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;
+            }
+        }
+        
+        // step each element along the mass shell geodesic that matches its
+        // velocity in the deformation found above
+        /* KLUDGE */
+        // since our test assemblies only include spheres, we assume that every
+        // 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 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
+        self.realize();
+    }
 }

app-proto/src/display.rs

@@ -1,5 +1,5 @@
 use core::array;
-use nalgebra::{DMatrix, Rotation3, Vector3};
+use nalgebra::{DMatrix, DVector, Rotation3, Vector3};
 use sycamore::{prelude::*, motion::create_raf};
 use web_sys::{
     console,
@@ -14,7 +14,7 @@ use web_sys::{
     wasm_bindgen::{JsCast, JsValue}
 };
 
-use crate::{AppState, assembly::ElementKey};
+use crate::{AppState, assembly::{ElementKey, ElementMotion}};
 
 fn compile_shader(
     context: &WebGl2RenderingContext,
@@ -123,6 +123,14 @@ pub fn Display() -> View {
     let zoom_out = create_signal(0.0);
     let turntable = create_signal(false); /* BENCHMARKING */
     
+    // manipulation
+    let translate_neg_x = create_signal(0.0);
+    let translate_pos_x = create_signal(0.0);
+    let translate_neg_y = create_signal(0.0);
+    let translate_pos_y = create_signal(0.0);
+    let translate_neg_z = create_signal(0.0);
+    let translate_pos_z = create_signal(0.0);
+    
     // change listener
     let scene_changed = create_signal(true);
     create_effect(move || {
@@ -141,6 +149,7 @@ pub fn Display() -> View {
     let mut frames_since_last_sample = 0;
     let mean_frame_interval = create_signal(0.0);
     
+    let assembly_for_raf = state.assembly.clone();
     on_mount(move || {
         // timing
         let mut last_time = 0.0;
@@ -153,6 +162,9 @@ pub fn Display() -> View {
         let mut rotation = DMatrix::<f64>::identity(5, 5);
         let mut location_z: f64 = 5.0;
         
+        // manipulation
+        const TRANSLATION_SPEED: f64 = 0.15; // in length units per second
+        
         // display parameters
         const OPACITY: f32 = 0.5; /* SCAFFOLDING */
         const HIGHLIGHT: f32 = 0.2; /* SCAFFOLDING */
@@ -273,6 +285,14 @@ pub fn Display() -> View {
             let zoom_out_val = zoom_out.get();
             let turntable_val = turntable.get(); /* BENCHMARKING */
             
+            // get the manipulation state
+            let translate_neg_x_val = translate_neg_x.get();
+            let translate_pos_x_val = translate_pos_x.get();
+            let translate_neg_y_val = translate_neg_y.get();
+            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();
+            
             // update the assembly's orientation
             let ang_vel = {
                 let pitch = pitch_up_val - pitch_down_val;
@@ -298,6 +318,41 @@ pub fn Display() -> View {
             let zoom = zoom_out_val - zoom_in_val;
             location_z *= (time_step * ZOOM_SPEED * zoom).exp();
             
+            // manipulate the assembly
+            if state.selection.with(|sel| sel.len() == 1) {
+                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;
+                if translate_x != 0.0 || translate_y != 0.0 || translate_z != 0.0 {
+                    let vel_field = {
+                        let u = Vector3::new(translate_x, translate_y, translate_z).normalize();
+                        DMatrix::from_column_slice(5, 5, &[
+                                 0.0,      0.0,      0.0, 0.0, u[0],
+                                 0.0,      0.0,      0.0, 0.0, u[1],
+                                 0.0,      0.0,      0.0, 0.0, u[2],
+                            2.0*u[0], 2.0*u[1], 2.0*u[2], 0.0,  0.0,
+                                 0.0,      0.0,      0.0, 0.0,  0.0
+                        ])
+                    };
+                    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()
+                            }
+                        ]
+                    );
+                    scene_changed.set(true);
+                }
+            }
+            
             if scene_changed.get() {
                 /* INSTRUMENTS */
                 // measure mean frame interval
@@ -416,7 +471,7 @@ pub fn Display() -> View {
         start_animation_loop();
     });
     
-    let set_nav_signal = move |event: KeyboardEvent, value: f64| {
+    let set_nav_signal = move |event: &KeyboardEvent, value: f64| {
         let mut navigating = true;
         let shift = event.shift_key();
         match event.key().as_str() {
@@ -436,6 +491,23 @@ pub fn Display() -> View {
         }
     };
     
+    let set_manip_signal = move |event: &KeyboardEvent, value: f64| {
+        let mut manipulating = true;
+        let shift = event.shift_key();
+        match event.key().as_str() {
+            "d" | "D" => translate_pos_x.set(value),
+            "a" | "A" => translate_neg_x.set(value),
+            "w" | "W" if shift => translate_neg_z.set(value),
+            "s" | "S" if shift => translate_pos_z.set(value),
+            "w" | "W" => translate_pos_y.set(value),
+            "s" | "S" => translate_neg_y.set(value),
+            _ => manipulating = false
+        };
+        if manipulating {
+            event.prevent_default();
+        }
+    };
+    
     view! {
         /* TO DO */
         // switch back to integer-valued parameters when that becomes possible
@@ -447,6 +519,7 @@ pub fn Display() -> View {
             tabindex="0",
             on:keydown=move |event: KeyboardEvent| {
                 if event.key() == "Shift" {
+                    // swap navigation inputs
                     roll_cw.set(yaw_right.get());
                     roll_ccw.set(yaw_left.get());
                     zoom_in.set(pitch_up.get());
@@ -455,16 +528,24 @@ pub fn Display() -> View {
                     yaw_left.set(0.0);
                     pitch_up.set(0.0);
                     pitch_down.set(0.0);
+                    
+                    // swap manipulation inputs
+                    translate_pos_z.set(translate_neg_y.get());
+                    translate_neg_z.set(translate_pos_y.get());
+                    translate_pos_y.set(0.0);
+                    translate_neg_y.set(0.0);
                 } else {
                     if event.key() == "Enter" { /* BENCHMARKING */
                         turntable.set_fn(|turn| !turn);
                         scene_changed.set(true);
                     }
-                    set_nav_signal(event, 1.0);
+                    set_nav_signal(&event, 1.0);
+                    set_manip_signal(&event, 1.0);
                 }
             },
             on:keyup=move |event: KeyboardEvent| {
                 if event.key() == "Shift" {
+                    // swap navigation inputs
                     yaw_right.set(roll_cw.get());
                     yaw_left.set(roll_ccw.get());
                     pitch_up.set(zoom_in.get());
@@ -473,8 +554,15 @@ pub fn Display() -> View {
                     roll_ccw.set(0.0);
                     zoom_in.set(0.0);
                     zoom_out.set(0.0);
+                    
+                    // swap manipulation inputs
+                    translate_pos_y.set(translate_neg_z.get());
+                    translate_neg_y.set(translate_pos_z.get());
+                    translate_pos_z.set(0.0);
+                    translate_neg_z.set(0.0);
                 } else {
-                    set_nav_signal(event, 0.0);
+                    set_nav_signal(&event, 0.0);
+                    set_manip_signal(&event, 0.0);
                 }
             },
             on:blur=move |_| {

app-proto/src/engine.rs

@@ -1,5 +1,5 @@
 use lazy_static::lazy_static;
-use nalgebra::{Const, DMatrix, DVector, Dyn};
+use nalgebra::{Const, DMatrix, DVector, DVectorView, Dyn, SymmetricEigen};
 use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */
 
 // --- elements ---
@@ -85,6 +85,75 @@ impl PartialMatrix {
     }
 }
 
+// --- configuration subspaces ---
+
+#[derive(Clone)]
+pub struct ConfigSubspace {
+    assembly_dim: usize,
+    basis: Vec<DMatrix<f64>>
+}
+
+impl ConfigSubspace {
+    pub fn zero(assembly_dim: usize) -> ConfigSubspace {
+        ConfigSubspace {
+            assembly_dim: assembly_dim,
+            basis: 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 {
+        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))
+                )
+            )
+        );
+        
+        /* DEBUG */
+        // print the eigenvalues
+        #[cfg(all(target_family = "wasm", target_os = "unknown"))]
+        console::log_1(&JsValue::from(
+            format!("Eigenvalues used to find kernel:{}", eig.eigenvalues)
+        ));
+        
+        ConfigSubspace {
+            assembly_dim: assembly_dim,
+            basis: basis.collect()
+        }
+    }
+    
+    pub fn dim(&self) -> usize {
+        self.basis.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`
+    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
+            ).sum()
+        }
+    }
+}
+
 // --- descent history ---
 
 pub struct DescentHistory {
@@ -181,7 +250,7 @@ pub fn realize_gram(
     reg_scale: f64,
     max_descent_steps: i32,
     max_backoff_steps: i32
-) -> (DMatrix<f64>, bool, DescentHistory) {
+) -> (DMatrix<f64>, ConfigSubspace, bool, DescentHistory) {
     // start the descent history
     let mut history = DescentHistory::new();
     
@@ -201,12 +270,8 @@ pub fn realize_gram(
     
     // use Newton's method with backtracking and gradient descent backup
     let mut state = SearchState::from_config(gram, guess);
+    let mut hess = DMatrix::zeros(element_dim, assembly_dim);
     for _ in 0..max_descent_steps {
-        // stop if the loss is tolerably low
-        history.config.push(state.config.clone());
-        history.scaled_loss.push(state.loss / scale_adjustment);
-        if state.loss < tol { break; }
-        
         // find the negative gradient of the loss function
         let neg_grad = 4.0 * &*Q * &state.config * &state.err_proj;
         let mut neg_grad_stacked = neg_grad.clone().reshape_generic(Dyn(total_dim), Const::<1>);
@@ -229,7 +294,7 @@ pub fn realize_gram(
                 hess_cols.push(deriv_grad.reshape_generic(Dyn(total_dim), Const::<1>));
             }
         }
-        let mut hess = DMatrix::from_columns(hess_cols.as_slice());
+        hess = DMatrix::from_columns(hess_cols.as_slice());
         
         // regularize the Hessian
         let min_eigval = hess.symmetric_eigenvalues().min();
@@ -249,6 +314,11 @@ pub fn realize_gram(
             hess[(k, k)] = 1.0;
         }
         
+        // stop if the loss is tolerably low
+        history.config.push(state.config.clone());
+        history.scaled_loss.push(state.loss / scale_adjustment);
+        if state.loss < tol { break; }
+        
         // compute the Newton step
         /*
           we need to either handle or eliminate the case where the minimum
@@ -256,7 +326,7 @@ pub fn realize_gram(
           singular. right now, this causes the Cholesky decomposition to return
           `None`, leading to a panic when we unrap
         */
-        let base_step_stacked = hess.cholesky().unwrap().solve(&neg_grad_stacked);
+        let base_step_stacked = hess.clone().cholesky().unwrap().solve(&neg_grad_stacked);
         let base_step = base_step_stacked.reshape_generic(Dyn(element_dim), Dyn(assembly_dim));
         history.base_step.push(base_step.clone());
         
@@ -269,10 +339,16 @@ pub fn realize_gram(
                 state = better_state;
                 history.backoff_steps.push(backoff_steps);
             },
-            None => return (state.config, false, history)
+            None => return (state.config, ConfigSubspace::zero(assembly_dim), false, history)
         };
     }
-    (state.config, state.loss < tol, history)
+    let success = state.loss < tol;
+    let tangent = if success {
+        ConfigSubspace::symmetric_kernel(hess, assembly_dim)
+    } else {
+        ConfigSubspace::zero(assembly_dim)
+    };
+    (state.config, tangent, success, history)
 }
 
 // --- tests ---
@@ -291,7 +367,7 @@ pub mod irisawa {
     
     use super::*;
     
-    pub fn realize_irisawa_hexlet(scaled_tol: f64) -> (DMatrix<f64>, bool, DescentHistory) {
+    pub fn realize_irisawa_hexlet(scaled_tol: f64) -> (DMatrix<f64>, ConfigSubspace, bool, DescentHistory) {
         let gram = {
             let mut gram_to_be = PartialMatrix::new();
             for s in 0..9 {
@@ -399,7 +475,7 @@ mod tests {
     fn irisawa_hexlet_test() {
         // solve Irisawa's problem
         const SCALED_TOL: f64 = 1.0e-12;
-        let (config, _, _) = realize_irisawa_hexlet(SCALED_TOL);
+        let (config, _, _, _) = realize_irisawa_hexlet(SCALED_TOL);
         
         // check against Irisawa's solution
         let entry_tol = SCALED_TOL.sqrt();
@@ -409,6 +485,61 @@ mod tests {
         }
     }
     
+    #[test]
+    fn tangent_test() {
+        const SCALED_TOL: f64 = 1.0e-12;
+        const ELEMENT_DIM: usize = 5;
+        const ASSEMBLY_DIM: usize = 3;
+        let gram = {
+            let mut gram_to_be = PartialMatrix::new();
+            for j in 0..3 {
+                for k in j..3 {
+                    gram_to_be.push_sym(j, k, if j == k { 1.0 } else { -1.0 });
+                }
+            }
+            gram_to_be
+        };
+        let guess = DMatrix::from_columns(&[
+            sphere(0.0, 0.0, 0.0, -2.0),
+            sphere(0.0, 0.0, 1.0, 1.0),
+            sphere(0.0, 0.0, -1.0, 1.0)
+        ]);
+        let frozen: [_; 5] = std::array::from_fn(|k| (k, 0));
+        let (config, tangent, success, history) = realize_gram(
+            &gram, guess.clone(), &frozen,
+            SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
+        );
+        assert_eq!(config, guess);
+        assert_eq!(success, true);
+        assert_eq!(history.scaled_loss.len(), 1);
+        
+        // confirm that the tangent space has dimension five or less
+        let ConfigSubspace(ref tangent_basis) = tangent;
+        assert_eq!(tangent_basis.len(), 5);
+        
+        // confirm that the tangent space contains all the motions we expect it
+        // to. since we've already bounded the dimension of the tangent space,
+        // this confirms that the tangent space is what we expect it to be
+        let tangent_motions = vec![
+            basis_matrix((0, 1), ELEMENT_DIM, ASSEMBLY_DIM),
+            basis_matrix((1, 1), ELEMENT_DIM, ASSEMBLY_DIM),
+            basis_matrix((0, 2), ELEMENT_DIM, ASSEMBLY_DIM),
+            basis_matrix((1, 2), ELEMENT_DIM, ASSEMBLY_DIM),
+            DMatrix::<f64>::from_column_slice(ELEMENT_DIM, 3, &[
+                0.0, 0.0, 0.0, 0.0, 0.0,
+                0.0, 0.0, -1.0, -0.25, -1.0,
+                0.0, 0.0, -1.0, 0.25, 1.0
+            ])
+        ];
+        let tol_sq = ((ELEMENT_DIM * ASSEMBLY_DIM) as f64) * SCALED_TOL * SCALED_TOL;
+        for motion in tangent_motions {
+            let motion_proj: DMatrix<_> = motion.column_iter().enumerate().map(
+                |(k, v)| tangent.proj(&v, k)
+            ).sum();
+            assert!((motion - motion_proj).norm_squared() < tol_sq);
+        }
+    }
+    
     // at the frozen indices, the optimization steps should have exact zeros,
     // and the realized configuration should match the initial guess
     #[test]
@@ -428,7 +559,7 @@ mod tests {
         ]);
         let frozen = [(3, 0), (3, 1)];
         println!();
-        let (config, success, history) = realize_gram(
+        let (config, _, success, history) = realize_gram(
             &gram, guess.clone(), &frozen,
             1.0e-12, 0.5, 0.9, 1.1, 200, 110
         );

app-proto/src/main.rs

@@ -46,6 +46,10 @@ impl AppState {
 }
 
 fn main() {
+    // set the console error panic hook
+    #[cfg(feature = "console_error_panic_hook")]
+    console_error_panic_hook::set_once();
+    
     sycamore::render(|| {
         provide_context(AppState::new());