diff --git a/app-proto/examples/point-on-sphere.rs b/app-proto/examples/point-on-sphere.rs index 13040e5..880d7b0 100644 --- a/app-proto/examples/point-on-sphere.rs +++ b/app-proto/examples/point-on-sphere.rs @@ -1,26 +1,19 @@ -use nalgebra::DMatrix; - -use dyna3::engine::{Q, point, realize_gram, sphere, PartialMatrix}; +use dyna3::engine::{Q, point, realize_gram, sphere, ConstraintProblem}; fn main() { - let gram = { - let mut gram_to_be = PartialMatrix::new(); - for j in 0..2 { - for k in j..2 { - gram_to_be.push_sym(j, k, if (j, k) == (1, 1) { 1.0 } else { 0.0 }); - } - } - gram_to_be - }; - let guess = DMatrix::from_columns(&[ + let mut problem = ConstraintProblem::from_guess(&[ point(0.0, 0.0, 2.0), sphere(0.0, 0.0, 0.0, 1.0) ]); - let frozen = [(3, 0)]; + for j in 0..2 { + for k in j..2 { + problem.gram.push_sym(j, k, if (j, k) == (1, 1) { 1.0 } else { 0.0 }); + } + } + problem.frozen.push(3, 0, problem.guess[(3, 0)]); println!(); let (config, _, success, history) = realize_gram( - &gram, guess, &frozen, - 1.0e-12, 0.5, 0.9, 1.1, 200, 110 + &problem, 1.0e-12, 0.5, 0.9, 1.1, 200, 110 ); print!("\nCompleted Gram matrix:{}", config.tr_mul(&*Q) * &config); print!("Configuration:{}", config); diff --git a/app-proto/examples/three-spheres.rs b/app-proto/examples/three-spheres.rs index 19acfd1..3f3cc44 100644 --- a/app-proto/examples/three-spheres.rs +++ b/app-proto/examples/three-spheres.rs @@ -1,29 +1,22 @@ -use nalgebra::DMatrix; - -use dyna3::engine::{Q, realize_gram, sphere, PartialMatrix}; +use dyna3::engine::{Q, realize_gram, sphere, ConstraintProblem}; fn main() { - 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 = { + let mut problem = ConstraintProblem::from_guess({ let a: f64 = 0.75_f64.sqrt(); - DMatrix::from_columns(&[ + &[ sphere(1.0, 0.0, 0.0, 1.0), sphere(-0.5, a, 0.0, 1.0), sphere(-0.5, -a, 0.0, 1.0) - ]) - }; + ] + }); + for j in 0..3 { + for k in j..3 { + problem.gram.push_sym(j, k, if j == k { 1.0 } else { -1.0 }); + } + } println!(); let (config, _, success, history) = realize_gram( - &gram, guess, &[], - 1.0e-12, 0.5, 0.9, 1.1, 200, 110 + &problem, 1.0e-12, 0.5, 0.9, 1.1, 200, 110 ); print!("\nCompleted Gram matrix:{}", config.tr_mul(&*Q) * &config); if success { diff --git a/app-proto/src/add_remove.rs b/app-proto/src/add_remove.rs index 5fed411..14fcd41 100644 --- a/app-proto/src/add_remove.rs +++ b/app-proto/src/add_remove.rs @@ -4,14 +4,14 @@ use web_sys::{console, wasm_bindgen::JsValue}; use crate::{ engine, AppState, - assembly::{Assembly, Element} + assembly::{Assembly, Element, InversiveDistanceRegulator} }; /* DEBUG */ // load an example assembly for testing. this code will be removed once we've // built a more formal test assembly system fn load_gen_assemb(assembly: &Assembly) { - let _ = assembly.try_insert_element( + let _ = assembly.try_insert_sphere( Element::new( String::from("gemini_a"), String::from("Castor"), @@ -19,7 +19,7 @@ fn load_gen_assemb(assembly: &Assembly) { engine::sphere(0.5, 0.5, 0.0, 1.0) ) ); - let _ = assembly.try_insert_element( + let _ = assembly.try_insert_sphere( Element::new( String::from("gemini_b"), String::from("Pollux"), @@ -27,7 +27,7 @@ fn load_gen_assemb(assembly: &Assembly) { engine::sphere(-0.5, -0.5, 0.0, 1.0) ) ); - let _ = assembly.try_insert_element( + let _ = assembly.try_insert_sphere( Element::new( String::from("ursa_major"), String::from("Ursa major"), @@ -35,7 +35,7 @@ fn load_gen_assemb(assembly: &Assembly) { engine::sphere(-0.5, 0.5, 0.0, 0.75) ) ); - let _ = assembly.try_insert_element( + let _ = assembly.try_insert_sphere( Element::new( String::from("ursa_minor"), String::from("Ursa minor"), @@ -43,7 +43,7 @@ fn load_gen_assemb(assembly: &Assembly) { engine::sphere(0.5, -0.5, 0.0, 0.5) ) ); - let _ = assembly.try_insert_element( + let _ = assembly.try_insert_sphere( Element::new( String::from("moon_deimos"), String::from("Deimos"), @@ -51,7 +51,7 @@ fn load_gen_assemb(assembly: &Assembly) { engine::sphere(0.0, 0.15, 1.0, 0.25) ) ); - let _ = assembly.try_insert_element( + let _ = assembly.try_insert_sphere( Element::new( String::from("moon_phobos"), String::from("Phobos"), @@ -66,7 +66,7 @@ fn load_gen_assemb(assembly: &Assembly) { // built a more formal test assembly system fn load_low_curv_assemb(assembly: &Assembly) { let a = 0.75_f64.sqrt(); - let _ = assembly.try_insert_element( + let _ = assembly.try_insert_sphere( Element::new( "central".to_string(), "Central".to_string(), @@ -74,7 +74,7 @@ fn load_low_curv_assemb(assembly: &Assembly) { engine::sphere(0.0, 0.0, 0.0, 1.0) ) ); - let _ = assembly.try_insert_element( + let _ = assembly.try_insert_sphere( Element::new( "assemb_plane".to_string(), "Assembly plane".to_string(), @@ -82,7 +82,7 @@ fn load_low_curv_assemb(assembly: &Assembly) { engine::sphere_with_offset(0.0, 0.0, 1.0, 0.0, 0.0) ) ); - let _ = assembly.try_insert_element( + let _ = assembly.try_insert_sphere( Element::new( "side1".to_string(), "Side 1".to_string(), @@ -90,7 +90,7 @@ fn load_low_curv_assemb(assembly: &Assembly) { engine::sphere_with_offset(1.0, 0.0, 0.0, 1.0, 0.0) ) ); - let _ = assembly.try_insert_element( + let _ = assembly.try_insert_sphere( Element::new( "side2".to_string(), "Side 2".to_string(), @@ -98,7 +98,7 @@ fn load_low_curv_assemb(assembly: &Assembly) { engine::sphere_with_offset(-0.5, a, 0.0, 1.0, 0.0) ) ); - let _ = assembly.try_insert_element( + let _ = assembly.try_insert_sphere( Element::new( "side3".to_string(), "Side 3".to_string(), @@ -106,7 +106,7 @@ fn load_low_curv_assemb(assembly: &Assembly) { engine::sphere_with_offset(-0.5, -a, 0.0, 1.0, 0.0) ) ); - let _ = assembly.try_insert_element( + let _ = assembly.try_insert_sphere( Element::new( "corner1".to_string(), "Corner 1".to_string(), @@ -114,7 +114,7 @@ fn load_low_curv_assemb(assembly: &Assembly) { engine::sphere(-4.0/3.0, 0.0, 0.0, 1.0/3.0) ) ); - let _ = assembly.try_insert_element( + let _ = assembly.try_insert_sphere( Element::new( "corner2".to_string(), "Corner 2".to_string(), @@ -122,7 +122,7 @@ fn load_low_curv_assemb(assembly: &Assembly) { engine::sphere(2.0/3.0, -4.0/3.0 * a, 0.0, 1.0/3.0) ) ); - let _ = assembly.try_insert_element( + let _ = assembly.try_insert_sphere( Element::new( String::from("corner3"), String::from("Corner 3"), @@ -148,6 +148,7 @@ pub fn AddRemove() -> View { let assembly = &state.assembly; // clear state + assembly.regulators.update(|regs| regs.clear()); assembly.elements.update(|elts| elts.clear()); assembly.elements_by_id.update(|elts_by_id| elts_by_id.clear()); state.selection.update(|sel| sel.clear()); @@ -166,18 +167,7 @@ pub fn AddRemove() -> View { button( on:click=|_| { let state = use_context::(); - state.assembly.insert_new_element(); - - /* DEBUG */ - // print updated list of elements by identifier - console::log_1(&JsValue::from("elements by identifier:")); - for (id, key) in state.assembly.elements_by_id.get_clone().iter() { - console::log_3( - &JsValue::from(" "), - &JsValue::from(id), - &JsValue::from(*key) - ); - } + state.assembly.insert_new_sphere(); } ) { "+" } button( @@ -188,13 +178,20 @@ pub fn AddRemove() -> View { }, on:click=|_| { let state = use_context::(); - let subjects = state.selection.with( - |sel| { - let subject_vec: Vec<_> = sel.into_iter().collect(); - (subject_vec[0].clone(), subject_vec[1].clone()) - } + let subjects: [_; 2] = state.selection.with( + // the button is only enabled when two elements are + // selected, so we know the cast to a two-element array + // will succeed + |sel| sel + .clone() + .into_iter() + .collect::>() + .try_into() + .unwrap() + ); + state.assembly.insert_regulator( + InversiveDistanceRegulator::new(subjects, &state.assembly) ); - state.assembly.insert_new_regulator(subjects); state.selection.update(|sel| sel.clear()); } ) { "🔗" } diff --git a/app-proto/src/assembly.rs b/app-proto/src/assembly.rs index 18176df..5c926ca 100644 --- a/app-proto/src/assembly.rs +++ b/app-proto/src/assembly.rs @@ -1,12 +1,21 @@ use nalgebra::{DMatrix, DVector, DVectorView, Vector3}; use rustc_hash::FxHashMap; use slab::Slab; -use std::{collections::BTreeSet, sync::atomic::{AtomicU64, Ordering}}; +use std::{collections::BTreeSet, rc::Rc, sync::atomic::{AtomicU64, Ordering}}; use sycamore::prelude::*; use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */ use crate::{ - engine::{Q, local_unif_to_std, realize_gram, ConfigSubspace, PartialMatrix}, + engine::{ + Q, + change_half_curvature, + local_unif_to_std, + realize_gram, + sphere, + ConfigSubspace, + ConstraintProblem + }, + outline::OutlineItem, specified::SpecifiedValue }; @@ -23,6 +32,10 @@ pub type ElementColor = [f32; 3]; // each assembly has a key that identifies it within the sesssion static NEXT_ELEMENT_SERIAL: AtomicU64 = AtomicU64::new(0); +pub trait ProblemPoser { + fn pose(&self, problem: &mut ConstraintProblem, elts: &Slab); +} + #[derive(Clone, PartialEq)] pub struct Element { pub id: String, @@ -30,8 +43,8 @@ pub struct Element { pub color: ElementColor, pub representation: Signal>, - // All regulators with this element as a subject. The assembly owning - // this element is responsible for keeping this set up to date. + // the regulators this element is subject to. the assembly that owns the + // element is responsible for keeping this set up to date pub regulators: Signal>, // a serial number, assigned by `Element::new`, that uniquely identifies @@ -45,6 +58,8 @@ pub struct Element { } impl Element { + const CURVATURE_COMPONENT: usize = 3; + pub fn new( id: String, label: String, @@ -117,13 +132,148 @@ impl Element { } } -#[derive(Clone, Copy)] -pub struct Regulator { - pub subjects: (ElementKey, ElementKey), +impl ProblemPoser for Element { + fn pose(&self, problem: &mut ConstraintProblem, _elts: &Slab) { + let index = self.column_index.expect( + format!("Element \"{}\" should be indexed before writing problem data", self.id).as_str() + ); + problem.gram.push_sym(index, index, 1.0); + problem.guess.set_column(index, &self.representation.get_clone_untracked()); + } +} + +pub trait Regulator: ProblemPoser + OutlineItem { + fn subjects(&self) -> Vec; + fn measurement(&self) -> ReadSignal; + fn set_point(&self) -> Signal; + + // this method is used to responsively precondition the assembly for + // realization when the regulator becomes a constraint, or is edited while + // acting as a constraint. it should track the set point, do any desired + // preconditioning when the set point is present, and use its return value + // to report whether the set is present. the default implementation does no + // preconditioning + fn try_activate(&self, _assembly: &Assembly) -> bool { + self.set_point().with(|set_pt| set_pt.is_present()) + } +} + +pub struct InversiveDistanceRegulator { + pub subjects: [ElementKey; 2], pub measurement: ReadSignal, pub set_point: Signal } +impl InversiveDistanceRegulator { + pub fn new(subjects: [ElementKey; 2], assembly: &Assembly) -> InversiveDistanceRegulator { + let measurement = assembly.elements.map( + move |elts| { + let representations = subjects.map(|subj| elts[subj].representation); + representations[0].with(|rep_0| + representations[1].with(|rep_1| + rep_0.dot(&(&*Q * rep_1)) + ) + ) + } + ); + + let set_point = create_signal(SpecifiedValue::from_empty_spec()); + + InversiveDistanceRegulator { subjects, measurement, set_point } + } +} + +impl Regulator for InversiveDistanceRegulator { + fn subjects(&self) -> Vec { + self.subjects.into() + } + + fn measurement(&self) -> ReadSignal { + self.measurement + } + + fn set_point(&self) -> Signal { + self.set_point + } +} + +impl ProblemPoser for InversiveDistanceRegulator { + fn pose(&self, problem: &mut ConstraintProblem, elts: &Slab) { + self.set_point.with_untracked(|set_pt| { + if let Some(val) = set_pt.value { + let [row, col] = self.subjects.map( + |subj| elts[subj].column_index.expect( + "Subjects should be indexed before inversive distance regulator writes problem data" + ) + ); + problem.gram.push_sym(row, col, val); + } + }); + } +} + +pub struct HalfCurvatureRegulator { + pub subject: ElementKey, + pub measurement: ReadSignal, + pub set_point: Signal +} + +impl HalfCurvatureRegulator { + pub fn new(subject: ElementKey, assembly: &Assembly) -> HalfCurvatureRegulator { + let measurement = assembly.elements.map( + move |elts| elts[subject].representation.with( + |rep| rep[Element::CURVATURE_COMPONENT] + ) + ); + + let set_point = create_signal(SpecifiedValue::from_empty_spec()); + + HalfCurvatureRegulator { subject, measurement, set_point } + } +} + +impl Regulator for HalfCurvatureRegulator { + fn subjects(&self) -> Vec { + vec![self.subject] + } + + fn measurement(&self) -> ReadSignal { + self.measurement + } + + fn set_point(&self) -> Signal { + self.set_point + } + + fn try_activate(&self, assembly: &Assembly) -> bool { + match self.set_point.with(|set_pt| set_pt.value) { + Some(half_curv) => { + let representation = assembly.elements.with_untracked( + |elts| elts[self.subject].representation + ); + representation.update( + |rep| change_half_curvature(rep, half_curv) + ); + true + } + None => false + } + } +} + +impl ProblemPoser for HalfCurvatureRegulator { + fn pose(&self, problem: &mut ConstraintProblem, elts: &Slab) { + self.set_point.with_untracked(|set_pt| { + if let Some(val) = set_pt.value { + let col = elts[self.subject].column_index.expect( + "Subject should be indexed before half-curvature regulator writes problem data" + ); + problem.frozen.push(Element::CURVATURE_COMPONENT, col, val); + } + }); + } +} + // the velocity is expressed in uniform coordinates pub struct ElementMotion<'a> { pub key: ElementKey, @@ -137,7 +287,7 @@ type AssemblyMotion<'a> = Vec>; pub struct Assembly { // elements and regulators pub elements: Signal>, - pub regulators: Signal>, + pub regulators: Signal>>, // solution variety tangent space. the basis vectors are stored in // configuration matrix format, ordered according to the elements' column @@ -167,26 +317,33 @@ impl Assembly { // --- inserting elements and regulators --- - // insert an element into the assembly without checking whether we already + // insert a sphere into the assembly without checking whether we already // have an element with the same identifier. any element that does have the // same identifier will get kicked out of the `elements_by_id` index - fn insert_element_unchecked(&self, elt: Element) { + fn insert_sphere_unchecked(&self, elt: Element) -> ElementKey { + // insert the sphere 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)); + + // regulate the sphere's curvature + self.insert_regulator(HalfCurvatureRegulator::new(key, &self)); + + key } - pub fn try_insert_element(&self, elt: Element) -> bool { + pub fn try_insert_sphere(&self, elt: Element) -> Option { let can_insert = self.elements_by_id.with_untracked( |elts_by_id| !elts_by_id.contains_key(&elt.id) ); if can_insert { - self.insert_element_unchecked(elt); + Some(self.insert_sphere_unchecked(elt)) + } else { + None } - can_insert } - pub fn insert_new_element(&self) { + pub fn insert_new_sphere(&self) { // find the next unused identifier in the default sequence let mut id_num = 1; let mut id = format!("sphere{}", id_num); @@ -197,70 +354,69 @@ impl Assembly { id = format!("sphere{}", id_num); } - // create and insert a new element - self.insert_element_unchecked( + // create and insert a sphere + let _ = self.insert_sphere_unchecked( Element::new( id, format!("Sphere {}", id_num), [0.75_f32, 0.75_f32, 0.75_f32], - DVector::::from_column_slice(&[0.0, 0.0, 0.0, 0.5, -0.5]) + sphere(0.0, 0.0, 0.0, 1.0) ) ); } - fn insert_regulator(&self, regulator: Regulator) { - let subjects = regulator.subjects; - let key = self.regulators.update(|regs| regs.insert(regulator)); - let subject_regulators = self.elements.with( - |elts| (elts[subjects.0].regulators, elts[subjects.1].regulators) + pub fn insert_regulator(&self, regulator: T) { + // add the regulator to the assembly's regulator list + let regulator_rc = Rc::new(regulator); + let key = self.regulators.update( + |regs| regs.insert(regulator_rc.clone()) ); - subject_regulators.0.update(|regs| regs.insert(key)); - subject_regulators.1.update(|regs| regs.insert(key)); - } - - pub fn insert_new_regulator(self, subjects: (ElementKey, ElementKey)) { - // create and insert a new regulator - let measurement = self.elements.map( - move |elts| { - let reps = ( - elts[subjects.0].representation.get_clone(), - elts[subjects.1].representation.get_clone() - ); - reps.0.dot(&(&*Q * reps.1)) + + // add the regulator to each subject's regulator list + let subjects = regulator_rc.subjects(); + let subject_regulators: Vec<_> = self.elements.with_untracked( + |elts| subjects.into_iter().map( + |subj| elts[subj].regulators + ).collect() + ); + for regulators in subject_regulators { + regulators.update(|regs| regs.insert(key)); + } + + // update the realization when the regulator becomes a constraint, or is + // edited while acting as a constraint + let self_for_effect = self.clone(); + create_effect(move || { + /* DEBUG */ + // log the regulator update + console::log_1(&JsValue::from( + format!("Updated regulator with subjects {:?}", regulator_rc.subjects()) + )); + + if regulator_rc.try_activate(&self_for_effect) { + self_for_effect.realize(); } - ); - let set_point = create_signal(SpecifiedValue::from_empty_spec()); - self.insert_regulator(Regulator { - subjects: subjects, - measurement: measurement, - set_point: set_point }); /* DEBUG */ // print an updated list of regulators console::log_1(&JsValue::from("Regulators:")); - self.regulators.with(|regs| { + self.regulators.with_untracked(|regs| { for (_, reg) in regs.into_iter() { - console::log_5( - &JsValue::from(" "), - &JsValue::from(reg.subjects.0), - &JsValue::from(reg.subjects.1), - &JsValue::from(":"), - ®.set_point.with_untracked( - |set_pt| JsValue::from(set_pt.spec.as_str()) + console::log_1(&JsValue::from(format!( + " {:?}: {}", + reg.subjects(), + reg.set_point().with_untracked( + |set_pt| { + let spec = &set_pt.spec; + if spec.is_empty() { + "__".to_string() + } else { + spec.clone() + } + } ) - ); - } - }); - - // update the realization when the regulator becomes a constraint, or is - // edited while acting as a constraint - create_effect(move || { - console::log_1(&JsValue::from( - format!("Updated constraint with subjects ({}, {})", subjects.0, subjects.1) - )); - if set_point.with(|set_pt| set_pt.is_present()) { - self.realize(); + ))); } }); } @@ -275,55 +431,39 @@ impl Assembly { } }); - // set up the Gram matrix and the initial configuration matrix - let (gram, guess) = self.elements.with_untracked(|elts| { - // set up the off-diagonal part of the Gram matrix - let mut gram_to_be = PartialMatrix::new(); + // set up the constraint problem + let problem = self.elements.with_untracked(|elts| { + let mut problem = ConstraintProblem::new(elts.len()); + for (_, elt) in elts { + elt.pose(&mut problem, elts); + } self.regulators.with_untracked(|regs| { for (_, reg) in regs { - reg.set_point.with_untracked(|set_pt| { - if let Some(val) = set_pt.value { - let subjects = reg.subjects; - let row = elts[subjects.0].column_index.unwrap(); - let col = elts[subjects.1].column_index.unwrap(); - gram_to_be.push_sym(row, col, val); - } - }); + reg.pose(&mut problem, elts); } }); - - // set up the initial configuration matrix and the diagonal of the - // Gram matrix - let mut guess_to_be = DMatrix::::zeros(5, elts.len()); - for (_, elt) in elts { - 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()); - } - - (gram_to_be, guess_to_be) + problem }); /* DEBUG */ // log the Gram matrix console::log_1(&JsValue::from("Gram matrix:")); - gram.log_to_console(); + problem.gram.log_to_console(); /* DEBUG */ // log the initial configuration matrix console::log_1(&JsValue::from("Old configuration:")); - for j in 0..guess.nrows() { + for j in 0..problem.guess.nrows() { let mut row_str = String::new(); - for k in 0..guess.ncols() { - row_str.push_str(format!(" {:>8.3}", guess[(j, k)]).as_str()); + 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( - &gram, guess, &[], - 1.0e-12, 0.5, 0.9, 1.1, 200, 110 + &problem, 1.0e-12, 0.5, 0.9, 1.1, 200, 110 ); /* DEBUG */ @@ -458,4 +598,48 @@ impl Assembly { // sync self.realize(); } +} + +#[cfg(test)] +mod tests { + use crate::engine; + + use super::*; + + #[test] + #[should_panic(expected = "Element \"sphere\" should be indexed before writing problem data")] + fn unindexed_element_test() { + let _ = create_root(|| { + Element::new( + "sphere".to_string(), + "Sphere".to_string(), + [1.0_f32, 1.0_f32, 1.0_f32], + engine::sphere(0.0, 0.0, 0.0, 1.0) + ).pose(&mut ConstraintProblem::new(1), &Slab::new()); + }); + } + + #[test] + #[should_panic(expected = "Subjects should be indexed before inversive distance regulator writes problem data")] + fn unindexed_subject_test_inversive_distance() { + let _ = create_root(|| { + let mut elts = Slab::new(); + let subjects = [0, 1].map(|k| { + elts.insert( + Element::new( + format!("sphere{k}"), + format!("Sphere {k}"), + [1.0_f32, 1.0_f32, 1.0_f32], + engine::sphere(0.0, 0.0, 0.0, 1.0) + ) + ) + }); + elts[subjects[0]].column_index = Some(0); + InversiveDistanceRegulator { + subjects: subjects, + measurement: create_memo(|| 0.0), + set_point: create_signal(SpecifiedValue::try_from("0.0".to_string()).unwrap()) + }.pose(&mut ConstraintProblem::new(2), &elts); + }); + } } \ No newline at end of file diff --git a/app-proto/src/engine.rs b/app-proto/src/engine.rs index 35f898c..869a7de 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -35,9 +35,43 @@ pub fn sphere_with_offset(dir_x: f64, dir_y: f64, dir_z: f64, off: f64, curv: f6 ]) } +// 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 +pub fn change_half_curvature(rep: &mut DVector, half_curv: f64) { + // set the sphere's half-curvature to the desired value + rep[3] = half_curv; + + // restore normalization by contracting toward the curvature axis + const SIZE_THRESHOLD: f64 = 1e-9; + let half_q_lt = -2.0 * half_curv * rep[4]; + let half_q_lt_sq = half_q_lt * half_q_lt; + let mut spatial = rep.fixed_rows_mut::<3>(0); + let q_sp = spatial.norm_squared(); + if q_sp < SIZE_THRESHOLD && half_q_lt_sq < SIZE_THRESHOLD { + spatial.copy_from_slice( + &[0.0, 0.0, (1.0 - 2.0 * half_q_lt).sqrt()] + ); + } else { + let scaling = half_q_lt + (q_sp + half_q_lt_sq).sqrt(); + spatial.scale_mut(1.0 / scaling); + rep[4] /= scaling; + } + + /* DEBUG */ + // verify normalization + let rep_for_debug = rep.clone(); + console::log_1(&JsValue::from( + format!( + "Sphere self-product after curvature change: {}", + rep_for_debug.dot(&(&*Q * &rep_for_debug)) + ) + )); +} + // --- partial matrices --- -struct MatrixEntry { +pub struct MatrixEntry { index: (usize, usize), value: f64 } @@ -49,42 +83,72 @@ impl PartialMatrix { PartialMatrix(Vec::::new()) } - pub fn push_sym(&mut self, row: usize, col: usize, value: f64) { + pub fn push(&mut self, row: usize, col: usize, value: f64) { let PartialMatrix(entries) = self; entries.push(MatrixEntry { index: (row, col), value: value }); + } + + pub fn push_sym(&mut self, row: usize, col: usize, value: f64) { + self.push(row, col, value); if row != col { - entries.push(MatrixEntry { index: (col, row), value: value }); + self.push(col, row, value); } } /* DEBUG */ pub fn log_to_console(&self) { - let PartialMatrix(entries) = self; - for ent in entries { - let ent_str = format!(" {} {} {}", ent.index.0, ent.index.1, ent.value); - console::log_1(&JsValue::from(ent_str.as_str())); + for &MatrixEntry { index: (row, col), value } in self { + console::log_1(&JsValue::from( + format!(" {} {} {}", row, col, value) + )); } } + fn freeze(&self, a: &DMatrix) -> DMatrix { + let mut result = a.clone(); + for &MatrixEntry { index, value } in self { + result[index] = value; + } + result + } + fn proj(&self, a: &DMatrix) -> DMatrix { let mut result = DMatrix::::zeros(a.nrows(), a.ncols()); - let PartialMatrix(entries) = self; - for ent in entries { - result[ent.index] = a[ent.index]; + for &MatrixEntry { index, .. } in self { + result[index] = a[index]; } result } fn sub_proj(&self, rhs: &DMatrix) -> DMatrix { let mut result = DMatrix::::zeros(rhs.nrows(), rhs.ncols()); - let PartialMatrix(entries) = self; - for ent in entries { - result[ent.index] = ent.value - rhs[ent.index]; + for &MatrixEntry { index, value } in self { + result[index] = value - rhs[index]; } result } } +impl IntoIterator for PartialMatrix { + type Item = MatrixEntry; + type IntoIter = std::vec::IntoIter; + + fn into_iter(self) -> Self::IntoIter { + let PartialMatrix(entries) = self; + entries.into_iter() + } +} + +impl<'a> IntoIterator for &'a PartialMatrix { + type Item = &'a MatrixEntry; + type IntoIter = std::slice::Iter<'a, MatrixEntry>; + + fn into_iter(self) -> Self::IntoIter { + let PartialMatrix(entries) = self; + entries.into_iter() + } +} + // --- configuration subspaces --- #[derive(Clone)] @@ -195,6 +259,34 @@ impl DescentHistory { } } +// --- constraint problems --- + +pub struct ConstraintProblem { + pub gram: PartialMatrix, + pub frozen: PartialMatrix, + pub guess: DMatrix, +} + +impl ConstraintProblem { + pub fn new(element_count: usize) -> ConstraintProblem { + const ELEMENT_DIM: usize = 5; + ConstraintProblem { + gram: PartialMatrix::new(), + frozen: PartialMatrix::new(), + guess: DMatrix::::zeros(ELEMENT_DIM, element_count) + } + } + + #[cfg(feature = "dev")] + pub fn from_guess(guess_columns: &[DVector]) -> ConstraintProblem { + ConstraintProblem { + gram: PartialMatrix::new(), + frozen: PartialMatrix::new(), + guess: DMatrix::from_columns(guess_columns) + } + } +} + // --- gram matrix realization --- // the Lorentz form @@ -286,12 +378,12 @@ fn seek_better_config( None } -// seek a matrix `config` for which `config' * Q * config` matches the partial -// matrix `gram`. use gradient descent starting from `guess` +// seek a matrix `config` that matches the partial matrix `problem.frozen` and +// has `config' * Q * config` matching the partial matrix `problem.gram`. start +// at `problem.guess`, set the frozen entries to their desired values, and then +// use a regularized Newton's method to seek the desired Gram matrix pub fn realize_gram( - gram: &PartialMatrix, - guess: DMatrix, - frozen: &[(usize, usize)], + problem: &ConstraintProblem, scaled_tol: f64, min_efficiency: f64, backoff: f64, @@ -299,6 +391,11 @@ pub fn realize_gram( max_descent_steps: i32, max_backoff_steps: i32 ) -> (DMatrix, ConfigSubspace, bool, DescentHistory) { + // destructure the problem data + let ConstraintProblem { + gram, guess, frozen + } = problem; + // start the descent history let mut history = DescentHistory::new(); @@ -313,11 +410,11 @@ pub fn realize_gram( // convert the frozen indices to stacked format let frozen_stacked: Vec = frozen.into_iter().map( - |index| index.1*element_dim + index.0 + |MatrixEntry { index: (row, col), .. }| col*element_dim + row ).collect(); - // use Newton's method with backtracking and gradient descent backup - let mut state = SearchState::from_config(gram, guess); + // use a regularized Newton's method with backtracking + let mut state = SearchState::from_config(gram, frozen.freeze(guess)); let mut hess = DMatrix::zeros(element_dim, assembly_dim); for _ in 0..max_descent_steps { // find the negative gradient of the loss function @@ -415,7 +512,7 @@ pub fn realize_gram( #[cfg(feature = "dev")] pub mod examples { - use std::{array, f64::consts::PI}; + use std::f64::consts::PI; use super::*; @@ -428,35 +525,7 @@ pub mod examples { // https://www.nippon.com/en/japan-topics/c12801/ // pub fn realize_irisawa_hexlet(scaled_tol: f64) -> (DMatrix, ConfigSubspace, bool, DescentHistory) { - let gram = { - let mut gram_to_be = PartialMatrix::new(); - for s in 0..9 { - // each sphere is represented by a spacelike vector - gram_to_be.push_sym(s, s, 1.0); - - // the circumscribing sphere is tangent to all of the other - // spheres, with matching orientation - if s > 0 { - gram_to_be.push_sym(0, s, 1.0); - } - - if s > 2 { - // each chain sphere is tangent to the "sun" and "moon" - // spheres, with opposing orientation - for n in 1..3 { - gram_to_be.push_sym(s, n, -1.0); - } - - // each chain sphere is tangent to the next chain sphere, - // with opposing orientation - let s_next = 3 + (s-2) % 6; - gram_to_be.push_sym(s, s_next, -1.0); - } - } - gram_to_be - }; - - let guess = DMatrix::from_columns( + let mut problem = ConstraintProblem::from_guess( [ sphere(0.0, 0.0, 0.0, 15.0), sphere(0.0, 0.0, -9.0, 5.0), @@ -471,42 +540,45 @@ pub mod examples { ).collect::>().as_slice() ); + for s in 0..9 { + // each sphere is represented by a spacelike vector + problem.gram.push_sym(s, s, 1.0); + + // the circumscribing sphere is tangent to all of the other + // spheres, with matching orientation + if s > 0 { + problem.gram.push_sym(0, s, 1.0); + } + + if s > 2 { + // each chain sphere is tangent to the "sun" and "moon" + // spheres, with opposing orientation + for n in 1..3 { + problem.gram.push_sym(s, n, -1.0); + } + + // each chain sphere is tangent to the next chain sphere, + // with opposing orientation + let s_next = 3 + (s-2) % 6; + problem.gram.push_sym(s, s_next, -1.0); + } + } + // the frozen entries fix the radii of the circumscribing sphere, the // "sun" and "moon" spheres, and one of the chain spheres - let frozen: [(usize, usize); 4] = array::from_fn(|k| (3, k)); + for k in 0..4 { + problem.frozen.push(3, k, problem.guess[(3, k)]); + } - realize_gram( - &gram, guess, &frozen, - scaled_tol, 0.5, 0.9, 1.1, 200, 110 - ) + realize_gram(&problem, scaled_tol, 0.5, 0.9, 1.1, 200, 110) } // set up a kaleidocycle, made of points with fixed distances between them, // and find its tangent space pub fn realize_kaleidocycle(scaled_tol: f64) -> (DMatrix, ConfigSubspace, bool, DescentHistory) { - 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( + const N_HINGES: usize = 6; + let mut problem = ConstraintProblem::from_guess( + (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; @@ -519,16 +591,30 @@ pub mod examples { point(x_vert, y_vert, 0.5) ] } - ).collect::>(); - DMatrix::from_columns(&guess_elts) - }; + ).collect::>().as_slice() + ); - let frozen: [_; N_POINTS] = array::from_fn(|k| (3, k)); + const N_POINTS: usize = 2 * N_HINGES; + 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 { + problem.gram.push_sym(block + j, block + k, if j == k { 0.0 } else { -0.5 }); + } + + // non-hinge edges + for k in 0..2 { + problem.gram.push_sym(block + j, block_next + k, -0.625); + } + } + } - realize_gram( - &gram, guess, &frozen, - scaled_tol, 0.5, 0.9, 1.1, 200, 110 - ) + for k in 0..N_POINTS { + problem.frozen.push(3, k, problem.guess[(3, k)]) + } + + realize_gram(&problem, scaled_tol, 0.5, 0.9, 1.1, 200, 110) } } @@ -539,6 +625,25 @@ mod tests { use super::{*, examples::*}; + #[test] + fn freeze_test() { + let frozen = PartialMatrix(vec![ + MatrixEntry { index: (0, 0), value: 14.0 }, + MatrixEntry { index: (0, 2), value: 28.0 }, + MatrixEntry { index: (1, 1), value: 42.0 }, + MatrixEntry { index: (1, 2), value: 49.0 } + ]); + let config = DMatrix::::from_row_slice(2, 3, &[ + 1.0, 2.0, 3.0, + 4.0, 5.0, 6.0 + ]); + let expected_result = DMatrix::::from_row_slice(2, 3, &[ + 14.0, 2.0, 28.0, + 4.0, 42.0, 49.0 + ]); + assert_eq!(frozen.freeze(&config), expected_result); + } + #[test] fn sub_proj_test() { let target = PartialMatrix(vec![ @@ -560,18 +665,12 @@ mod tests { #[test] fn zero_loss_test() { - let gram = PartialMatrix({ - let mut entries = Vec::::new(); - for j in 0..3 { - for k in 0..3 { - entries.push(MatrixEntry { - index: (j, k), - value: if j == k { 1.0 } else { -1.0 } - }); - } + let mut gram = PartialMatrix::new(); + for j in 0..3 { + for k in 0..3 { + gram.push(j, k, if j == k { 1.0 } else { -1.0 }); } - entries - }); + } let config = { let a = 0.75_f64.sqrt(); DMatrix::from_columns(&[ @@ -584,37 +683,33 @@ mod tests { assert!(state.loss.abs() < f64::EPSILON); } + /* TO DO */ // at the frozen indices, the optimization steps should have exact zeros, - // and the realized configuration should match the initial guess + // and the realized configuration should have the desired values #[test] fn frozen_entry_test() { - let gram = { - let mut gram_to_be = PartialMatrix::new(); - for j in 0..2 { - for k in j..2 { - gram_to_be.push_sym(j, k, if (j, k) == (1, 1) { 1.0 } else { 0.0 }); - } - } - gram_to_be - }; - let guess = DMatrix::from_columns(&[ + let mut problem = ConstraintProblem::from_guess(&[ point(0.0, 0.0, 2.0), - sphere(0.0, 0.0, 0.0, 1.0) + sphere(0.0, 0.0, 0.0, 0.95) ]); - let frozen = [(3, 0), (3, 1)]; - println!(); + for j in 0..2 { + for k in j..2 { + problem.gram.push_sym(j, k, if (j, k) == (1, 1) { 1.0 } else { 0.0 }); + } + } + problem.frozen.push(3, 0, problem.guess[(3, 0)]); + problem.frozen.push(3, 1, 0.5); let (config, _, success, history) = realize_gram( - &gram, guess.clone(), &frozen, - 1.0e-12, 0.5, 0.9, 1.1, 200, 110 + &problem, 1.0e-12, 0.5, 0.9, 1.1, 200, 110 ); assert_eq!(success, true); for base_step in history.base_step.into_iter() { - for index in frozen { + for &MatrixEntry { index, .. } in &problem.frozen { assert_eq!(base_step[index], 0.0); } } - for index in frozen { - assert_eq!(config[index], guess[index]); + for MatrixEntry { index, value } in problem.frozen { + assert_eq!(config[index], value); } } @@ -635,34 +730,32 @@ mod tests { #[test] fn tangent_test_three_spheres() { const SCALED_TOL: f64 = 1.0e-12; - 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(&[ + const ELEMENT_DIM: usize = 5; + let mut problem = ConstraintProblem::from_guess(&[ 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)); + for j in 0..3 { + for k in j..3 { + problem.gram.push_sym(j, k, if j == k { 1.0 } else { -1.0 }); + } + } + for n in 0..ELEMENT_DIM { + problem.frozen.push(n, 0, problem.guess[(n, 0)]); + } let (config, tangent, success, history) = realize_gram( - &gram, guess.clone(), &frozen, - SCALED_TOL, 0.5, 0.9, 1.1, 200, 110 + &problem, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110 ); - assert_eq!(config, guess); + assert_eq!(config, problem.guess); assert_eq!(success, true); assert_eq!(history.scaled_loss.len(), 1); // list some motions that should form a basis for the tangent space of // the solution variety const UNIFORM_DIM: usize = 4; - let element_dim = guess.nrows(); - let assembly_dim = guess.ncols(); + let element_dim = problem.guess.nrows(); + let assembly_dim = problem.guess.ncols(); let tangent_motions_unif = vec![ basis_matrix((0, 1), UNIFORM_DIM, assembly_dim), basis_matrix((1, 1), UNIFORM_DIM, assembly_dim), @@ -805,22 +898,17 @@ mod tests { fn proj_equivar_test() { // find a pair of spheres that meet at 120° const SCALED_TOL: f64 = 1.0e-12; - let gram = { - let mut gram_to_be = PartialMatrix::new(); - gram_to_be.push_sym(0, 0, 1.0); - gram_to_be.push_sym(1, 1, 1.0); - gram_to_be.push_sym(0, 1, 0.5); - gram_to_be - }; - let guess_orig = DMatrix::from_columns(&[ + let mut problem_orig = ConstraintProblem::from_guess(&[ sphere(0.0, 0.0, 0.5, 1.0), sphere(0.0, 0.0, -0.5, 1.0) ]); + problem_orig.gram.push_sym(0, 0, 1.0); + problem_orig.gram.push_sym(1, 1, 1.0); + problem_orig.gram.push_sym(0, 1, 0.5); let (config_orig, tangent_orig, success_orig, history_orig) = realize_gram( - &gram, guess_orig.clone(), &[], - SCALED_TOL, 0.5, 0.9, 1.1, 200, 110 + &problem_orig, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110 ); - assert_eq!(config_orig, guess_orig); + assert_eq!(config_orig, problem_orig.guess); assert_eq!(success_orig, true); assert_eq!(history_orig.scaled_loss.len(), 1); @@ -833,11 +921,15 @@ mod tests { sphere(-a, 0.0, 7.0 - a, 1.0) ]) }; + let problem_tfm = ConstraintProblem { + gram: problem_orig.gram, + guess: guess_tfm, + frozen: problem_orig.frozen + }; let (config_tfm, tangent_tfm, success_tfm, history_tfm) = realize_gram( - &gram, guess_tfm.clone(), &[], - SCALED_TOL, 0.5, 0.9, 1.1, 200, 110 + &problem_tfm, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110 ); - assert_eq!(config_tfm, guess_tfm); + assert_eq!(config_tfm, problem_tfm.guess); assert_eq!(success_tfm, true); assert_eq!(history_tfm.scaled_loss.len(), 1); @@ -869,7 +961,7 @@ mod tests { // the comparison tolerance because the transformation seems to // introduce some numerical error const SCALED_TOL_TFM: f64 = 1.0e-9; - let tol_sq = ((guess_orig.nrows() * guess_orig.ncols()) as f64) * SCALED_TOL_TFM * SCALED_TOL_TFM; + let tol_sq = ((problem_orig.guess.nrows() * problem_orig.guess.ncols()) as f64) * SCALED_TOL_TFM * SCALED_TOL_TFM; assert!((motion_proj_tfm - motion_tfm_proj).norm_squared() < tol_sq); } } \ No newline at end of file diff --git a/app-proto/src/outline.rs b/app-proto/src/outline.rs index 002baea..2446337 100644 --- a/app-proto/src/outline.rs +++ b/app-proto/src/outline.rs @@ -1,4 +1,5 @@ use itertools::Itertools; +use std::rc::Rc; use sycamore::prelude::*; use web_sys::{ KeyboardEvent, @@ -9,24 +10,38 @@ use web_sys::{ use crate::{ AppState, assembly, - assembly::{ElementKey, Regulator, RegulatorKey}, + assembly::{ + ElementKey, + HalfCurvatureRegulator, + InversiveDistanceRegulator, + Regulator, + RegulatorKey + }, specified::SpecifiedValue }; // an editable view of a regulator #[component(inline_props)] -fn RegulatorInput(regulator: Regulator) -> View { +fn RegulatorInput(regulator: Rc) -> View { + // get the regulator's measurement and set point signals + let measurement = regulator.measurement(); + let set_point = regulator.set_point(); + + // the `valid` signal tracks whether the last entered value is a valid set + // point specification let valid = create_signal(true); + + // the `value` signal holds the current set point specification let value = create_signal( - regulator.set_point.with_untracked(|set_pt| set_pt.spec.clone()) + set_point.with_untracked(|set_pt| set_pt.spec.clone()) ); - // this closure resets the input value to the regulator's set point - // specification + // this `reset_value` closure resets the input value to the regulator's set + // point specification let reset_value = move || { batch(|| { valid.set(true); - value.set(regulator.set_point.with(|set_pt| set_pt.spec.clone())); + value.set(set_point.with(|set_pt| set_pt.spec.clone())); }) }; @@ -39,7 +54,7 @@ fn RegulatorInput(regulator: Regulator) -> View { r#type="text", class=move || { if valid.get() { - regulator.set_point.with(|set_pt| { + set_point.with(|set_pt| { if set_pt.is_present() { "regulator-input constraint" } else { @@ -50,13 +65,13 @@ fn RegulatorInput(regulator: Regulator) -> View { "regulator-input invalid" } }, - placeholder=regulator.measurement.with(|result| result.to_string()), + placeholder=measurement.with(|result| result.to_string()), bind:value=value, on:change=move |_| { valid.set( match SpecifiedValue::try_from(value.get_clone_untracked()) { Ok(set_pt) => { - regulator.set_point.set(set_pt); + set_point.set(set_pt); true } Err(_) => false @@ -75,26 +90,53 @@ fn RegulatorInput(regulator: Regulator) -> View { } } +pub trait OutlineItem { + fn outline_item(self: Rc, element_key: ElementKey) -> View; +} + +impl OutlineItem for InversiveDistanceRegulator { + fn outline_item(self: Rc, element_key: ElementKey) -> View { + let state = use_context::(); + let other_subject = if self.subjects[0] == element_key { + self.subjects[1] + } else { + self.subjects[0] + }; + let other_subject_label = state.assembly.elements.with( + |elts| elts[other_subject].label.clone() + ); + view! { + li(class="regulator") { + div(class="regulator-label") { (other_subject_label) } + div(class="regulator-type") { "Inversive distance" } + RegulatorInput(regulator=self) + div(class="status") + } + } + } +} + +impl OutlineItem for HalfCurvatureRegulator { + fn outline_item(self: Rc, _element_key: ElementKey) -> View { + view! { + li(class="regulator") { + div(class="regulator-label") // for spacing + div(class="regulator-type") { "Half-curvature" } + RegulatorInput(regulator=self) + div(class="status") + } + } + } +} + // a list item that shows a regulator in an outline view of an element #[component(inline_props)] fn RegulatorOutlineItem(regulator_key: RegulatorKey, element_key: ElementKey) -> View { let state = use_context::(); - let assembly = &state.assembly; - let regulator = assembly.regulators.with(|regs| regs[regulator_key]); - let other_subject = if regulator.subjects.0 == element_key { - regulator.subjects.1 - } else { - regulator.subjects.0 - }; - let other_subject_label = assembly.elements.with(|elts| elts[other_subject].label.clone()); - view! { - li(class="regulator") { - div(class="regulator-label") { (other_subject_label) } - div(class="regulator-type") { "Inversive distance" } - RegulatorInput(regulator=regulator) - div(class="status") - } - } + let regulator = state.assembly.regulators.with( + |regs| regs[regulator_key].clone() + ); + regulator.outline_item(element_key) } // a list item that shows an element in an outline view of an assembly @@ -117,7 +159,15 @@ fn ElementOutlineItem(key: ElementKey, element: assembly::Element) -> View { }; let regulated = element.regulators.map(|regs| regs.len() > 0); let regulator_list = element.regulators.map( - |regs| regs.clone().into_iter().collect() + move |elt_reg_keys| elt_reg_keys + .clone() + .into_iter() + .sorted_by_key( + |®_key| state.assembly.regulators.with( + |regs| regs[reg_key].subjects().len() + ) + ) + .collect() ); let details_node = create_node_ref(); view! {