Separate test and example for Irisawa hexlet

Put shared code in the conditionally compiled `engine::irisawa` module.

app-proto/Cargo.toml

@@ -1,15 +1,17 @@
 [package]
-name = "sketch-outline"
+name = "dyna3"
 version = "0.1.0"
-authors = ["Aaron"]
+authors = ["Aaron Fenyes", "Glen Whitney"]
 edition = "2021"
 
 [features]
 default = ["console_error_panic_hook"]
+irisawa = []
 
 [dependencies]
 itertools = "0.13.0"
 js-sys = "0.3.70"
+lazy_static = "1.5.0"
 nalgebra = "0.33.0"
 rustc-hash = "2.0.0"
 slab = "0.4.9"
@@ -25,6 +27,7 @@ console_error_panic_hook = { version = "0.1.7", optional = true }
 version = "0.3.69"
 features = [
   'HtmlCanvasElement',
+  'HtmlInputElement',
   'Performance',
   'WebGl2RenderingContext',
   'WebGlBuffer',
@@ -34,7 +37,12 @@ features = [
   'WebGlVertexArrayObject'
 ]
 
+# the self-dependency specifies features to use for tests and examples
+#
+#   https://github.com/rust-lang/cargo/issues/2911#issuecomment-1483256987
+#
 [dev-dependencies]
+dyna3 = { path = ".", default-features = false, features = ["irisawa"] }
 wasm-bindgen-test = "0.3.34"
 
 [profile.release]

app-proto/examples/irisawa-hexlet.rs

@@ -0,0 +1,25 @@
+use dyna3::engine::{Q, irisawa::realize_irisawa_hexlet};
+
+fn main() {
+    const SCALED_TOL: f64 = 1.0e-12;
+    let (config, success, history) = realize_irisawa_hexlet(SCALED_TOL);
+    print!("\nCompleted Gram matrix:{}", config.tr_mul(&*Q) * &config);
+    if success {
+        println!("Target accuracy achieved!");
+    } else {
+        println!("Failed to reach target accuracy");
+    }
+    println!("Steps: {}", history.scaled_loss.len() - 1);
+    println!("Loss: {}", history.scaled_loss.last().unwrap());
+    if success {
+        println!("\nChain diameters:");
+        println!("  {} sun (given)", 1.0 / config[(3, 3)]);
+        for k in 4..9 {
+            println!("  {} sun", 1.0 / config[(3, k)]);
+        }
+    }
+    println!("\nStep │ Loss\n─────┼────────────────────────────────");
+    for (step, scaled_loss) in history.scaled_loss.into_iter().enumerate() {
+        println!("{:<4} │ {}", step, scaled_loss);
+    }
+}

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

@@ -0,0 +1,38 @@
+use nalgebra::DMatrix;
+
+use dyna3::engine::{Q, point, realize_gram, sphere, PartialMatrix};
+
+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(&[
+        point(0.0, 0.0, 2.0),
+        sphere(0.0, 0.0, 0.0, 1.0)
+    ]);
+    let frozen = [(3, 0)];
+    println!();
+    let (config, success, history) = realize_gram(
+        &gram, guess, &frozen,
+        1.0e-12, 0.5, 0.9, 1.1, 200, 110
+    );
+    print!("\nCompleted Gram matrix:{}", config.tr_mul(&*Q) * &config);
+    print!("Configuration:{}", config);
+    if success {
+        println!("Target accuracy achieved!");
+    } else {
+        println!("Failed to reach target accuracy");
+    }
+    println!("Steps: {}", history.scaled_loss.len() - 1);
+    println!("Loss: {}", history.scaled_loss.last().unwrap());
+    println!("\nStep │ Loss\n─────┼────────────────────────────────");
+    for (step, scaled_loss) in history.scaled_loss.into_iter().enumerate() {
+        println!("{:<4} │ {}", step, scaled_loss);
+    }
+}

app-proto/examples/three-spheres.rs

@@ -0,0 +1,40 @@
+use nalgebra::DMatrix;
+
+use dyna3::engine::{Q, realize_gram, sphere, PartialMatrix};
+
+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 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)
+        ])
+    };
+    println!();
+    let (config, success, history) = realize_gram(
+        &gram, guess, &[],
+        1.0e-12, 0.5, 0.9, 1.1, 200, 110
+    );
+    print!("\nCompleted Gram matrix:{}", config.tr_mul(&*Q) * &config);
+    if success {
+        println!("Target accuracy achieved!");
+    } else {
+        println!("Failed to reach target accuracy");
+    }
+    println!("Steps: {}", history.scaled_loss.len() - 1);
+    println!("Loss: {}", history.scaled_loss.last().unwrap());
+    println!("\nStep │ Loss\n─────┼────────────────────────────────");
+    for (step, scaled_loss) in history.scaled_loss.into_iter().enumerate() {
+        println!("{:<4} │ {}", step, scaled_loss);
+    }
+}

app-proto/index.html

@@ -2,8 +2,10 @@
 <html>
   <head>
     <meta charset="utf-8"/>
-    <title>Sketch outline</title>
+    <title>dyna3</title>
     <link data-trunk rel="css" href="main.css"/>
+    <link href="https://fonts.bunny.net/css?family=fira-sans:ital,wght@0,400;1,400&display=swap" rel="stylesheet">
+    <link href="https://fonts.bunny.net/css?family=noto-emoji:wght@400&text=%f0%9f%94%97%e2%9a%a0&display=swap" rel="stylesheet">
   </head>
   <body></body>
 </html>

app-proto/main.css

@@ -2,6 +2,7 @@ body {
   margin: 0px;
   color: #fcfcfc;
   background-color: #222;
+  font-family: 'Fira Sans', sans-serif;
 }
 
 /* sidebar */
@@ -33,6 +34,11 @@ body {
   font-size: large;
 }
 
+/* KLUDGE */
+#add-remove > button.emoji {
+  font-family: 'Noto Emoji', sans-serif;
+}
+
 /* outline */
 
 #outline {
@@ -93,10 +99,11 @@ details[open]:has(li) .elt-switch::after {
   display: flex;
 }
 
