Change conditional panic to expect

Use `expect` to communicate that every element should have a column
index when `pose` is called. Rewrite the panic messages in the style
recommended by the `expect` documentation.

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);

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 {

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::<AppState>();
-                    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::<AppState>();
-                    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::<Vec<_>>()
+                            .try_into()
+                            .unwrap()
+                    );
+                    state.assembly.insert_regulator(
+                        InversiveDistanceRegulator::new(subjects, &state.assembly)
                     );
-                    state.assembly.insert_new_regulator(subjects);
                     state.selection.update(|sel| sel.clear());
                 }
             ) { "🔗" }

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<Element>);
+}
+
 #[derive(Clone, PartialEq)]
 pub struct Element {
     pub id: String,
@@ -30,8 +43,8 @@ pub struct Element {
     pub color: ElementColor,
     pub representation: Signal<DVector<f64>>,
     
-    // 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<BTreeSet<RegulatorKey>>,
     
     // 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<Element>) {
+        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<ElementKey>;
+    fn measurement(&self) -> ReadSignal<f64>;
+    fn set_point(&self) -> Signal<SpecifiedValue>;
+    
+    // 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<f64>,
     pub set_point: Signal<SpecifiedValue>
 }
 
+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<ElementKey> {
+        self.subjects.into()
+    }
+    
+    fn measurement(&self) -> ReadSignal<f64> {
+        self.measurement
+    }
+    
+    fn set_point(&self) -> Signal<SpecifiedValue> {
+        self.set_point
+    }
+}
+
+impl ProblemPoser for InversiveDistanceRegulator {
+    fn pose(&self, problem: &mut ConstraintProblem, elts: &Slab<Element>) {
+        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<f64>,
+    pub set_point: Signal<SpecifiedValue>
+}
+
+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<ElementKey> {
+        vec![self.subject]
+    }
+    
+    fn measurement(&self) -> ReadSignal<f64> {
+        self.measurement
+    }
+    
+    fn set_point(&self) -> Signal<SpecifiedValue> {
+        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<Element>) {
+        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<ElementMotion<'a>>;
 pub struct Assembly {
     // elements and regulators
     pub elements: Signal<Slab<Element>>,
-    pub regulators: Signal<Slab<Regulator>>,
+    pub regulators: Signal<Slab<Rc<dyn Regulator>>>,
     
     // 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<ElementKey> {
         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::<f64>::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<T: Regulator + 'static>(&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(":"),
-                    &reg.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::<f64>::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);
+        });
+    }
 }

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<f64>, 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::<MatrixEntry>::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<f64>) -> DMatrix<f64> {
+        let mut result = a.clone();
+        for &MatrixEntry { index, value } in self {
+            result[index] = value;
+        }
+        result
+    }
+    
     fn proj(&self, a: &DMatrix<f64>) -> DMatrix<f64> {
         let mut result = DMatrix::<f64>::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<f64>) -> DMatrix<f64> {
         let mut result = DMatrix::<f64>::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<Self::Item>;
+    
+    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<f64>,
+}
+
+impl ConstraintProblem {
+    pub fn new(element_count: usize) -> ConstraintProblem {
+        const ELEMENT_DIM: usize = 5;
+        ConstraintProblem {
+            gram: PartialMatrix::new(),
+            frozen: PartialMatrix::new(),
+            guess: DMatrix::<f64>::zeros(ELEMENT_DIM, element_count)
+        }
+    }
+    
+    #[cfg(feature = "dev")]
+    pub fn from_guess(guess_columns: &[DVector<f64>]) -> 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<f64>,
-    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<f64>, 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<usize> = 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<f64>, 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::<Vec<_>>().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<f64>, 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::<Vec<_>>();
-            DMatrix::from_columns(&guess_elts)
-        };
+            ).collect::<Vec<_>>().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::<f64>::from_row_slice(2, 3, &[
+            1.0, 2.0, 3.0,
+            4.0, 5.0, 6.0
+        ]);
+        let expected_result = DMatrix::<f64>::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::<MatrixEntry>::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);
     }
 }

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<dyn Regulator>) -> 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<Self>, element_key: ElementKey) -> View;
+}
+
+impl OutlineItem for InversiveDistanceRegulator {
+    fn outline_item(self: Rc<Self>, element_key: ElementKey) -> View {
+        let state = use_context::<AppState>();
+        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<Self>, _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::<AppState>();
-    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(
+                |&reg_key| state.assembly.regulators.with(
+                    |regs| regs[reg_key].subjects().len()
+                )
+            )
+            .collect()
     );
     let details_node = create_node_ref();
     view! {