Enforce constraints in the editor

app-proto/src/add_remove.rs

@@ -12,7 +12,8 @@ fn load_gen_assemb(assembly: &Assembly) {
             label: String::from("Castor"),
             color: [1.00_f32, 0.25_f32, 0.00_f32],
             rep: engine::sphere(0.5, 0.5, 0.0, 1.0),
-            constraints: BTreeSet::default()
+            constraints: BTreeSet::default(),
+            index: 0
         }
     );
     let _ = assembly.try_insert_element(
@@ -21,7 +22,8 @@ fn load_gen_assemb(assembly: &Assembly) {
             label: String::from("Pollux"),
             color: [0.00_f32, 0.25_f32, 1.00_f32],
             rep: engine::sphere(-0.5, -0.5, 0.0, 1.0),
-            constraints: BTreeSet::default()
+            constraints: BTreeSet::default(),
+            index: 0
         }
     );
     let _ = assembly.try_insert_element(
@@ -30,7 +32,8 @@ fn load_gen_assemb(assembly: &Assembly) {
             label: String::from("Ursa major"),
             color: [0.25_f32, 0.00_f32, 1.00_f32],
             rep: engine::sphere(-0.5, 0.5, 0.0, 0.75),
-            constraints: BTreeSet::default()
+            constraints: BTreeSet::default(),
+            index: 0
         }
     );
     let _ = assembly.try_insert_element(
@@ -39,7 +42,8 @@ fn load_gen_assemb(assembly: &Assembly) {
             label: String::from("Ursa minor"),
             color: [0.25_f32, 1.00_f32, 0.00_f32],
             rep: engine::sphere(0.5, -0.5, 0.0, 0.5),
-            constraints: BTreeSet::default()
+            constraints: BTreeSet::default(),
+            index: 0
         }
     );
     let _ = assembly.try_insert_element(
@@ -48,7 +52,8 @@ fn load_gen_assemb(assembly: &Assembly) {
             label: String::from("Deimos"),
             color: [0.75_f32, 0.75_f32, 0.00_f32],
             rep: engine::sphere(0.0, 0.15, 1.0, 0.25),
-            constraints: BTreeSet::default()
+            constraints: BTreeSet::default(),
+            index: 0
         }
     );
     let _ = assembly.try_insert_element(
@@ -57,17 +62,8 @@ fn load_gen_assemb(assembly: &Assembly) {
             label: String::from("Phobos"),
             color: [0.00_f32, 0.75_f32, 0.50_f32],
             rep: engine::sphere(0.0, -0.15, -1.0, 0.25),
-            constraints: BTreeSet::default()
+            constraints: BTreeSet::default(),
-        }
+            index: 0
-    );
-    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)
         }
     );
 }
@@ -81,7 +77,8 @@ fn load_low_curv_assemb(assembly: &Assembly) {
             label: "Central".to_string(),
             color: [0.75_f32, 0.75_f32, 0.75_f32],
             rep: engine::sphere(0.0, 0.0, 0.0, 1.0),
-            constraints: BTreeSet::default()
+            constraints: BTreeSet::default(),
+            index: 0
         }
     );
     let _ = assembly.try_insert_element(
@@ -90,7 +87,8 @@ fn load_low_curv_assemb(assembly: &Assembly) {
             label: "Assembly plane".to_string(),
             color: [0.75_f32, 0.75_f32, 0.75_f32],
             rep: engine::sphere_with_offset(0.0, 0.0, 1.0, 0.0, 0.0),
-            constraints: BTreeSet::default()
+            constraints: BTreeSet::default(),
+            index: 0
         }
     );
     let _ = assembly.try_insert_element(
@@ -99,7 +97,8 @@ fn load_low_curv_assemb(assembly: &Assembly) {
             label: "Side 1".to_string(),
             color: [1.00_f32, 0.00_f32, 0.25_f32],
             rep: engine::sphere_with_offset(1.0, 0.0, 0.0, 1.0, 0.0),
-            constraints: BTreeSet::default()
+            constraints: BTreeSet::default(),
+            index: 0
         }
     );
     let _ = assembly.try_insert_element(
@@ -108,7 +107,8 @@ fn load_low_curv_assemb(assembly: &Assembly) {
             label: "Side 2".to_string(),
             color: [0.25_f32, 1.00_f32, 0.00_f32],
             rep: engine::sphere_with_offset(-0.5, a, 0.0, 1.0, 0.0),
-            constraints: BTreeSet::default()
+            constraints: BTreeSet::default(),
+            index: 0
         }
     );
     let _ = assembly.try_insert_element(
@@ -117,7 +117,8 @@ fn load_low_curv_assemb(assembly: &Assembly) {
             label: "Side 3".to_string(),
             color: [0.00_f32, 0.25_f32, 1.00_f32],
             rep: engine::sphere_with_offset(-0.5, -a, 0.0, 1.0, 0.0),
-            constraints: BTreeSet::default()
+            constraints: BTreeSet::default(),
+            index: 0
         }
     );
     let _ = assembly.try_insert_element(
@@ -126,7 +127,8 @@ fn load_low_curv_assemb(assembly: &Assembly) {
             label: "Corner 1".to_string(),
             color: [0.75_f32, 0.75_f32, 0.75_f32],
             rep: engine::sphere(-4.0/3.0, 0.0, 0.0, 1.0/3.0),
-            constraints: BTreeSet::default()
+            constraints: BTreeSet::default(),
+            index: 0
         }
     );
     let _ = assembly.try_insert_element(
@@ -135,7 +137,8 @@ fn load_low_curv_assemb(assembly: &Assembly) {
             label: "Corner 2".to_string(),
             color: [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),
-            constraints: BTreeSet::default()
+            constraints: BTreeSet::default(),
+            index: 0
         }
     );
     let _ = assembly.try_insert_element(
@@ -144,7 +147,8 @@ fn load_low_curv_assemb(assembly: &Assembly) {
             label: String::from("Corner 3"),
             color: [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),
-            constraints: BTreeSet::default()
+            constraints: BTreeSet::default(),
+            index: 0
         }
     );
 }