-.elt-rep > div, .cst-rep {
+.elt-rep > div {
   padding: 2px 0px 0px 0px;
   font-size: 10pt;
-  text-align: center;
+  font-variant-numeric: tabular-nums;
+  text-align: right;
   width: 56px;
 }
 
@@ -104,10 +111,38 @@ details[open]:has(li) .elt-switch::after {
   font-style: italic;
 }
 
-.cst > input {
+.cst.invalid {
+  color: #f58fc2;
+}
+
+.cst > input[type=checkbox] {
   margin: 0px 8px 0px 0px;
 }
 
+.cst > input[type=text] {
+  color: inherit;
+  background-color: inherit;
+  border: 1px solid #555;
+  border-radius: 2px;
+}
+
+.cst.invalid > input[type=text] {
+  border-color: #70495c;
+}
+
+.status {
+  width: 20px;
+  padding-left: 4px;
+  text-align: center;
+  font-family: 'Noto Emoji';
+  font-style: normal;
+}
+
+.invalid > .status::after, details:has(.invalid):not([open]) .status::after {
+  content: '⚠';
+  color: #f58fc2;
+}
+
 /* display */
 
 canvas {

app-proto/run-examples

@@ -0,0 +1,9 @@
+# run all Cargo examples, as described here:
+#
+#   Karol Kuczmarski. "Add examples to your Rust libraries"
+#   http://xion.io/post/code/rust-examples.html
+#
+
+cargo run --example irisawa-hexlet
+cargo run --example three-spheres
+cargo run --example point-on-sphere

app-proto/src/add_remove.rs

@@ -1,4 +1,3 @@
-use std::collections::BTreeSet; /* DEBUG */
 use sycamore::prelude::*;
 use web_sys::{console, wasm_bindgen::JsValue};
 
@@ -7,68 +6,52 @@ use crate::{engine, AppState, assembly::{Assembly, Constraint, Element}};
 /* DEBUG */
 fn load_gen_assemb(assembly: &Assembly) {
     let _ = assembly.try_insert_element(
-        Element {
+        Element::new(
-            id: String::from("gemini_a"),
+            String::from("gemini_a"),
-            label: String::from("Castor"),
+            String::from("Castor"),
-            color: [1.00_f32, 0.25_f32, 0.00_f32],
+            [1.00_f32, 0.25_f32, 0.00_f32],
-            rep: engine::sphere(0.5, 0.5, 0.0, 1.0),
+            engine::sphere(0.5, 0.5, 0.0, 1.0)
-            constraints: BTreeSet::default()
+        )
-        }
     );
     let _ = assembly.try_insert_element(
-        Element {
+        Element::new(
-            id: String::from("gemini_b"),
+            String::from("gemini_b"),
-            label: String::from("Pollux"),
+            String::from("Pollux"),
-            color: [0.00_f32, 0.25_f32, 1.00_f32],
+            [0.00_f32, 0.25_f32, 1.00_f32],
-            rep: engine::sphere(-0.5, -0.5, 0.0, 1.0),
+            engine::sphere(-0.5, -0.5, 0.0, 1.0)
-            constraints: BTreeSet::default()
+        )
-        }
     );
     let _ = assembly.try_insert_element(
-        Element {
+        Element::new(
-            id: String::from("ursa_major"),
+            String::from("ursa_major"),
-            label: String::from("Ursa major"),
+            String::from("Ursa major"),
-            color: [0.25_f32, 0.00_f32, 1.00_f32],
+            [0.25_f32, 0.00_f32, 1.00_f32],
-            rep: engine::sphere(-0.5, 0.5, 0.0, 0.75),
+            engine::sphere(-0.5, 0.5, 0.0, 0.75)
-            constraints: BTreeSet::default()
+        )
-        }
     );
     let _ = assembly.try_insert_element(
-        Element {
+        Element::new(
-            id: String::from("ursa_minor"),
+            String::from("ursa_minor"),
-            label: String::from("Ursa minor"),
+            String::from("Ursa minor"),
-            color: [0.25_f32, 1.00_f32, 0.00_f32],
+            [0.25_f32, 1.00_f32, 0.00_f32],
-            rep: engine::sphere(0.5, -0.5, 0.0, 0.5),
+            engine::sphere(0.5, -0.5, 0.0, 0.5)
-            constraints: BTreeSet::default()
+        )
-        }
     );
     let _ = assembly.try_insert_element(
-        Element {
+        Element::new(
-            id: String::from("moon_deimos"),
+            String::from("moon_deimos"),
-            label: String::from("Deimos"),
+            String::from("Deimos"),
-            color: [0.75_f32, 0.75_f32, 0.00_f32],
+            [0.75_f32, 0.75_f32, 0.00_f32],
-            rep: engine::sphere(0.0, 0.15, 1.0, 0.25),
+            engine::sphere(0.0, 0.15, 1.0, 0.25)
-            constraints: BTreeSet::default()
+        )
-        }
     );
     let _ = assembly.try_insert_element(
-        Element {
+        Element::new(
-            id: String::from("moon_phobos"),
+            String::from("moon_phobos"),
-            label: String::from("Phobos"),
+            String::from("Phobos"),
-            color: [0.00_f32, 0.75_f32, 0.50_f32],
+            [0.00_f32, 0.75_f32, 0.50_f32],
-            rep: engine::sphere(0.0, -0.15, -1.0, 0.25),
+            engine::sphere(0.0, -0.15, -1.0, 0.25)
-            constraints: BTreeSet::default()
+        )
-        }
-    );
-    assembly.insert_constraint(
-        Constraint {
-            args: (
-                assembly.elements_by_id.with_untracked(|elts_by_id| elts_by_id["gemini_a"]),
-                assembly.elements_by_id.with_untracked(|elts_by_id| elts_by_id["gemini_b"])
-            ),
-            rep: 0.5,
-            active: create_signal(true)
-        }
     );
 }
 
@@ -76,76 +59,68 @@ fn load_gen_assemb(assembly: &Assembly) {
 fn load_low_curv_assemb(assembly: &Assembly) {
     let a = 0.75_f64.sqrt();
     let _ = assembly.try_insert_element(
-        Element {
+        Element::new(
-            id: "central".to_string(),
+            "central".to_string(),
-            label: "Central".to_string(),
+            "Central".to_string(),
-            color: [0.75_f32, 0.75_f32, 0.75_f32],
+            [0.75_f32, 0.75_f32, 0.75_f32],
-            rep: engine::sphere(0.0, 0.0, 0.0, 1.0),
+            engine::sphere(0.0, 0.0, 0.0, 1.0)
-            constraints: BTreeSet::default()
+        )
-        }
     );
     let _ = assembly.try_insert_element(
-        Element {
+        Element::new(
-            id: "assemb_plane".to_string(),
+            "assemb_plane".to_string(),
-            label: "Assembly plane".to_string(),
+            "Assembly plane".to_string(),
-            color: [0.75_f32, 0.75_f32, 0.75_f32],
+            [0.75_f32, 0.75_f32, 0.75_f32],
-            rep: engine::sphere_with_offset(0.0, 0.0, 1.0, 0.0, 0.0),
+            engine::sphere_with_offset(0.0, 0.0, 1.0, 0.0, 0.0)
-            constraints: BTreeSet::default()
+        )
-        }
     );
     let _ = assembly.try_insert_element(
-        Element {
+        Element::new(
-            id: "side1".to_string(),
+            "side1".to_string(),
-            label: "Side 1".to_string(),
+            "Side 1".to_string(),
-            color: [1.00_f32, 0.00_f32, 0.25_f32],
+            [1.00_f32, 0.00_f32, 0.25_f32],
-            rep: engine::sphere_with_offset(1.0, 0.0, 0.0, 1.0, 0.0),
+            engine::sphere_with_offset(1.0, 0.0, 0.0, 1.0, 0.0)
-            constraints: BTreeSet::default()
+        )
-        }
     );
     let _ = assembly.try_insert_element(
-        Element {
+        Element::new(
-            id: "side2".to_string(),
+            "side2".to_string(),
-            label: "Side 2".to_string(),
+            "Side 2".to_string(),
-            color: [0.25_f32, 1.00_f32, 0.00_f32],
+            [0.25_f32, 1.00_f32, 0.00_f32],
-            rep: engine::sphere_with_offset(-0.5, a, 0.0, 1.0, 0.0),
+            engine::sphere_with_offset(-0.5, a, 0.0, 1.0, 0.0)
-            constraints: BTreeSet::default()
+        )
-        }
     );
     let _ = assembly.try_insert_element(
-        Element {
+        Element::new(
-            id: "side3".to_string(),
+            "side3".to_string(),
-            label: "Side 3".to_string(),
+            "Side 3".to_string(),
-            color: [0.00_f32, 0.25_f32, 1.00_f32],
+            [0.00_f32, 0.25_f32, 1.00_f32],
-            rep: engine::sphere_with_offset(-0.5, -a, 0.0, 1.0, 0.0),
+            engine::sphere_with_offset(-0.5, -a, 0.0, 1.0, 0.0)
-            constraints: BTreeSet::default()
+        )
-        }
     );
     let _ = assembly.try_insert_element(
-        Element {
+        Element::new(
-            id: "corner1".to_string(),
+            "corner1".to_string(),
-            label: "Corner 1".to_string(),
+            "Corner 1".to_string(),
-            color: [0.75_f32, 0.75_f32, 0.75_f32],
+            [0.75_f32, 0.75_f32, 0.75_f32],
-            rep: engine::sphere(-4.0/3.0, 0.0, 0.0, 1.0/3.0),
+            engine::sphere(-4.0/3.0, 0.0, 0.0, 1.0/3.0)
-            constraints: BTreeSet::default()
+        )
-        }
     );
     let _ = assembly.try_insert_element(
-        Element {
+        Element::new(
-            id: "corner2".to_string(),
+            "corner2".to_string(),
-            label: "Corner 2".to_string(),
+            "Corner 2".to_string(),
-            color: [0.75_f32, 0.75_f32, 0.75_f32],
+            [0.75_f32, 0.75_f32, 0.75_f32],
-            rep: engine::sphere(2.0/3.0, -4.0/3.0 * a, 0.0, 1.0/3.0),
+            engine::sphere(2.0/3.0, -4.0/3.0 * a, 0.0, 1.0/3.0)
-            constraints: BTreeSet::default()
+        )
-        }
     );
     let _ = assembly.try_insert_element(
-        Element {
+        Element::new(
-            id: String::from("corner3"),
+            String::from("corner3"),
-            label: String::from("Corner 3"),
+            String::from("Corner 3"),
-            color: [0.75_f32, 0.75_f32, 0.75_f32],
+            [0.75_f32, 0.75_f32, 0.75_f32],
-            rep: engine::sphere(2.0/3.0, 4.0/3.0 * a, 0.0, 1.0/3.0),
+            engine::sphere(2.0/3.0, 4.0/3.0 * a, 0.0, 1.0/3.0)
-            constraints: BTreeSet::default()
+        )
-        }
     );
 }
 
@@ -198,6 +173,7 @@ pub fn AddRemove() -> View {
                 }
             ) { "+" }
             button(
+                class="emoji", /* KLUDGE */
                 disabled={
                     let state = use_context::<AppState>();
                     state.selection.with(|sel| sel.len() != 2)
@@ -210,16 +186,21 @@ pub fn AddRemove() -> View {
                             (arg_vec[0].clone(), arg_vec[1].clone())
                         }
                     );
+                    let rep = create_signal(0.0);
+                    let rep_valid = create_signal(false);
+                    let active = create_signal(true);
                     state.assembly.insert_constraint(Constraint {
                         args: args,
-                        rep: 0.0,
+                        rep: rep,
-                        active: create_signal(true)
+                        rep_text: create_signal(String::new()),
+                        rep_valid: rep_valid,
+                        active: active,
                     });
                     state.selection.update(|sel| sel.clear());
                     
                     /* DEBUG */
                     // print updated constraint list
-                    console::log_1(&JsValue::from("constraints:"));
+                    console::log_1(&JsValue::from("Constraints:"));
                     state.assembly.constraints.with(|csts| {
                         for (_, cst) in csts.into_iter() {
                             console::log_5(
@@ -227,10 +208,22 @@ pub fn AddRemove() -> View {
                                 &JsValue::from(cst.args.0),
                                 &JsValue::from(cst.args.1),
                                 &JsValue::from(":"),
-                                &JsValue::from(cst.rep)
+                                &JsValue::from(cst.rep.get_untracked())
                             );
                         }
                     });
+                    
+                    // update the realization when the constraint becomes active
+                    // and valid, or is edited while active and valid
+                    create_effect(move || {
+                        console::log_1(&JsValue::from(
+                            format!("Constraint ({}, {}) updated", args.0, args.1)
+                        ));
+                        rep.track();
+                        if active.get() && rep_valid.get() {
+                            state.assembly.realize();
+                        }
+                    });
                 }
             ) { "🔗" }
             select(bind:value=assembly_name) { /* DEBUG */

app-proto/src/assembly.rs

@@ -1,22 +1,49 @@
-use nalgebra::DVector;
+use nalgebra::{DMatrix, DVector};
 use rustc_hash::FxHashMap;
 use slab::Slab;
 use std::collections::BTreeSet;
 use sycamore::prelude::*;
+use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */
+
+use crate::engine::{realize_gram, PartialMatrix};
 
 #[derive(Clone, PartialEq)]
 pub struct Element {
     pub id: String,
     pub label: String,
     pub color: [f32; 3],
-    pub rep: DVector<f64>,
+    pub rep: Signal<DVector<f64>>,
-    pub constraints: BTreeSet<usize>
+    pub constraints: Signal<BTreeSet<usize>>,
+    
+    // internal properties, not reflected in any view
+    pub index: usize
 }
 
+impl Element {
+    pub fn new(
+        id: String,
+        label: String,
+        color: [f32; 3],
+        rep: DVector<f64>
+    ) -> Element {
+        Element {
+            id: id,
+            label: label,
+            color: color,
+            rep: create_signal(rep),
+            constraints: create_signal(BTreeSet::default()),
+            index: 0
+        }
+    }
+}
+            
+
 #[derive(Clone)]
 pub struct Constraint {
     pub args: (usize, usize),
-    pub rep: f64,
+    pub rep: Signal<f64>,
+    pub rep_text: Signal<String>,
+    pub rep_valid: Signal<bool>,
     pub active: Signal<bool>
 }
 
@@ -40,6 +67,8 @@ impl Assembly {
         }
     }
     
+    // --- inserting elements and constraints ---
+    
     // 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
@@ -72,22 +101,103 @@ impl Assembly {
         
         // create and insert a new element
         self.insert_element_unchecked(
-            Element {
+            Element::new(
-                id: id,
+                id,
-                label: format!("Sphere {}", id_num),
+                format!("Sphere {}", id_num),
-                color: [0.75_f32, 0.75_f32, 0.75_f32],
+                [0.75_f32, 0.75_f32, 0.75_f32],
-                rep: DVector::<f64>::from_column_slice(&[0.0, 0.0, 0.0, 0.5, -0.5]),
+                DVector::<f64>::from_column_slice(&[0.0, 0.0, 0.0, 0.5, -0.5])
-                constraints: BTreeSet::default()
+            )
-            }
         );
     }
     
     pub fn insert_constraint(&self, constraint: Constraint) {
         let args = constraint.args;
         let key = self.constraints.update(|csts| csts.insert(constraint));
-        self.elements.update(|elts| {
+        let arg_constraints = self.elements.with(
-            elts[args.0].constraints.insert(key);
+            |elts| (elts[args.0].constraints, elts[args.1].constraints)
-            elts[args.1].constraints.insert(key);
+        );
-        })
+        arg_constraints.0.update(|csts| csts.insert(key));
+        arg_constraints.1.update(|csts| csts.insert(key));
+    }
+    
+    // --- realization ---
+    
+    pub fn realize(&self) {
+        // index the elements
+        self.elements.update_silent(|elts| {
+            for (index, (_, elt)) in elts.into_iter().enumerate() {
+                elt.index = index;
+            }
+        });
+        
+        // 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();
+            self.constraints.with_untracked(|csts| {
+                for (_, cst) in csts {
+                    if cst.active.get_untracked() && cst.rep_valid.get_untracked() {
+                        let args = cst.args;
+                        let row = elts[args.0].index;
+                        let col = elts[args.1].index;
+                        gram_to_be.push_sym(row, col, cst.rep.get_untracked());
+                    }
+                }
+            });
+            
+            // 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.index;
+                gram_to_be.push_sym(index, index, 1.0);
+                guess_to_be.set_column(index, &elt.rep.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();
+        
+        /* DEBUG */
+        // log the initial configuration matrix
+        console::log_1(&JsValue::from("Old configuration:"));
+        for j in 0..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());
+            }
+            console::log_1(&JsValue::from(row_str));
+        }
+        
+        // look for a configuration with the given Gram matrix
+        let (config, success, history) = realize_gram(
+            &gram, guess, &[],
+            1.0e-12, 0.5, 0.9, 1.1, 200, 110
+        );
+        
+        /* DEBUG */
+        // report the outcome of the search
+        console::log_1(&JsValue::from(
+            if success {
+                "Target accuracy achieved!"
+            } else {
+                "Failed to reach target accuracy"
+            }
+        ));
+        console::log_2(&JsValue::from("Steps:"), &JsValue::from(history.scaled_loss.len() - 1));
+        console::log_2(&JsValue::from("Loss:"), &JsValue::from(*history.scaled_loss.last().unwrap()));
+        
+        if success {
+            // read out the solution
+            for (_, elt) in self.elements.get_clone_untracked() {
+                elt.rep.update(
+                    |rep| rep.set_column(0, &config.column(elt.index))
+                );
+            }
+        }
     }
 }

