diff --git a/app-proto/examples/irisawa-hexlet.rs b/app-proto/examples/irisawa-hexlet.rs index 2bc94c0..fc14f91 100644 --- a/app-proto/examples/irisawa-hexlet.rs +++ b/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!"); diff --git a/app-proto/examples/point-on-sphere.rs b/app-proto/examples/point-on-sphere.rs index 13040e5..0e4765a 100644 --- a/app-proto/examples/point-on-sphere.rs +++ b/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 ); diff --git a/app-proto/examples/three-spheres.rs b/app-proto/examples/three-spheres.rs index 19acfd1..d348b18 100644 --- a/app-proto/examples/three-spheres.rs +++ b/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 ); diff --git a/app-proto/src/assembly.rs b/app-proto/src/assembly.rs index 24f9b93..fb5bbf7 100644 --- a/app-proto/src/assembly.rs +++ b/app-proto/src/assembly.rs @@ -1,11 +1,11 @@ -use nalgebra::{DMatrix, DVector, DVectorView, Vector3}; +use nalgebra::{DMatrix, DVector, 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, ConfigSubspace, PartialMatrix, Q}; +use crate::engine::{realize_gram, PartialMatrix}; // the types of the keys we use to access an assembly's elements and constraints pub type ElementKey = usize; @@ -109,6 +109,7 @@ impl Element { } } } + #[derive(Clone)] pub struct Constraint { @@ -126,9 +127,6 @@ pub struct Assembly { pub elements: Signal>, pub constraints: Signal>, - // solution variety tangent space - pub tangent: Signal, - // indexing pub elements_by_id: Signal> } @@ -138,7 +136,6 @@ impl Assembly { Assembly { elements: create_signal(Slab::new()), constraints: create_signal(Slab::new()), - tangent: create_signal(ConfigSubspace::zero()), elements_by_id: create_signal(FxHashMap::default()) } } @@ -152,13 +149,6 @@ 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 { @@ -257,7 +247,7 @@ impl Assembly { } // look for a configuration with the given Gram matrix - let (config, tangent, success, history) = realize_gram( + let (config, success, history) = realize_gram( &gram, guess, &[], 1.0e-12, 0.5, 0.9, 1.1, 200, 110 ); @@ -273,7 +263,6 @@ 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 @@ -282,50 +271,6 @@ impl Assembly { |rep| rep.set_column(0, &config.column(elt.column_index)) ); } - - // save the tangent space - self.tangent.set_silent(tangent); } } - - // --- deformation --- - - pub fn deform(&self, element_motions: Vec<(ElementKey, DVectorView)>) { - /* KLUDGE */ - // 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; - } - - 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 - 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 - // velocity in the deformation found above - /* KLUDGE */ - // since our test assemblies only involve 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| { - 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)); - }); - } - - // bring the configuration back onto the solution variety - self.realize(); - } } \ No newline at end of file diff --git a/app-proto/src/display.rs b/app-proto/src/display.rs index 184b29e..c39e575 100644 --- a/app-proto/src/display.rs +++ b/app-proto/src/display.rs @@ -1,5 +1,5 @@ use core::array; -use nalgebra::{DMatrix, DVector, Rotation3, Vector3}; +use nalgebra::{DMatrix, Rotation3, Vector3}; use sycamore::{prelude::*, motion::create_raf}; use web_sys::{ console, @@ -123,10 +123,6 @@ 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); - // change listener let scene_changed = create_signal(true); create_effect(move || { @@ -145,7 +141,6 @@ 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; @@ -158,9 +153,6 @@ pub fn Display() -> View { let mut rotation = DMatrix::::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 */ @@ -281,10 +273,6 @@ 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(); - // update the assembly's orientation let ang_vel = { let pitch = pitch_up_val - pitch_down_val; @@ -310,30 +298,6 @@ 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 vel_field_x = DMatrix::from_column_slice(5, 5, &[ - 0.0, 0.0, 0.0, 0.0, 1.0, - 0.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 0.0, - 2.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 0.0 - ]); - let translate_x = translate_pos_x_val - translate_neg_x_val; - if translate_x != 0.0 { - let vel = translate_x * vel_field_x * rep; - let elt_motion: DVector = time_step * TRANSLATION_SPEED * vel; - assembly_for_raf.deform(vec![(sel, elt_motion.as_view())]); - scene_changed.set(true); - } - } - if scene_changed.get() { /* INSTRUMENTS */ // measure mean frame interval @@ -452,7 +416,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() { @@ -472,18 +436,6 @@ pub fn Display() -> View { } }; - let set_manip_signal = move |event: &KeyboardEvent, value: f64| { - let mut manipulating = true; - match event.key().as_str() { - "d" => translate_pos_x.set(value), - "a" => translate_neg_x.set(value), - _ => manipulating = false - }; - if manipulating { - event.prevent_default(); - } - }; - view! { /* TO DO */ // switch back to integer-valued parameters when that becomes possible @@ -508,8 +460,7 @@ pub fn Display() -> View { turntable.set_fn(|turn| !turn); scene_changed.set(true); } - set_nav_signal(&event, 1.0); - set_manip_signal(&event, 1.0); + set_nav_signal(event, 1.0); } }, on:keyup=move |event: KeyboardEvent| { @@ -523,8 +474,7 @@ pub fn Display() -> View { zoom_in.set(0.0); zoom_out.set(0.0); } else { - set_nav_signal(&event, 0.0); - set_manip_signal(&event, 0.0); + set_nav_signal(event, 0.0); } }, on:blur=move |_| { diff --git a/app-proto/src/engine.rs b/app-proto/src/engine.rs index a48a1af..285a9c6 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -1,5 +1,5 @@ use lazy_static::lazy_static; -use nalgebra::{Const, DMatrix, DVector, DVectorView, Dyn, SymmetricEigen}; +use nalgebra::{Const, DMatrix, DVector, Dyn}; use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */ // --- elements --- @@ -85,56 +85,6 @@ impl PartialMatrix { } } -// --- configuration subspaces --- - -#[derive(Clone)] -pub struct ConfigSubspace(Vec>); - -impl ConfigSubspace { - pub fn zero() -> ConfigSubspace { - ConfigSubspace(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, 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::>::into( - v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim)) - ) - ) - ); - console::log_1(&JsValue::from( - format!("Hessian eigenvalues: {}", eig.eigenvalues) - )); /* DEBUG */ - ConfigSubspace(basis.collect()) - } - - pub fn dim(&self) -> usize { - 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, column_index: usize) -> DMatrix { - let ConfigSubspace(basis) = self; - basis.into_iter().map( - |b| b.column(column_index).dot(&v) * b - ).sum() - } -} - // --- descent history --- pub struct DescentHistory { @@ -231,7 +181,7 @@ pub fn realize_gram( reg_scale: f64, max_descent_steps: i32, max_backoff_steps: i32 -) -> (DMatrix, ConfigSubspace, bool, DescentHistory) { +) -> (DMatrix, bool, DescentHistory) { // start the descent history let mut history = DescentHistory::new(); @@ -251,8 +201,12 @@ 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>); @@ -275,7 +229,7 @@ pub fn realize_gram( hess_cols.push(deriv_grad.reshape_generic(Dyn(total_dim), Const::<1>)); } } - hess = DMatrix::from_columns(hess_cols.as_slice()); + let mut hess = DMatrix::from_columns(hess_cols.as_slice()); // regularize the Hessian let min_eigval = hess.symmetric_eigenvalues().min(); @@ -295,11 +249,6 @@ 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 @@ -307,7 +256,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.clone().cholesky().unwrap().solve(&neg_grad_stacked); + let base_step_stacked = hess.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()); @@ -320,16 +269,10 @@ pub fn realize_gram( state = better_state; history.backoff_steps.push(backoff_steps); }, - None => return (state.config, ConfigSubspace::zero(), false, history) + None => return (state.config, false, history) }; } - let success = state.loss < tol; - let tangent = if success { - ConfigSubspace::symmetric_kernel(hess, assembly_dim) - } else { - ConfigSubspace::zero() - }; - (state.config, tangent, success, history) + (state.config, state.loss < tol, history) } // --- tests --- @@ -348,7 +291,7 @@ pub mod irisawa { use super::*; - pub fn realize_irisawa_hexlet(scaled_tol: f64) -> (DMatrix, ConfigSubspace, bool, DescentHistory) { + pub fn realize_irisawa_hexlet(scaled_tol: f64) -> (DMatrix, bool, DescentHistory) { let gram = { let mut gram_to_be = PartialMatrix::new(); for s in 0..9 { @@ -456,7 +399,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(); @@ -466,61 +409,6 @@ 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::::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] @@ -540,7 +428,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 ); diff --git a/app-proto/src/main.rs b/app-proto/src/main.rs index f961504..8a012d3 100644 --- a/app-proto/src/main.rs +++ b/app-proto/src/main.rs @@ -46,10 +46,6 @@ 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());