Dispatch normalization routines correctly (#87)

Addresses issue #86 by correctly dispatching the routine used to normalize spheres during nudging. Adds a test that would have detected the issue.

Since the tests aren't built for WebAssembly, we have to replace `console::log` with `console_log!` in all of the functions used by `assembly::curvature_drift_test`. We'll eventually want to do this replacement everywhere.

Co-authored-by: Aaron Fenyes <aaron.fenyes@fareycircles.ooo>
Reviewed-on: #87
Co-authored-by: Vectornaut <vectornaut@nobody@nowhere.net>
Co-committed-by: Vectornaut <vectornaut@nobody@nowhere.net>

app-proto/Cargo.toml

@@ -46,6 +46,13 @@ features = [
 dyna3 = { path = ".", default-features = false, features = ["dev"] }
 wasm-bindgen-test = "0.3.34"
 
+# turn off spurious warnings about the custom config that Sycamore uses
+#
+#   https://sycamore.dev/book/troubleshooting#unexpected-cfg-condition-name--sycamore-force-ssr
+#
+[lints.rust]
+unexpected_cfgs = { level = "warn", check-cfg = ["cfg(sycamore_force_ssr)"] }
+
 [profile.release]
 opt-level = "s" # optimize for small code size
 debug = true # include debug symbols

app-proto/src/assembly.rs

@@ -1,6 +1,5 @@
 use nalgebra::{DMatrix, DVector, DVectorView};
 use std::{
-    any::{Any, TypeId},
     cell::Cell,
     collections::{BTreeMap, BTreeSet},
     cmp::Ordering,
@@ -20,6 +19,8 @@ use crate::{
         change_half_curvature,
         local_unif_to_std,
         point,
+        project_point_to_normalized,
+        project_sphere_to_normalized,
         realize_gram,
         sphere,
         ConfigSubspace,
@@ -107,6 +108,10 @@ pub trait Element: Serial + ProblemPoser + DisplayItem {
     // element is responsible for keeping this set up to date
     fn regulators(&self) -> Signal<BTreeSet<Rc<dyn Regulator>>>;
     
+    // project a representation vector for this kind of element onto its
+    // normalization variety
+    fn project_to_normalized(&self, rep: &mut DVector<f64>);
+    
     // the configuration matrix column index that was assigned to the element
     // last time the assembly was realized, or `None` if the element has never
     // been through a realization
@@ -221,6 +226,10 @@ impl Element for Sphere {
         self.regulators
     }
     
+    fn project_to_normalized(&self, rep: &mut DVector<f64>) {
+        project_sphere_to_normalized(rep);
+    }
+    
     fn column_index(&self) -> Option<usize> {
         self.column_index.get()
     }
@@ -313,6 +322,10 @@ impl Element for Point {
         self.regulators
     }
     
+    fn project_to_normalized(&self, rep: &mut DVector<f64>) {
+        project_point_to_normalized(rep);
+    }
+    
     fn column_index(&self) -> Option<usize> {
         self.column_index.get()
     }
@@ -611,9 +624,7 @@ impl Assembly {
         create_effect(move || {
             /* DEBUG */
             // log the regulator update
-            console::log_1(&JsValue::from(
+            console_log!("Updated regulator with subjects {:?}", regulator.subjects());
-                format!("Updated regulator with subjects {:?}", regulator.subjects())
-            ));
             
             if regulator.try_activate() {
                 self_for_effect.realize();
@@ -622,10 +633,10 @@ impl Assembly {
         
         /* DEBUG */
         // print an updated list of regulators
-        console::log_1(&JsValue::from("Regulators:"));
+        console_log!("Regulators:");
         self.regulators.with_untracked(|regs| {
             for reg in regs.into_iter() {
-                console::log_1(&JsValue::from(format!(
+                console_log!(
                     "  {:?}: {}",
                     reg.subjects(),
                     reg.set_point().with_untracked(
@@ -638,7 +649,7 @@ impl Assembly {
                             }
                         }
                     )
-                )));
+                );
             }
         });
     }
@@ -669,19 +680,11 @@ impl Assembly {
         
         /* DEBUG */
         // log the Gram matrix
-        console::log_1(&JsValue::from("Gram matrix:"));
+        console_log!("Gram matrix:\n{}", problem.gram);
-        problem.gram.log_to_console();
         
         /* DEBUG */
         // log the initial configuration matrix
-        console::log_1(&JsValue::from("Old configuration:"));
+        console_log!("Old configuration:{:>8.3}", problem.guess);
-        for j in 0..problem.guess.nrows() {
-            let mut row_str = String::new();
-            for k in 0..problem.guess.ncols() {
-                row_str.push_str(format!(" {:>8.3}", problem.guess[(j, k)]).as_str());
-            }
-            console::log_1(&JsValue::from(row_str));
-        }
         
         // look for a configuration with the given Gram matrix
         let (config, tangent, success, history) = realize_gram(
@@ -690,16 +693,14 @@ impl Assembly {
         
         /* DEBUG */
         // report the outcome of the search
-        console::log_1(&JsValue::from(
+        if success {
-            if success {
+            console_log!("Target accuracy achieved!")
-                "Target accuracy achieved!"
+        } else {
-            } else {
+            console_log!("Failed to reach target accuracy")
-                "Failed to reach target accuracy"
+        }
-            }
+        console_log!("Steps: {}", history.scaled_loss.len() - 1);
-        ));
+        console_log!("Loss: {}", *history.scaled_loss.last().unwrap());
-        console::log_2(&JsValue::from("Steps:"), &JsValue::from(history.scaled_loss.len() - 1));
+        console_log!("Tangent dimension: {}", tangent.dim());
-        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
@@ -782,29 +783,17 @@ impl Assembly {
         
         // step the assembly along the deformation. this changes the elements'
         // normalizations, so we restore those afterward
-        /* KLUDGE */
-        // for now, we only restore the normalizations of spheres
         for elt in self.elements.get_clone_untracked() {
             elt.representation().update_silent(|rep| {
                 match elt.column_index() {
                     Some(column_index) => {
-                        // step the assembly along the deformation
+                        // step the element along the deformation and then
+                        // restore its normalization
                         *rep += motion_proj.column(column_index);
-                        
+                        elt.project_to_normalized(rep);
-                        if elt.type_id() == TypeId::of::<Sphere>() {
-                            // restore normalization by contracting toward the
-                            // last coordinate axis
-                            let q_sp = rep.fixed_rows::<3>(0).norm_squared();
-                            let half_q_lt = -2.0 * rep[3] * rep[4];
-                            let half_q_lt_sq = half_q_lt * half_q_lt;
-                            let scaling = half_q_lt + (q_sp + half_q_lt_sq).sqrt();
-                            rep.fixed_rows_mut::<4>(0).scale_mut(1.0 / scaling);
-                        }
                     },
                     None => {
-                        console::log_1(&JsValue::from(
+                        console_log!("No velocity to unpack for fresh element \"{}\"", elt.id())
-                            format!("No velocity to unpack for fresh element \"{}\"", elt.id())
-                        ))
                     }
                 };
             });
@@ -821,6 +810,8 @@ impl Assembly {
 mod tests {
     use super::*;
     
+    use crate::engine;
+    
     #[test]
     #[should_panic(expected = "Sphere \"sphere\" should be indexed before writing problem data")]
     fn unindexed_element_test() {
@@ -846,4 +837,50 @@ mod tests {
             }.pose(&mut ConstraintProblem::new(2));
         });
     }
+    
+    #[test]
+    fn curvature_drift_test() {
+        const INITIAL_RADIUS: f64 = 0.25;
+        let _ = create_root(|| {
+            // set up an assembly containing a single sphere centered at the
+            // origin
+            let assembly = Assembly::new();
+            let sphere_id = "sphere0";
+            let _ = assembly.try_insert_element(
+                // we create the sphere by hand for two reasons: to choose the
+                // curvature (which can affect drift rate) and to make the test
+                // independent of `Sphere::default`
+                Sphere::new(
+                    String::from(sphere_id),
+                    String::from("Sphere 0"),
+                    [0.75_f32, 0.75_f32, 0.75_f32],
+                    engine::sphere(0.0, 0.0, 0.0, INITIAL_RADIUS)
+                )
+            );
+            
+            // nudge the sphere repeatedly along the `z` axis
+            const STEP_SIZE: f64 = 0.0025;
+            const STEP_CNT: usize = 400;
+            let sphere = assembly.elements_by_id.with(|elts_by_id| elts_by_id[sphere_id].clone());
+            let velocity = DVector::from_column_slice(&[0.0, 0.0, STEP_SIZE, 0.0]);
+            for _ in 0..STEP_CNT {
+                assembly.deform(
+                    vec![
+                        ElementMotion {
+                            element: sphere.clone(),
+                            velocity: velocity.as_view()
+                        }
+                    ]
+                );
+            }
+            
+            // check how much the sphere's curvature has drifted
+            const INITIAL_HALF_CURV: f64 = 0.5 / INITIAL_RADIUS;
+            const DRIFT_TOL: f64 = 0.015;
+            let final_half_curv = sphere.representation().with_untracked(
+                |rep| rep[Sphere::CURVATURE_COMPONENT]
+            );
+            assert!((final_half_curv / INITIAL_HALF_CURV - 1.0).abs() < DRIFT_TOL);
+        });
+    }
 }

app-proto/src/engine.rs

@@ -1,5 +1,6 @@
 use lazy_static::lazy_static;
 use nalgebra::{Const, DMatrix, DVector, DVectorView, Dyn, SymmetricEigen};
+use std::fmt::{Display, Error, Formatter};
 use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */
 
 // --- elements ---
@@ -34,6 +35,21 @@ pub fn sphere_with_offset(dir_x: f64, dir_y: f64, dir_z: f64, off: f64, curv: f6
     ])
 }
 
+// project a sphere's representation vector to the normalization variety by
+// contracting toward the last coordinate axis
+pub fn project_sphere_to_normalized(rep: &mut DVector<f64>) {
+    let q_sp = rep.fixed_rows::<3>(0).norm_squared();
+    let half_q_lt = -2.0 * rep[3] * rep[4];
+    let half_q_lt_sq = half_q_lt * half_q_lt;
+    let scaling = half_q_lt + (q_sp + half_q_lt_sq).sqrt();
+    rep.fixed_rows_mut::<4>(0).scale_mut(1.0 / scaling);
+}
+
+// normalize a point's representation vector by scaling
+pub fn project_point_to_normalized(rep: &mut DVector<f64>) {
+    rep.scale_mut(0.5 / rep[3]);
+}
+
 // given a sphere's representation vector, change the sphere's half-curvature to
 // `half-curv` and then restore normalization by contracting the representation
 // vector toward the curvature axis
@@ -94,15 +110,6 @@ impl PartialMatrix {
         }
     }
     
-    /* DEBUG */
-    pub fn log_to_console(&self) {
-        for &MatrixEntry { index: (row, col), value } in self {
-            console::log_1(&JsValue::from(
-                format!("  {} {} {}", row, col, value)
-            ));
-        }
-    }
-    
     fn freeze(&self, a: &DMatrix<f64>) -> DMatrix<f64> {
         let mut result = a.clone();
         for &MatrixEntry { index, value } in self {
@@ -128,6 +135,15 @@ impl PartialMatrix {
     }
 }
 
+impl Display for PartialMatrix {
+    fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), Error> {
+        for &MatrixEntry { index: (row, col), value } in self {
+            writeln!(f, "  {row} {col} {value}")?;
+        }
+        Ok(())
+    }
+}
+
 impl IntoIterator for PartialMatrix {
     type Item = MatrixEntry;
     type IntoIter = std::vec::IntoIter<Self::Item>;