app-proto/Cargo.toml

@@ -10,7 +10,6 @@ default = ["console_error_panic_hook"]
 [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"

app-proto/run-examples

@@ -1,8 +0,0 @@
-# based on "Enabling print statements in Cargo tests", by Jon Almeida
-#
-#   https://jonalmeida.com/posts/2015/01/23/print-cargo/
-#
-
-cargo test -- --nocapture engine::tests::irisawa_hexlet_test
-cargo test -- --nocapture engine::tests::three_spheres_example
-cargo test -- --nocapture engine::tests::point_on_sphere_example

app-proto/src/add_remove.rs

@@ -12,8 +12,7 @@ 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(
@@ -22,8 +21,7 @@ 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(
@@ -32,8 +30,7 @@ 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(
@@ -42,8 +39,7 @@ 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(
@@ -52,8 +48,7 @@ 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(
@@ -62,8 +57,17 @@ 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)
         }
     );
 }
@@ -77,8 +81,7 @@ 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(
@@ -87,8 +90,7 @@ 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(
@@ -97,8 +99,7 @@ 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(
@@ -107,8 +108,7 @@ 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,8 +117,7 @@ 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(
@@ -127,8 +126,7 @@ 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(
@@ -137,8 +135,7 @@ 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(
@@ -147,8 +144,7 @@ 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
         }
     );
 }
@@ -219,7 +215,6 @@ 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,11 +1,8 @@
-use nalgebra::{DMatrix, DVector};
+use nalgebra::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 {
@@ -13,10 +10,7 @@ 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)]
@@ -46,8 +40,6 @@ 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
@@ -85,8 +77,7 @@ 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
             }
         );
     }
@@ -99,83 +90,4 @@ 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,12 +1,4 @@
-use lazy_static::lazy_static;
+use nalgebra::DVector;
-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> {
@@ -33,470 +25,3 @@ pub fn sphere_with_offset(dir_x: f64, dir_y: f64, dir_z: f64, off: f64, curv: f6
         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! {
-    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 ---
-
-#[cfg(test)]
-mod tests {
-    use std::{array, f64::consts::PI};
-    
-    use super::*;
-    
-    #[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);
-    }
-    
-    // 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/
-    //
-    #[test]
-    fn irisawa_hexlet_test() {
-        let gram = PartialMatrix({
-            let mut entries = Vec::<MatrixEntry>::new();
-            for s in 0..9 {
-                // each sphere is represented by a spacelike vector
-                entries.push(MatrixEntry { index: (s, s), value: 1.0 });
-                
-                // the circumscribing sphere is tangent to all of the other
-                // spheres, with matching orientation
-                if s > 0 {
-                    entries.push(MatrixEntry { index: (0, s), value: 1.0 });
-                    entries.push(MatrixEntry { index: (s, 0), value: 1.0 });
-                }
-                
-                if s > 2 {
-                    // each chain sphere is tangent to the "sun" and "moon"
-                    // spheres, with opposing orientation
-                    for n in 1..3 {
-                        entries.push(MatrixEntry { index: (s, n), value: -1.0 });
-                        entries.push(MatrixEntry { index: (n, s), value: -1.0 });
-                    }
-                    
-                    // each chain sphere is tangent to the next chain sphere,
-                    // with opposing orientation
-                    let s_next = 3 + (s-2) % 6;
-                    entries.push(MatrixEntry { index: (s, s_next), value: -1.0 });
-                    entries.push(MatrixEntry { index: (s_next, s), value: -1.0 });
-                }
-            }
-            entries
-        });
-        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()
-        );
-        let frozen: [(usize, usize); 4] = array::from_fn(|k| (3, k));
-        const SCALED_TOL: f64 = 1.0e-12;
-        let (config, success, history) = realize_gram(
-            &gram, guess, &frozen,
-            SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
-        );
-        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);
-        }
-        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);
-        }
-    }
-    
-    // --- process inspection examples ---
-    
-    // these tests are meant for human inspection, not automated use. run them
-    // one at a time in `--nocapture` mode and read through the results and
-    // optimization histories that they print out. the `run-examples` script
-    // will run all of them
-    
-    #[test]
-    fn three_spheres_example() {
-        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 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);
-        }
-    }
-    
-    #[test]
-    fn point_on_sphere_example() {
-        let gram = PartialMatrix({
-            let mut entries = Vec::<MatrixEntry>::new();
-            for j in 0..2 {
-                for k in 0..2 {
-                    entries.push(MatrixEntry {
-                        index: (j, k),
-                        value: if (j, k) == (1, 1) { 1.0 } else { 0.0 }
-                    });
-                }
-            }
-            entries
-        });
-        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);
-        }
-    }
-}

engine-proto/gram-test/Engine.jl

@@ -8,8 +8,7 @@ using Optim
 
 export
   rand_on_shell, Q, DescentHistory,
-  realize_gram_gradient, realize_gram_newton, realize_gram_optim,
+  realize_gram_gradient, realize_gram_newton, realize_gram_optim, realize_gram
-  realize_gram_alt_proj, realize_gram
 
 # === guessing ===
 
@@ -144,7 +143,7 @@ function realize_gram_gradient(
       break
     end
     
-    # find the negative gradient of the loss function
+    # find negative gradient of loss function
     neg_grad = 4*Q*L*Δ_proj
     slope = norm(neg_grad)
     dir = neg_grad / slope
@@ -233,7 +232,7 @@ function realize_gram_newton(
       break
     end
     
-    # find the negative gradient of the loss function
+    # find the negative gradient of loss function
     neg_grad = 4*Q*L*Δ_proj
     
     # find the negative Hessian of the loss function
@@ -314,129 +313,6 @@ 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(
@@ -445,6 +321,7 @@ 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,
@@ -475,19 +352,20 @@ function realize_gram(
     unfrozen_stacked = reshape(is_unfrozen, total_dim)
   end
   
-  # initialize search state
+  # initialize variables
+  grad_rate = init_rate
   L = copy(guess)
-  Δ_proj = proj_diff(gram, L'*Q*L)
-  loss = dot(Δ_proj, Δ_proj)
   
   # use Newton's method with backtracking and gradient descent backup
+  Δ_proj = proj_diff(gram, L'*Q*L)
+  loss = dot(Δ_proj, Δ_proj)
   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
+    # find the negative gradient of loss function
     neg_grad = 4*Q*L*Δ_proj
     
     # find the negative Hessian of the loss function
@@ -542,7 +420,6 @@ 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
@@ -551,7 +428,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 * base_target_improvement
+      if improvement >= min_efficiency * rate * dot(neg_grad, base_step)
         history.backoff_steps[end] = backoff_steps
         step_success = true
         break

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

@@ -75,12 +75,3 @@ if success
     println("  ", 1 / L[4,k], " sun")
   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

@@ -65,12 +65,3 @@ else
 end
 println("Steps: ", size(history.scaled_loss, 1))
 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

@@ -94,12 +94,3 @@ if success
   radius_ratio = dot(infty, Engine.Q * L[:,5]) / dot(infty, Engine.Q * L[:,6])
   println("\nCircumradius / inradius: ", radius_ratio)
 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")