Add a half-curvature regulator

In the process, add the `OutlineItem` trait so that each regulator can
implement its own outline item view.

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

@@ -1,26 +1,19 @@
-use nalgebra::DMatrix;
+use dyna3::engine::{Q, point, realize_gram, sphere, ConstraintProblem};
-
-use dyna3::engine::{Q, point, realize_gram, sphere, PartialMatrix};
 
 fn main() {
-    let gram = {
+    let mut problem = ConstraintProblem::from_guess(&[
-        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(&[
         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,
+        &problem, 1.0e-12, 0.5, 0.9, 1.1, 200, 110
-        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, ConstraintProblem};
-
-use dyna3::engine::{Q, realize_gram, sphere, PartialMatrix};
 
 fn main() {
-    let gram = {
+    let mut problem = ConstraintProblem::from_guess({
-        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 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, &[],
+        &problem, 1.0e-12, 0.5, 0.9, 1.1, 200, 110
-        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

@@ -166,18 +166,7 @@ pub fn AddRemove() -> View {
             button(
                 on:click=|_| {
                     let state = use_context::<AppState>();
-                    state.assembly.insert_new_element();
+                    state.assembly.insert_new_sphere();
-                    
-                    /* 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)
-                        );
-                    }
                 }
             ) { "+" }
             button(
@@ -188,13 +177,18 @@ pub fn AddRemove() -> View {
                 },
                 on:click=|_| {
                     let state = use_context::<AppState>();
-                    let subjects = state.selection.with(
+                    let subjects: [_; 2] = state.selection.with(
-                        |sel| {
+                        // the button is only enabled when two elements are
-                            let subject_vec: Vec<_> = sel.into_iter().collect();
+                        // selected, so we know the cast to a two-element array
-                            (subject_vec[0].clone(), subject_vec[1].clone())
+                        // will succeed
-                        }
+                        |sel| sel
+                            .clone()
+                            .into_iter()
+                            .collect::<Vec<_>>()
+                            .try_into()
+                            .unwrap()
                     );
-                    state.assembly.insert_new_regulator(subjects);
+                    state.assembly.insert_new_product_regulator(subjects);
                     state.selection.update(|sel| sel.clear());
                 }
             ) { "🔗" }

app-proto/src/assembly.rs

@@ -1,12 +1,20 @@
 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,
+        local_unif_to_std,
+        realize_gram,
+        sphere,
+        ConfigSubspace,
+        ConstraintProblem
+    },
+    outline::OutlineItem,
     specified::SpecifiedValue
 };
 
@@ -30,8 +38,8 @@ pub struct Element {
     pub color: ElementColor,
     pub representation: Signal<DVector<f64>>,
     
-    // All regulators with this element as a subject. The assembly owning
+    // the regulators this element is subject to. the assembly that owns the
-    // this element is responsible for keeping this set up to date.
+    // 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
@@ -115,15 +123,95 @@ impl Element {
             None
         }
     }
+    
+    fn write_to_problem(&self, problem: &mut ConstraintProblem) {
+        if let Some(index) = self.column_index {
+            problem.gram.push_sym(index, index, 1.0);
+            problem.guess.set_column(index, &self.representation.get_clone_untracked());
+        } else {
+            panic!("Tried to write problem data from an unindexed element: \"{}\"", self.id);
+        }
+    }
 }
 
-#[derive(Clone, Copy)]
+pub trait Regulator: OutlineItem {
-pub struct Regulator {
+    // get information
-    pub subjects: (ElementKey, ElementKey),
+    fn subjects(&self) -> Vec<ElementKey>;
+    fn measurement(&self) -> ReadSignal<f64>;
+    fn set_point(&self) -> Signal<SpecifiedValue>;
+    
+    // write problem data
+    fn write_to_problem(&self, problem: &mut ConstraintProblem, elts: &Slab<Element>);
+}
+
+pub struct ProductRegulator {
+    pub subjects: [ElementKey; 2],
     pub measurement: ReadSignal<f64>,
     pub set_point: Signal<SpecifiedValue>
 }
 
+impl Regulator for ProductRegulator {
+    fn subjects(&self) -> Vec<ElementKey> {
+        self.subjects.into()
+    }
+    
+    fn measurement(&self) -> ReadSignal<f64> {
+        self.measurement
+    }
+    
+    fn set_point(&self) -> Signal<SpecifiedValue> {
+        self.set_point
+    }
+    
+    fn write_to_problem(&self, problem: &mut ConstraintProblem, elts: &Slab<Element>) {
+        self.set_point.with_untracked(|set_pt| {
+            if let Some(val) = set_pt.value {
+                let subject_column_indices = self.subjects.map(
+                    |subj| elts[subj].column_index
+                );
+                if let [Some(row), Some(col)] = subject_column_indices {
+                    problem.gram.push_sym(row, col, val);
+                } else {
+                    panic!("Tried to write problem data from a regulator with an unindexed subject");
+                }
+            }
+        });
+    }
+}
+
+pub struct HalfCurvatureRegulator {
+    pub subject: ElementKey,
+    pub measurement: ReadSignal<f64>,
+    pub set_point: Signal<SpecifiedValue>
+}
+
+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 write_to_problem(&self, problem: &mut ConstraintProblem, elts: &Slab<Element>) {
+        self.set_point.with_untracked(|set_pt| {
+            if let Some(val) = set_pt.value {
+                if let Some(col) = elts[self.subject].column_index {
+                    const CURVATURE_COMPONENT: usize = 3;
+                    problem.frozen.push(CURVATURE_COMPONENT, col, val);
+                } else {
+                    panic!("Tried to write problem data from a regulator with an unindexed subject");
+                }
+            }
+        });
+    }
+}
+
 // the velocity is expressed in uniform coordinates
 pub struct ElementMotion<'a> {
     pub key: ElementKey,
@@ -137,7 +225,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
@@ -170,23 +258,25 @@ impl Assembly {
     // insert an element 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_element_unchecked(&self, elt: Element) -> ElementKey {
         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));
+        key
     }
     
-    pub fn try_insert_element(&self, elt: Element) -> bool {
+    pub fn try_insert_element(&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_element_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,67 +287,81 @@ impl Assembly {
             id = format!("sphere{}", id_num);
         }
         
-        // create and insert a new element
+        // create and insert a sphere
-        self.insert_element_unchecked(
+        let key = self.insert_element_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)
             )
         );
+        
+        // create and insert a curvature regulator
+        self.insert_new_half_curvature_regulator(key);
     }
     
-    fn insert_regulator(&self, regulator: Regulator) {
+    fn insert_regulator(&self, regulator: Rc<dyn Regulator>) {
-        let subjects = regulator.subjects;
+        let subjects = regulator.subjects();
-        let key = self.regulators.update(|regs| regs.insert(regulator));
+        let key = self.regulators.update(
-        let subject_regulators = self.elements.with(
+            |regs| regs.insert(regulator)
-            |elts| (elts[subjects.0].regulators, elts[subjects.1].regulators)
         );
-        subject_regulators.0.update(|regs| regs.insert(key));
+        let subject_regulators: Vec<_> = self.elements.with(
-        subject_regulators.1.update(|regs| regs.insert(key));
+            |elts| subjects.into_iter().map(
-    }
+                |subj| elts[subj].regulators
-    
+            ).collect()
-    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))
-            }
         );
-        let set_point = create_signal(SpecifiedValue::from_empty_spec());
+        for regulators in subject_regulators {
-        self.insert_regulator(Regulator {
+            regulators.update(|regs| regs.insert(key));
-            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| {
             for (_, reg) in regs.into_iter() {
-                console::log_5(
+                console::log_1(&JsValue::from(format!(
-                    &JsValue::from(" "),
+                    "  {:?}: {}",
-                    &JsValue::from(reg.subjects.0),
+                    reg.subjects(),
-                    &JsValue::from(reg.subjects.1),
+                    reg.set_point().with_untracked(
-                    &JsValue::from(":"),
+                        |set_pt| {
-                    &reg.set_point.with_untracked(
+                            let spec = &set_pt.spec;
-                        |set_pt| JsValue::from(set_pt.spec.as_str())
+                            if spec.is_empty() {
+                                "__".to_string()
+                            } else {
+                                spec.clone()
+                            }
+                        }
                     )
-                );
+                )));
             }
         });
+    }
+    
+    pub fn insert_new_product_regulator(self, subjects: [ElementKey; 2]) {
+        // create and insert a new product regulator
+        let measurement = self.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());
+        self.insert_regulator(Rc::new(ProductRegulator {
+            subjects: subjects,
+            measurement: measurement,
+            set_point: set_point
+        }));
         
         // 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)
+                format!("Updated regulator with subjects {:?}", subjects)
             ));
             if set_point.with(|set_pt| set_pt.is_present()) {
                 self.realize();
@@ -265,6 +369,64 @@ impl Assembly {
         });
     }
     
+    pub fn insert_new_half_curvature_regulator(self, subject: ElementKey) {
+        // create and insert a new half-curvature regulator
+        let measurement = self.elements.map(
+            move |elts| elts[subject].representation.with(|rep| rep[3])
+        );
+        let set_point = create_signal(SpecifiedValue::from_empty_spec());
+        self.insert_regulator(Rc::new(HalfCurvatureRegulator {
+            subject: subject,
+            measurement: measurement,
+            set_point: set_point
+        }));
+        
+        // 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 regulator with subjects [{}]", subject)
+            ));
+            if let Some(half_curv) = set_point.with(|set_pt| set_pt.value) {
+                let representation = self.elements.with(
+                    |elts| elts[subject].representation
+                );
+                representation.update(|rep| {
+                    // 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))
+                        )
+                    ));
+                });
+                self.realize();
+            }
+        });
+    }
+    
     // --- realization ---
     
     pub fn realize(&self) {
@@ -275,55 +437,39 @@ impl Assembly {
             }
         });
         
-        // set up the Gram matrix and the initial configuration matrix
+        // set up the constraint problem
-        let (gram, guess) = self.elements.with_untracked(|elts| {
+        let problem = self.elements.with_untracked(|elts| {
-            // set up the off-diagonal part of the Gram matrix
+            let mut problem_to_be = ConstraintProblem::new(elts.len());
-            let mut gram_to_be = PartialMatrix::new();
+            for (_, elt) in elts {
+                elt.write_to_problem(&mut problem_to_be);
+            }
             self.regulators.with_untracked(|regs| {
                 for (_, reg) in regs {
-                    reg.set_point.with_untracked(|set_pt| {
+                    reg.write_to_problem(&mut problem_to_be, elts);
-                        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);
-                        }
-                    });
                 }
             });
-            
+            problem_to_be
-            // 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)
         });
         
         /* 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() {
+            for k in 0..problem.guess.ncols() {
-                row_str.push_str(format!(" {:>8.3}", guess[(j, k)]).as_str());
+                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, &[],
+            &problem, 1.0e-12, 0.5, 0.9, 1.1, 200, 110
-            1.0e-12, 0.5, 0.9, 1.1, 200, 110
         );
         
         /* DEBUG */
@@ -458,4 +604,48 @@ impl Assembly {
         // sync
         self.realize();
     }
+}
+
+#[cfg(test)]
+mod tests {
+    use crate::engine;
+    
+    use super::*;
+    
+    #[test]
+    #[should_panic(expected = "Tried to write problem data from an unindexed element: \"sphere\"")]
+    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)
+            ).write_to_problem(&mut ConstraintProblem::new(1));
+        });
+    }
+    
+    #[test]
+    #[should_panic(expected = "Tried to write problem data from a regulator with an unindexed subject")]
+    fn unindexed_subject_test() {
+        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);
+            ProductRegulator {
+                subjects: subjects,
+                measurement: create_memo(|| 0.0),
+                set_point: create_signal(SpecifiedValue::try_from("0.0".to_string()).unwrap())
+            }.write_to_problem(&mut ConstraintProblem::new(2), &elts);
+        });
+    }
 }

app-proto/src/engine.rs

@@ -37,7 +37,7 @@ pub fn sphere_with_offset(dir_x: f64, dir_y: f64, dir_z: f64, off: f64, curv: f6
 
 // --- partial matrices ---
 
-struct MatrixEntry {
+pub struct MatrixEntry {
     index: (usize, usize),
     value: f64
 }
@@ -49,42 +49,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 &MatrixEntry { index: (row, col), value } in self {
-        for ent in entries {
+            console::log_1(&JsValue::from(
-            let ent_str = format!("  {} {} {}", ent.index.0, ent.index.1, ent.value);
+                format!("  {} {} {}", row, col, value)
-            console::log_1(&JsValue::from(ent_str.as_str()));
+            ));
         }
     }
     
+    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 &MatrixEntry { index, .. } in self {
-        for ent in entries {
+            result[index] = a[index];
-            result[ent.index] = a[ent.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 &MatrixEntry { index, value } in self {
-        for ent in entries {
+            result[index] = value - rhs[index];
-            result[ent.index] = ent.value - rhs[ent.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 +225,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 +344,12 @@ fn seek_better_config(
     None
 }
 
-// seek a matrix `config` for which `config' * Q * config` matches the partial
+// seek a matrix `config` that matches the partial matrix `problem.frozen` and
-// matrix `gram`. use gradient descent starting from `guess`
+// 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,
+    problem: &ConstraintProblem,
-    guess: DMatrix<f64>,
-    frozen: &[(usize, usize)],
     scaled_tol: f64,
     min_efficiency: f64,
     backoff: f64,
@@ -299,6 +357,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 +376,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
+    // use a regularized Newton's method with backtracking
-    let mut state = SearchState::from_config(gram, guess);
+    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 +478,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 +491,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 problem = ConstraintProblem::from_guess(
-            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(
             [
                 sphere(0.0, 0.0, 0.0, 15.0),
                 sphere(0.0, 0.0, -9.0, 5.0),
@@ -471,42 +506,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(
+        realize_gram(&problem, scaled_tol, 0.5, 0.9, 1.1, 200, 110)
-            &gram, guess, &frozen,
-            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;
+        const N_HINGES: usize = 6;
-        let gram = {
+        let mut problem = ConstraintProblem::from_guess(
-            let mut gram_to_be = PartialMatrix::new();
+            (0..N_HINGES).step_by(2).flat_map(
-            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(
                 |n| {
                     let ang_hor = (n as f64) * PI/3.0;
                     let ang_vert = ((n + 1) as f64) * PI/3.0;
@@ -519,16 +557,30 @@ pub mod examples {
                         point(x_vert, y_vert, 0.5)
                     ]
                 }
-            ).collect::<Vec<_>>();
+            ).collect::<Vec<_>>().as_slice()
-            DMatrix::from_columns(&guess_elts)
+        );
-        };
         
-        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(
+        for k in 0..N_POINTS {
-            &gram, guess, &frozen,
+            problem.frozen.push(3, k, problem.guess[(3, k)])
-            scaled_tol, 0.5, 0.9, 1.1, 200, 110
+        }
-        )
+        
+        realize_gram(&problem, scaled_tol, 0.5, 0.9, 1.1, 200, 110)
     }
 }
 
@@ -539,6 +591,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 +631,12 @@ mod tests {
     
     #[test]
     fn zero_loss_test() {
-        let gram = PartialMatrix({
+        let mut gram = PartialMatrix::new();
-            let mut entries = Vec::<MatrixEntry>::new();
+        for j in 0..3 {
-            for j in 0..3 {
+            for k in 0..3 {
-                for k in 0..3 {
+                gram.push(j, k, if j == k { 1.0 } else { -1.0 });
-                    entries.push(MatrixEntry {
-                        index: (j, k),
-                        value: if j == k { 1.0 } else { -1.0 }
-                    });
-                }
             }
-            entries
+        }
-        });
         let config = {
             let a = 0.75_f64.sqrt();
             DMatrix::from_columns(&[
@@ -584,37 +649,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 problem = ConstraintProblem::from_guess(&[
-            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(&[
             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)];
+        for j in 0..2 {
-        println!();
+            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,
+            &problem, 1.0e-12, 0.5, 0.9, 1.1, 200, 110
-            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 {
+        for MatrixEntry { index, value } in problem.frozen {
-            assert_eq!(config[index], guess[index]);
+            assert_eq!(config[index], value);
         }
     }
     
@@ -635,34 +696,32 @@ mod tests {
     #[test]
     fn tangent_test_three_spheres() {
         const SCALED_TOL: f64 = 1.0e-12;
-        let gram = {
+        const ELEMENT_DIM: usize = 5;
-            let mut gram_to_be = PartialMatrix::new();
+        let mut problem = ConstraintProblem::from_guess(&[
-            for j in 0..3 {
-                for k in j..3 {
-                    gram_to_be.push_sym(j, k, if j == k { 1.0 } else { -1.0 });
-                }
-            }
-            gram_to_be
-        };
-        let guess = DMatrix::from_columns(&[
             sphere(0.0, 0.0, 0.0, -2.0),
             sphere(0.0, 0.0, 1.0, 1.0),
             sphere(0.0, 0.0, -1.0, 1.0)
         ]);
-        let frozen: [_; 5] = std::array::from_fn(|k| (k, 0));
+        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,
+            &problem, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
-            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 element_dim = problem.guess.nrows();
-        let assembly_dim = guess.ncols();
+        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 +864,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 problem_orig = ConstraintProblem::from_guess(&[
-            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(&[
             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(), &[],
+            &problem_orig, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
-            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 +887,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(), &[],
+            &problem_tfm, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
-            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 +927,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,
+        ProductRegulator,
+        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
+    // this `reset_value` closure resets the input value to the regulator's set
-    // specification
+    // 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 ProductRegulator {
+    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 = state.assembly.regulators.with(
-    let regulator = assembly.regulators.with(|regs| regs[regulator_key]);
+        |regs| regs[regulator_key].clone()
-    let other_subject = if regulator.subjects.0 == element_key {
+    );
-        regulator.subjects.1
+    regulator.outline_item(element_key)
-    } 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")
-        }
-    }
 }
 
 // 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! {