app-proto/src/display.rs

@@ -103,7 +103,11 @@ pub fn Display() -> View {
     // change listener
     let scene_changed = create_signal(true);
     create_effect(move || {
-        state.assembly.elements.track();
+        state.assembly.elements.with(|elts| {
+            for (_, elt) in elts {
+                elt.rep.track();
+            }
+        });
         state.selection.track();
         scene_changed.set(true);
     });
@@ -295,23 +299,40 @@ pub fn Display() -> View {
                 let assembly_to_world = &location * &orientation;
                 
                 // get the assembly
-                let elements = state.assembly.elements.get_clone();
+                let (
-                let element_iter = (&elements).into_iter();
+                    elt_cnt,
-                let reps_world: Vec<_> = element_iter.clone().map(|(_, elt)| &assembly_to_world * &elt.rep).collect();
+                    reps_world,
-                let colors: Vec<_> = element_iter.clone().map(|(key, elt)|
+                    colors,
-                    if state.selection.with(|sel| sel.contains(&key)) {
+                    highlights
-                        elt.color.map(|ch| 0.2 + 0.8*ch)
+                ) = state.assembly.elements.with(|elts| {
-                    } else {
+                    (
-                        elt.color
+                        // number of elements
-                    }
+                        elts.len() as i32,
-                ).collect();
+                        
-                let highlights: Vec<_> = element_iter.map(|(key, _)|
+                        // representation vectors in world coordinates
-                    if state.selection.with(|sel| sel.contains(&key)) {
+                        elts.iter().map(
-                        1.0_f32
+                            |(_, elt)| elt.rep.with(|rep| &assembly_to_world * rep)
-                    } else {
+                        ).collect::<Vec<_>>(),
-                        HIGHLIGHT
+                        
-                    }
+                        // colors
-                ).collect();
+                        elts.iter().map(|(key, elt)| {
+                            if state.selection.with(|sel| sel.contains(&key)) {
+                                elt.color.map(|ch| 0.2 + 0.8*ch)
+                            } else {
+                                elt.color
+                            }
+                        }).collect::<Vec<_>>(),
+                        
+                        // highlight levels
+                        elts.iter().map(|(key, _)| {
+                            if state.selection.with(|sel| sel.contains(&key)) {
+                                1.0_f32
+                            } else {
+                                HIGHLIGHT
+                            }
+                        }).collect::<Vec<_>>()
+                    )
+                });
                 
                 // set the resolution
                 let width = canvas.width() as f32;
@@ -320,7 +341,7 @@ pub fn Display() -> View {
                 ctx.uniform1f(shortdim_loc.as_ref(), width.min(height));
                 
                 // pass the assembly
-                ctx.uniform1i(sphere_cnt_loc.as_ref(), elements.len() as i32);
+                ctx.uniform1i(sphere_cnt_loc.as_ref(), elt_cnt);
                 for n in 0..reps_world.len() {
                     let v = &reps_world[n];
                     ctx.uniform3f(

app-proto/src/engine.rs

@@ -1,4 +1,12 @@
-use nalgebra::DVector;
+use lazy_static::lazy_static;
+use nalgebra::{Const, DMatrix, DVector, Dyn};
+use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */
+
+// --- elements ---
+
+pub fn point(x: f64, y: f64, z: f64) -> DVector<f64> {
+    DVector::from_column_slice(&[x, y, z, 0.5, 0.5*(x*x + y*y + z*z)])
+}
 
 // the sphere with the given center and radius, with inward-pointing normals
 pub fn sphere(center_x: f64, center_y: f64, center_z: f64, radius: f64) -> DVector<f64> {
@@ -24,4 +32,379 @@ pub fn sphere_with_offset(dir_x: f64, dir_y: f64, dir_z: f64, off: f64, curv: f6
         0.5 * curv,
         off * (1.0 + 0.5 * off * curv)
     ])
+}
+
+// --- partial matrices ---
+
+struct MatrixEntry {
+    index: (usize, usize),
+    value: f64
+}
+
+pub struct PartialMatrix(Vec<MatrixEntry>);
+
+impl PartialMatrix {
+    pub fn new() -> PartialMatrix {
+        PartialMatrix(Vec::<MatrixEntry>::new())
+    }
+    
+    pub fn push_sym(&mut self, row: usize, col: usize, value: f64) {
+        let PartialMatrix(entries) = self;
+        entries.push(MatrixEntry { index: (row, col), value: value });
+        if row != col {
+            entries.push(MatrixEntry { index: (col, row), value: 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()));
+        }
+    }
+    
+    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];
+        }
+        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];
+        }
+        result
+    }
+}
+
+// --- descent history ---
+
+pub struct DescentHistory {
+    pub config: Vec<DMatrix<f64>>,
+    pub scaled_loss: Vec<f64>,
+    pub neg_grad: Vec<DMatrix<f64>>,
+    pub min_eigval: Vec<f64>,
+    pub base_step: Vec<DMatrix<f64>>,
+    pub backoff_steps: Vec<i32>
+}
+
+impl DescentHistory {
+    fn new() -> DescentHistory {
+        DescentHistory {
+            config: Vec::<DMatrix<f64>>::new(),
+            scaled_loss: Vec::<f64>::new(),
+            neg_grad: Vec::<DMatrix<f64>>::new(),
+            min_eigval: Vec::<f64>::new(),
+            base_step: Vec::<DMatrix<f64>>::new(),
+            backoff_steps: Vec::<i32>::new(),
+        }
+    }
+}
+
+// --- gram matrix realization ---
+
+// the Lorentz form
+lazy_static! {
+    pub static ref Q: DMatrix<f64> = DMatrix::from_row_slice(5, 5, &[
+        1.0, 0.0, 0.0,  0.0,  0.0,
+        0.0, 1.0, 0.0,  0.0,  0.0,
+        0.0, 0.0, 1.0,  0.0,  0.0,
+        0.0, 0.0, 0.0,  0.0, -2.0,
+        0.0, 0.0, 0.0, -2.0,  0.0
+    ]);
+}
+
+struct SearchState {
+    config: DMatrix<f64>,
+    err_proj: DMatrix<f64>,
+    loss: f64
+}
+
+impl SearchState {
+    fn from_config(gram: &PartialMatrix, config: DMatrix<f64>) -> SearchState {
+        let err_proj = gram.sub_proj(&(config.tr_mul(&*Q) * &config));
+        let loss = err_proj.norm_squared();
+        SearchState {
+            config: config,
+            err_proj: err_proj,
+            loss: loss
+        }
+    }
+}
+
+fn basis_matrix(index: (usize, usize), nrows: usize, ncols: usize) -> DMatrix<f64> {
+    let mut result = DMatrix::<f64>::zeros(nrows, ncols);
+    result[index] = 1.0;
+    result
+}
+
+// use backtracking line search to find a better configuration
+fn seek_better_config(
+    gram: &PartialMatrix,
+    state: &SearchState,
+    base_step: &DMatrix<f64>,
+    base_target_improvement: f64,
+    min_efficiency: f64,
+    backoff: f64,
+    max_backoff_steps: i32
+) -> Option<(SearchState, i32)> {
+    let mut rate = 1.0;
+    for backoff_steps in 0..max_backoff_steps {
+        let trial_config = &state.config + rate * base_step;
+        let trial_state = SearchState::from_config(gram, trial_config);
+        let improvement = state.loss - trial_state.loss;
+        if improvement >= min_efficiency * rate * base_target_improvement {
+            return Some((trial_state, backoff_steps));
+        }
+        rate *= backoff;
+    }
+    None
+}
+
+// seek a matrix `config` for which `config' * Q * config` matches the partial
+// matrix `gram`. use gradient descent starting from `guess`
+pub fn realize_gram(
+    gram: &PartialMatrix,
+    guess: DMatrix<f64>,
+    frozen: &[(usize, usize)],
+    scaled_tol: f64,
+    min_efficiency: f64,
+    backoff: f64,
+    reg_scale: f64,
+    max_descent_steps: i32,
+    max_backoff_steps: i32
+) -> (DMatrix<f64>, bool, DescentHistory) {
+    // start the descent history
+    let mut history = DescentHistory::new();
+    
+    // find the dimension of the search space
+    let element_dim = guess.nrows();
+    let assembly_dim = guess.ncols();
+    let total_dim = element_dim * assembly_dim;
+    
+    // scale the tolerance
+    let scale_adjustment = (gram.0.len() as f64).sqrt();
+    let tol = scale_adjustment * scaled_tol;
+    
+    // convert the frozen indices to stacked format
+    let frozen_stacked: Vec<usize> = frozen.into_iter().map(
+        |index| index.1*element_dim + index.0
+    ).collect();
+    
+    // use Newton's method with backtracking and gradient descent backup
+    let mut state = SearchState::from_config(gram, guess);
+    for _ in 0..max_descent_steps {
+        // stop if the loss is tolerably low
+        history.config.push(state.config.clone());
+        history.scaled_loss.push(state.loss / scale_adjustment);
+        if state.loss < tol { break; }
+        
+        // find the negative gradient of the loss function
+        let neg_grad = 4.0 * &*Q * &state.config * &state.err_proj;
+        let mut neg_grad_stacked = neg_grad.clone().reshape_generic(Dyn(total_dim), Const::<1>);
+        history.neg_grad.push(neg_grad.clone());
+        
+        // find the negative Hessian of the loss function
+        let mut hess_cols = Vec::<DVector<f64>>::with_capacity(total_dim);
+        for col in 0..assembly_dim {
+            for row in 0..element_dim {
+                let index = (row, col);
+                let basis_mat = basis_matrix(index, element_dim, assembly_dim);
+                let neg_d_err =
+                    basis_mat.tr_mul(&*Q) * &state.config
+                    + state.config.tr_mul(&*Q) * &basis_mat;
+                let neg_d_err_proj = gram.proj(&neg_d_err);
+                let deriv_grad = 4.0 * &*Q * (
+                    -&basis_mat * &state.err_proj
+                    + &state.config * &neg_d_err_proj
+                );
+                hess_cols.push(deriv_grad.reshape_generic(Dyn(total_dim), Const::<1>));
+            }
+        }
+        let mut hess = DMatrix::from_columns(hess_cols.as_slice());
+        
+        // regularize the Hessian
+        let min_eigval = hess.symmetric_eigenvalues().min();
+        if min_eigval <= 0.0 {
+            hess -= reg_scale * min_eigval * DMatrix::identity(total_dim, total_dim);
+        }
+        history.min_eigval.push(min_eigval);
+        
+        // project the negative gradient and negative Hessian onto the
+        // orthogonal complement of the frozen subspace
+        let zero_col = DVector::zeros(total_dim);
+        let zero_row = zero_col.transpose();
+        for &k in &frozen_stacked {
+            neg_grad_stacked[k] = 0.0;
+            hess.set_row(k, &zero_row);
+            hess.set_column(k, &zero_col);
+            hess[(k, k)] = 1.0;
+        }
+        
+        // compute the Newton step
+        /*
+          we need to either handle or eliminate the case where the minimum
+          eigenvalue of the Hessian is zero, so the regularized Hessian is
+          singular. right now, this causes the Cholesky decomposition to return
+          `None`, leading to a panic when we unrap
+        */
+        let base_step_stacked = hess.cholesky().unwrap().solve(&neg_grad_stacked);
+        let base_step = base_step_stacked.reshape_generic(Dyn(element_dim), Dyn(assembly_dim));
+        history.base_step.push(base_step.clone());
+        
+        // use backtracking line search to find a better configuration
+        match seek_better_config(
+            gram, &state, &base_step, neg_grad.dot(&base_step),
+            min_efficiency, backoff, max_backoff_steps
+        ) {
+            Some((better_state, backoff_steps)) => {
+                state = better_state;
+                history.backoff_steps.push(backoff_steps);
+            },
+            None => return (state.config, false, history)
+        };
+    }
+    (state.config, state.loss < tol, history)
+}
+
+// --- tests ---
+
+// this problem is from a sangaku by Irisawa Shintarō Hiroatsu. the article
+// below includes a nice translation of the problem statement, which was
+// recorded in Uchida Itsumi's book _Kokon sankan_ (_Mathematics, Past and
+// Present_)
+//
+//   "Japan's 'Wasan' Mathematical Tradition", by Abe Haruki
+//   https://www.nippon.com/en/japan-topics/c12801/
+//
+#[cfg(feature = "irisawa")]
+pub mod irisawa {
+    use std::{array, f64::consts::PI};
+    
+    use super::*;
+    
+    pub fn realize_irisawa_hexlet(scaled_tol: f64) -> (DMatrix<f64>, 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(
+            [
+                sphere(0.0, 0.0, 0.0, 15.0),
+                sphere(0.0, 0.0, -9.0, 5.0),
+                sphere(0.0, 0.0, 11.0, 3.0)
+            ].into_iter().chain(
+                (1..=6).map(
+                    |k| {
+                        let ang = (k as f64) * PI/3.0;
+                        sphere(9.0 * ang.cos(), 9.0 * ang.sin(), 0.0, 2.5)
+                    }
+                )
+            ).collect::<Vec<_>>().as_slice()
+        );
+        
+        // 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));
+        
+        realize_gram(
+            &gram, guess, &frozen,
+            scaled_tol, 0.5, 0.9, 1.1, 200, 110
+        )
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::{*, irisawa::realize_irisawa_hexlet};
+    
+    #[test]
+    fn sub_proj_test() {
+        let target = PartialMatrix(vec![
+            MatrixEntry { index: (0, 0), value: 19.0 },
+            MatrixEntry { index: (0, 2), value: 39.0 },
+            MatrixEntry { index: (1, 1), value: 59.0 },
+            MatrixEntry { index: (1, 2), value: 69.0 }
+        ]);
+        let attempt = 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, &[
+            18.0, 0.0, 36.0,
+            0.0, 54.0, 63.0
+        ]);
+        assert_eq!(target.sub_proj(&attempt), expected_result);
+    }
+    
+    #[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 }
+                    });
+                }
+            }
+            entries
+        });
+        let config = {
+            let a: f64 = 0.75_f64.sqrt();
+            DMatrix::from_columns(&[
+                sphere(1.0, 0.0, 0.0, a),
+                sphere(-0.5, a, 0.0, a),
+                sphere(-0.5, -a, 0.0, a)
+            ])
+        };
+        let state = SearchState::from_config(&gram, config);
+        assert!(state.loss.abs() < f64::EPSILON);
+    }
+    
+    #[test]
+    fn irisawa_hexlet_test() {
+        // solve Irisawa's problem
+        const SCALED_TOL: f64 = 1.0e-12;
+        let (config, _, _) = realize_irisawa_hexlet(SCALED_TOL);
+        
+        // check against Irisawa's solution
+        let entry_tol = SCALED_TOL.sqrt();
+        let solution_diams = [30.0, 10.0, 6.0, 5.0, 15.0, 10.0, 3.75, 2.5, 2.0 + 8.0/11.0];
+        for (k, diam) in solution_diams.into_iter().enumerate() {
+            assert!((config[(3, k)] - 1.0 / diam).abs() < entry_tol);
+        }
+    }
 }