@@ -215,6 +219,7 @@ pub fn AddRemove() -> View {
                         rep: 0.0,
                         active: create_signal(true)
                     });
+                    state.assembly.realize();
                     state.selection.update(|sel| sel.clear());
                     
                     /* DEBUG */

app-proto/src/assembly.rs

@@ -1,8 +1,11 @@
-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 {
@@ -10,7 +13,10 @@ pub struct Element {
     pub label: String,
     pub color: [f32; 3],
     pub rep: DVector<f64>,
-    pub constraints: BTreeSet<usize>
+    pub constraints: BTreeSet<usize>,
+    
+    // internal properties, not reflected in any view
+    pub index: usize
 }
 
 #[derive(Clone)]
@@ -40,6 +46,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
@@ -77,7 +85,8 @@ impl Assembly {
                 label: format!("Sphere {}", id_num),
                 color: [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]),
-                constraints: BTreeSet::default()
+                constraints: BTreeSet::default(),
+                index: 0
             }
         );
     }
@@ -90,4 +99,83 @@ impl Assembly {
             elts[args.1].constraints.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 {
+                    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);
+                }
+            });
+            
+            // 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);
+            }
+            
+            (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
+            self.elements.update(|elts| {
+                for (_, elt) in elts.iter_mut() {
+                    elt.rep.set_column(0, &config.column(elt.index));
+                }
+            });
+        }
+    }
 }

app-proto/src/engine.rs

@@ -1,5 +1,6 @@
 use lazy_static::lazy_static;
 use nalgebra::{Const, DMatrix, DVector, Dyn};
+use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */
 
 // --- elements ---
 
@@ -40,9 +41,30 @@ struct MatrixEntry {
     value: f64
 }
 
-struct PartialMatrix(Vec<MatrixEntry>);
+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;
@@ -64,13 +86,13 @@ impl PartialMatrix {
 
 // --- descent history ---
 
-struct DescentHistory {
+pub struct DescentHistory {
-    config: Vec<DMatrix<f64>>,
+    pub config: Vec<DMatrix<f64>>,
-    scaled_loss: Vec<f64>,
+    pub scaled_loss: Vec<f64>,
-    neg_grad: Vec<DMatrix<f64>>,
+    pub neg_grad: Vec<DMatrix<f64>>,
-    min_eigval: Vec<f64>,
+    pub min_eigval: Vec<f64>,
-    base_step: Vec<DMatrix<f64>>,
+    pub base_step: Vec<DMatrix<f64>>,
-    backoff_steps: Vec<i32>
+    pub backoff_steps: Vec<i32>
 }
 
 impl DescentHistory {
@@ -148,7 +170,7 @@ fn seek_better_config(
 
 // seek a matrix `config` for which `config' * Q * config` matches the partial
 // matrix `gram`. use gradient descent starting from `guess`
-fn realize_gram(
+pub fn realize_gram(
     gram: &PartialMatrix,
     guess: DMatrix<f64>,
     frozen: &[(usize, usize)],