app-proto/src/lib.rs

@@ -0,0 +1 @@
+pub mod engine;

app-proto/src/outline.rs

@@ -1,26 +1,191 @@
 use itertools::Itertools;
-use sycamore::{prelude::*, web::tags::div};
+use sycamore::prelude::*;
-use web_sys::{Element, KeyboardEvent, MouseEvent, wasm_bindgen::JsCast};
+use web_sys::{
+    Event,
+    HtmlInputElement,
+    KeyboardEvent,
+    MouseEvent,
+    wasm_bindgen::JsCast
+};
 
-use crate::AppState;
+use crate::{AppState, assembly, assembly::Constraint};
 
-// this component lists the elements of the assembly, showing the constraints
+// an editable view of the Lorentz product representing a constraint
-// on each element as a collapsible sub-list. its implementation is based on
+#[component(inline_props)]
-// Kate Morley's HTML + CSS tree views:
+fn LorentzProductInput(constraint: Constraint) -> View {
+    view! {
+        input(
+            r#type="text",
+            bind:value=constraint.rep_text,
+            on:change=move |event: Event| {
+                let target: HtmlInputElement = event.target().unwrap().unchecked_into();
+                match target.value().parse::<f64>() {
+                    Ok(rep) => batch(|| {
+                        constraint.rep.set(rep);
+                        constraint.rep_valid.set(true);
+                    }),
+                    Err(_) => constraint.rep_valid.set(false)
+                };
+            }
+        )
+    }
+}
+
+// a list item that shows a constraint in an outline view of an element
+#[component(inline_props)]
+fn ConstraintOutlineItem(constraint_key: usize, element_key: usize) -> View {
+    let state = use_context::<AppState>();
+    let assembly = &state.assembly;
+    let constraint = assembly.constraints.with(|csts| csts[constraint_key].clone());
+    let other_arg = if constraint.args.0 == element_key {
+        constraint.args.1
+    } else {
+        constraint.args.0
+    };
+    let other_arg_label = assembly.elements.with(|elts| elts[other_arg].label.clone());
+    let class = constraint.rep_valid.map(
+        |&rep_valid| if rep_valid { "cst" } else { "cst invalid" }
+    );
+    view! {
+        li(class=class.get()) {
+            input(r#type="checkbox", bind:checked=constraint.active)
+            div(class="cst-label") { (other_arg_label) }
+            LorentzProductInput(constraint=constraint)
+            div(class="status")
+        }
+    }
+}
+
+// a list item that shows an element in an outline view of an assembly
+#[component(inline_props)]
+fn ElementOutlineItem(key: usize, element: assembly::Element) -> View {
+    let state = use_context::<AppState>();
+    let class = state.selection.map(
+        move |sel| if sel.contains(&key) { "selected" } else { "" }
+    );
+    let label = element.label.clone();
+    let rep_components = element.rep.map(
+        |rep| rep.iter().map(
+            |u| format!("{:.3}", u).replace("-", "\u{2212}")
+        ).collect()
+    );
+    let constrained = element.constraints.map(|csts| csts.len() > 0);
+    let constraint_list = element.constraints.map(
+        |csts| csts.clone().into_iter().collect()
+    );
+    let details_node = create_node_ref();
+    view! {
+        li {
+            details(ref=details_node) {
+                summary(
+                    class=class.get(),
+                    on:keydown={
+                        move |event: KeyboardEvent| {
+                            match event.key().as_str() {
+                                "Enter" => {
+                                    if event.shift_key() {
+                                        state.selection.update(|sel| {
+                                            if !sel.remove(&key) {
+                                                sel.insert(key);
+                                            }
+                                        });
+                                    } else {
+                                        state.selection.update(|sel| {
+                                            sel.clear();
+                                            sel.insert(key);
+                                        });
+                                    }
+                                    event.prevent_default();
+                                },
+                                "ArrowRight" if constrained.get() => {
+                                    let _ = details_node
+                                        .get()
+                                        .unchecked_into::<web_sys::Element>()
+                                        .set_attribute("open", "");
+                                },
+                                "ArrowLeft" => {
+                                    let _ = details_node
+                                        .get()
+                                        .unchecked_into::<web_sys::Element>()
+                                        .remove_attribute("open");
+                                },
+                                _ => ()
+                            }
+                        }
+                    }
+                ) {
+                    div(
+                        class="elt-switch",
+                        on:click=|event: MouseEvent| event.stop_propagation()
+                    )
+                    div(
+                        class="elt",
+                        on:click={
+                            move |event: MouseEvent| {
+                                if event.shift_key() {
+                                    state.selection.update(|sel| {
+                                        if !sel.remove(&key) {
+                                            sel.insert(key);
+                                        }
+                                    });
+                                } else {
+                                    state.selection.update(|sel| {
+                                        sel.clear();
+                                        sel.insert(key);
+                                    });
+                                }
+                                event.stop_propagation();
+                                event.prevent_default();
+                            }
+                        }
+                    ) {
+                        div(class="elt-label") { (label) }
+                        div(class="elt-rep") {
+                            Indexed(
+                                list=rep_components,
+                                view=|coord_str| view! {
+                                    div { (coord_str) }
+                                }
+                            )
+                        }
+                        div(class="status")
+                    }
+                }
+                ul(class="constraints") {
+                    Keyed(
+                        list=constraint_list,
+                        view=move |cst_key| view! {
+                            ConstraintOutlineItem(
+                                constraint_key=cst_key,
+                                element_key=key
+                            )
+                        },
+                        key=|cst_key| cst_key.clone()
+                    )
+                }
+            }
+        }
+    }
+}
+
+// a component that lists the elements of the current assembly, showing the
+// constraints on each element as a collapsible sub-list. its implementation
+// is based on Kate Morley's HTML + CSS tree views:
 //
 //   https://iamkate.com/code/tree-views/
 //
 #[component]
 pub fn Outline() -> View {
-    // sort the elements alphabetically by ID
+    let state = use_context::<AppState>();
-    let elements_sorted = create_memo(|| {
+    
-        let state = use_context::<AppState>();
+    // list the elements alphabetically by ID
-        state.assembly.elements
+    let element_list = state.assembly.elements.map(
-            .get_clone()
+        |elts| elts
+            .clone()
             .into_iter()
             .sorted_by_key(|(_, elt)| elt.id.clone())
             .collect()
-    });
+    );
     
     view! {
         ul(
@@ -31,130 +196,11 @@ pub fn Outline() -> View {
             }
         ) {
             Keyed(
-                list=elements_sorted,
+                list=element_list,
-                view=|(key, elt)| {
+                view=|(key, elt)| view! {
-                    let state = use_context::<AppState>();
+                    ElementOutlineItem(key=key, element=elt)
-                    let class = create_memo({
-                        move || {
-                            if state.selection.with(|sel| sel.contains(&key)) {
-                                "selected"
-                            } else {
-                                ""
-                            }
-                        }
-                    });
-                    let label = elt.label.clone();
-                    let rep_components = elt.rep.iter().map(|u| {
-                        let u_coord = u.to_string().replace("-", "\u{2212}");
-                        View::from(div().children(u_coord))
-                    }).collect::<Vec<_>>();
-                    let constrained = elt.constraints.len() > 0;
-                    let details_node = create_node_ref();
-                    view! {
-                        /* [TO DO] switch to integer-valued parameters whenever
-                        that becomes possible again */
-                        li {
-                            details(ref=details_node) {
-                                summary(
-                                    class=class.get(),
-                                    on:keydown={
-                                        move |event: KeyboardEvent| {
-                                            match event.key().as_str() {
-                                                "Enter" => {
-                                                    if event.shift_key() {
-                                                        state.selection.update(|sel| {
-                                                            if !sel.remove(&key) {
-                                                                sel.insert(key);
-                                                            }
-                                                        });
-                                                    } else {
-                                                        state.selection.update(|sel| {
-                                                            sel.clear();
-                                                            sel.insert(key);
-                                                        });
-                                                    }
-                                                    event.prevent_default();
-                                                },
-                                                "ArrowRight" if constrained => {
-                                                    let _ = details_node
-                                                        .get()
-                                                        .unchecked_into::<Element>()
-                                                        .set_attribute("open", "");
-                                                },
-                                                "ArrowLeft" => {
-                                                    let _ = details_node
-                                                        .get()
-                                                        .unchecked_into::<Element>()
-                                                        .remove_attribute("open");
-                                                },
-                                                _ => ()
-                                            }
-                                        }
-                                    }
-                                ) {
-                                    div(
-                                        class="elt-switch",
-                                        on:click=|event: MouseEvent| event.stop_propagation()
-                                    )
-                                    div(
-                                        class="elt",
-                                        on:click={
-                                            move |event: MouseEvent| {
-                                                if event.shift_key() {
-                                                    state.selection.update(|sel| {
-                                                        if !sel.remove(&key) {
-                                                            sel.insert(key);
-                                                        }
-                                                    });
-                                                } else {
-                                                    state.selection.update(|sel| {
-                                                        sel.clear();
-                                                        sel.insert(key);
-                                                    });
-                                                }
-                                                event.stop_propagation();
-                                                event.prevent_default();
-                                            }
-                                        }
-                                    ) {
-                                        div(class="elt-label") { (label) }
-                                        div(class="elt-rep") { (rep_components) }
-                                    }
-                                }
-                                ul(class="constraints") {
-                                    Keyed(
-                                        list=elt.constraints.into_iter().collect::<Vec<_>>(),
-                                        view=move |c_key: usize| {
-                                            let c_state = use_context::<AppState>();
-                                            let assembly = &c_state.assembly;
-                                            let cst = assembly.constraints.with(|csts| csts[c_key].clone());
-                                            let other_arg = if cst.args.0 == key {
-                                                cst.args.1
-                                            } else {
-                                                cst.args.0
-                                            };
-                                            let other_arg_label = assembly.elements.with(|elts| elts[other_arg].label.clone());
-                                            view! {
-                                                li(class="cst") {
-                                                    input(r#type="checkbox", bind:checked=cst.active)
-                                                    div(class="cst-label") { (other_arg_label) }
-                                                    div(class="cst-rep") { (cst.rep) }
-                                                }
-                                            }
-                                        },
-                                        key=|c_key| c_key.clone()
-                                    )
-                                }
-                            }
-                        }
-                    }
                 },
-                key=|(key, elt)| (
+                key=|(key, _)| key.clone()
-                    key.clone(),
-                    elt.id.clone(),
-                    elt.label.clone(),
-                    elt.constraints.clone()
-                )
             )
         }
     }

engine-proto/gram-test/Engine.jl

@@ -8,7 +8,8 @@ using Optim
 
 export
   rand_on_shell, Q, DescentHistory,
-  realize_gram_gradient, realize_gram_newton, realize_gram_optim, realize_gram
+  realize_gram_gradient, realize_gram_newton, realize_gram_optim,
+  realize_gram_alt_proj, realize_gram
 
 # === guessing ===
 
@@ -143,7 +144,7 @@ function realize_gram_gradient(
       break
     end
     
-    # find negative gradient of loss function
+    # find the negative gradient of the loss function
     neg_grad = 4*Q*L*Δ_proj
     slope = norm(neg_grad)
     dir = neg_grad / slope
@@ -232,7 +233,7 @@ function realize_gram_newton(
       break
     end
     
-    # find the negative gradient of loss function
+    # find the negative gradient of the loss function
     neg_grad = 4*Q*L*Δ_proj
     
     # find the negative Hessian of the loss function
@@ -313,6 +314,129 @@ function realize_gram_optim(
   )
 end
 
+# seek a matrix `L` for which `L'QL` matches the sparse matrix `gram` at every
+# explicit entry of `gram`. use gradient descent starting from `guess`, with an
+# alternate technique for finding the projected base step from the unprojected
+# Hessian
+function realize_gram_alt_proj(
+  gram::SparseMatrixCSC{T, <:Any},
+  guess::Matrix{T},
+  frozen = CartesianIndex[];
+  scaled_tol = 1e-30,
+  min_efficiency = 0.5,
+  backoff = 0.9,
+  reg_scale = 1.1,
+  max_descent_steps = 200,
+  max_backoff_steps = 110
+) where T <: Number
+  # start history
+  history = DescentHistory{T}()
+  
+  # find the dimension of the search space
+  dims = size(guess)
+  element_dim, construction_dim = dims
+  total_dim = element_dim * construction_dim
+  
+  # list the constrained entries of the gram matrix
+  J, K, _ = findnz(gram)
+  constrained = zip(J, K)
+  
+  # scale the tolerance
+  scale_adjustment = sqrt(T(length(constrained)))
+  tol = scale_adjustment * scaled_tol
+  
+  # convert the frozen indices to stacked format
+  frozen_stacked = [(index[2]-1)*element_dim + index[1] for index in frozen]
+  
+  # initialize search state
+  L = copy(guess)
+  Δ_proj = proj_diff(gram, L'*Q*L)
+  loss = dot(Δ_proj, Δ_proj)
+  
+  # use Newton's method with backtracking and gradient descent backup
+  for step in 1:max_descent_steps
+    # stop if the loss is tolerably low
+    if loss < tol
+      break
+    end
+    
+    # find the negative gradient of the loss function
+    neg_grad = 4*Q*L*Δ_proj
+    
+    # find the negative Hessian of the loss function
+    hess = Matrix{T}(undef, total_dim, total_dim)
+    indices = [(j, k) for k in 1:construction_dim for j in 1:element_dim]
+    for (j, k) in indices
+      basis_mat = basis_matrix(T, j, k, dims)
+      neg_dΔ = basis_mat'*Q*L + L'*Q*basis_mat
+      neg_dΔ_proj = proj_to_entries(neg_dΔ, constrained)
+      deriv_grad = 4*Q*(-basis_mat*Δ_proj + L*neg_dΔ_proj)
+      hess[:, (k-1)*element_dim + j] = reshape(deriv_grad, total_dim)
+    end
+    hess_sym = Hermitian(hess)
+    push!(history.hess, hess_sym)
+    
+    # regularize the Hessian
+    min_eigval = minimum(eigvals(hess_sym))
+    push!(history.positive, min_eigval > 0)
+    if min_eigval <= 0
+      hess -= reg_scale * min_eigval * I
+    end
+    
+    # compute the Newton step
+    neg_grad_stacked = reshape(neg_grad, total_dim)
+    for k in frozen_stacked
+      neg_grad_stacked[k] = 0
+      hess[k, :] .= 0
+      hess[:, k] .= 0
+      hess[k, k] = 1
+    end
+    base_step_stacked = Hermitian(hess) \ neg_grad_stacked
+    base_step = reshape(base_step_stacked, dims)
+    push!(history.base_step, base_step)
+    
+    # store the current position, loss, and slope
+    L_last = L
+    loss_last = loss
+    push!(history.scaled_loss, loss / scale_adjustment)
+    push!(history.neg_grad, neg_grad)
+    push!(history.slope, norm(neg_grad))
+    
+    # find a good step size using backtracking line search
+    push!(history.stepsize, 0)
+    push!(history.backoff_steps, max_backoff_steps)
+    empty!(history.last_line_L)
+    empty!(history.last_line_loss)
+    rate = one(T)
+    step_success = false
+    base_target_improvement = dot(neg_grad, base_step)
+    for backoff_steps in 0:max_backoff_steps
+      history.stepsize[end] = rate
+      L = L_last + rate * base_step
+      Δ_proj = proj_diff(gram, L'*Q*L)
+      loss = dot(Δ_proj, Δ_proj)
+      improvement = loss_last - loss
+      push!(history.last_line_L, L)
+      push!(history.last_line_loss, loss / scale_adjustment)
+      if improvement >= min_efficiency * rate * base_target_improvement
+        history.backoff_steps[end] = backoff_steps
+        step_success = true
+        break
+      end
+      rate *= backoff
+    end
+    
+    # if we've hit a wall, quit
+    if !step_success
+      return L_last, false, history
+    end
+  end
+  
+  # return the factorization and its history
+  push!(history.scaled_loss, loss / scale_adjustment)
+  L, loss < tol, history
+end
+
 # seek a matrix `L` for which `L'QL` matches the sparse matrix `gram` at every
 # explicit entry of `gram`. use gradient descent starting from `guess`
 function realize_gram(
@@ -321,7 +445,6 @@ function realize_gram(
   frozen = nothing;
   scaled_tol = 1e-30,
   min_efficiency = 0.5,
-  init_rate = 1.0,
   backoff = 0.9,
   reg_scale = 1.1,
   max_descent_steps = 200,
@@ -352,20 +475,19 @@ function realize_gram(
     unfrozen_stacked = reshape(is_unfrozen, total_dim)
   end
   
-  # initialize variables
+  # initialize search state
-  grad_rate = init_rate
   L = copy(guess)
-  
-  # use Newton's method with backtracking and gradient descent backup
   Δ_proj = proj_diff(gram, L'*Q*L)
   loss = dot(Δ_proj, Δ_proj)
+  
+  # use Newton's method with backtracking and gradient descent backup
   for step in 1:max_descent_steps
     # stop if the loss is tolerably low
     if loss < tol
       break
     end
     
-    # find the negative gradient of loss function
+    # find the negative gradient of the loss function
     neg_grad = 4*Q*L*Δ_proj
     
     # find the negative Hessian of the loss function
@@ -420,6 +542,7 @@ function realize_gram(
     empty!(history.last_line_loss)
     rate = one(T)
     step_success = false
+    base_target_improvement = dot(neg_grad, base_step)
     for backoff_steps in 0:max_backoff_steps
       history.stepsize[end] = rate
       L = L_last + rate * base_step
@@ -428,7 +551,7 @@ function realize_gram(
       improvement = loss_last - loss
       push!(history.last_line_L, L)
       push!(history.last_line_loss, loss / scale_adjustment)
-      if improvement >= min_efficiency * rate * dot(neg_grad, base_step)
+      if improvement >= min_efficiency * rate * base_target_improvement
         history.backoff_steps[end] = backoff_steps
         step_success = true
         break

engine-proto/gram-test/irisawa-hexlet.jl

@@ -74,4 +74,13 @@ if success
   for k in 5:9
     println("  ", 1 / L[4,k], " sun")
   end
-end
+end
+
+# test an alternate technique for finding the projected base step from the
+# unprojected Hessian
+L_alt, success_alt, history_alt = Engine.realize_gram_alt_proj(gram, guess, frozen)
+completed_gram_alt = L_alt'*Engine.Q*L_alt
+println("\nDifference in result using alternate projection:\n")
+display(completed_gram_alt - completed_gram)
+println("\nDifference in steps: ", size(history_alt.scaled_loss, 1) - size(history.scaled_loss, 1))
+println("Difference in loss: ", history_alt.scaled_loss[end] - history.scaled_loss[end], "\n")

engine-proto/gram-test/sphere-in-tetrahedron.jl

@@ -64,4 +64,13 @@ else
   println("\nFailed to reach target accuracy")
 end
 println("Steps: ", size(history.scaled_loss, 1))
-println("Loss: ", history.scaled_loss[end], "\n")
+println("Loss: ", history.scaled_loss[end], "\n")
+
+# test an alternate technique for finding the projected base step from the
+# unprojected Hessian
+L_alt, success_alt, history_alt = Engine.realize_gram_alt_proj(gram, guess, frozen)
+completed_gram_alt = L_alt'*Engine.Q*L_alt
+println("\nDifference in result using alternate projection:\n")
+display(completed_gram_alt - completed_gram)
+println("\nDifference in steps: ", size(history_alt.scaled_loss, 1) - size(history.scaled_loss, 1))
+println("Difference in loss: ", history_alt.scaled_loss[end] - history.scaled_loss[end], "\n")

engine-proto/gram-test/tetrahedron-radius-ratio.jl

@@ -93,4 +93,13 @@ if success
   infty = BigFloat[0, 0, 0, 0, 1]
   radius_ratio = dot(infty, Engine.Q * L[:,5]) / dot(infty, Engine.Q * L[:,6])
   println("\nCircumradius / inradius: ", radius_ratio)
-end
+end
+
+# test an alternate technique for finding the projected base step from the
+# unprojected Hessian
+L_alt, success_alt, history_alt = Engine.realize_gram_alt_proj(gram, guess, frozen)
+completed_gram_alt = L_alt'*Engine.Q*L_alt
+println("\nDifference in result using alternate projection:\n")
+display(completed_gram_alt - completed_gram)
+println("\nDifference in steps: ", size(history_alt.scaled_loss, 1) - size(history.scaled_loss, 1))
+println("Difference in loss: ", history_alt.scaled_loss[end] - history.scaled_loss[end], "\n")