feat: Curvature regulators (#80)

Prior to this commit, there's only one kind of regulator: the one that regulates the inversive distance between two spheres (or, more generally, the Lorentz product between two element representation vectors). Adds a new kind of regulator, which regulates the curvature of a sphere (issue #55). In the process, introduces a general framework based on new traits for organizing and sharing code between different kinds of regulators.

Co-authored-by: Aaron Fenyes <aaron.fenyes@fareycircles.ooo>
Reviewed-on: #80
Co-authored-by: Vectornaut <vectornaut@nobody@nowhere.net>
Co-committed-by: Vectornaut <vectornaut@nobody@nowhere.net>

app-proto/src/engine.rs

@@ -35,9 +35,43 @@ pub fn sphere_with_offset(dir_x: f64, dir_y: f64, dir_z: f64, off: f64, curv: f6
     ])
 }
 
+// given a sphere's representation vector, change the sphere's half-curvature to
+// `half-curv` and then restore normalization by contracting the representation
+// vector toward the curvature axis
+pub fn change_half_curvature(rep: &mut DVector<f64>, half_curv: f64) {
+    // set the sphere's half-curvature to the desired value
+    rep[3] = half_curv;
+    
+    // restore normalization by contracting toward the curvature axis
+    const SIZE_THRESHOLD: f64 = 1e-9;
+    let half_q_lt = -2.0 * half_curv * rep[4];
+    let half_q_lt_sq = half_q_lt * half_q_lt;
+    let mut spatial = rep.fixed_rows_mut::<3>(0);
+    let q_sp = spatial.norm_squared();
+    if q_sp < SIZE_THRESHOLD && half_q_lt_sq < SIZE_THRESHOLD {
+        spatial.copy_from_slice(
+            &[0.0, 0.0, (1.0 - 2.0 * half_q_lt).sqrt()]
+        );
+    } else {
+        let scaling = half_q_lt + (q_sp + half_q_lt_sq).sqrt();
+        spatial.scale_mut(1.0 / scaling);
+        rep[4] /= scaling;
+    }
+    
+    /* DEBUG */
+    // verify normalization
+    let rep_for_debug = rep.clone();
+    console::log_1(&JsValue::from(
+        format!(
+            "Sphere self-product after curvature change: {}",
+            rep_for_debug.dot(&(&*Q * &rep_for_debug))
+        )
+    ));
+}
+
 // --- partial matrices ---
 
-struct MatrixEntry {
+pub struct MatrixEntry {
     index: (usize, usize),
     value: f64
 }
@@ -49,42 +83,72 @@ impl PartialMatrix {
         PartialMatrix(Vec::<MatrixEntry>::new())
     }
     
-    pub fn push_sym(&mut self, row: usize, col: usize, value: f64) {
+    pub fn push(&mut self, row: usize, col: usize, value: f64) {
         let PartialMatrix(entries) = self;
         entries.push(MatrixEntry { index: (row, col), value: value });
+    }
+    
+    pub fn push_sym(&mut self, row: usize, col: usize, value: f64) {
+        self.push(row, col, value);
         if row != col {
-            entries.push(MatrixEntry { index: (col, row), value: value });
+            self.push(col, row, value);
         }
     }
     
     /* DEBUG */
     pub fn log_to_console(&self) {
-        let PartialMatrix(entries) = self;
-        for ent in entries {
-            let ent_str = format!("  {} {} {}", ent.index.0, ent.index.1, ent.value);
-            console::log_1(&JsValue::from(ent_str.as_str()));
+        for &MatrixEntry { index: (row, col), value } in self {
+            console::log_1(&JsValue::from(
+                format!("  {} {} {}", row, col, value)
+            ));
         }
     }
     
+    fn freeze(&self, a: &DMatrix<f64>) -> DMatrix<f64> {
+        let mut result = a.clone();
+        for &MatrixEntry { index, value } in self {
+            result[index] = value;
+        }
+        result
+    }
+    
     fn proj(&self, a: &DMatrix<f64>) -> DMatrix<f64> {
         let mut result = DMatrix::<f64>::zeros(a.nrows(), a.ncols());
-        let PartialMatrix(entries) = self;
-        for ent in entries {
-            result[ent.index] = a[ent.index];
+        for &MatrixEntry { index, .. } in self {
+            result[index] = a[index];
         }
         result
     }
     
     fn sub_proj(&self, rhs: &DMatrix<f64>) -> DMatrix<f64> {
         let mut result = DMatrix::<f64>::zeros(rhs.nrows(), rhs.ncols());
-        let PartialMatrix(entries) = self;
-        for ent in entries {
-            result[ent.index] = ent.value - rhs[ent.index];
+        for &MatrixEntry { index, value } in self {
+            result[index] = value - rhs[index];
         }
         result
     }
 }
 
+impl IntoIterator for PartialMatrix {
+    type Item = MatrixEntry;
+    type IntoIter = std::vec::IntoIter<Self::Item>;
+    
+    fn into_iter(self) -> Self::IntoIter {
+        let PartialMatrix(entries) = self;
+        entries.into_iter()
+    }
+}
+
+impl<'a> IntoIterator for &'a PartialMatrix {
+    type Item = &'a MatrixEntry;
+    type IntoIter = std::slice::Iter<'a, MatrixEntry>;
+    
+    fn into_iter(self) -> Self::IntoIter {
+        let PartialMatrix(entries) = self;
+        entries.into_iter()
+    }
+}
+
 // --- configuration subspaces ---
 
 #[derive(Clone)]
@@ -195,6 +259,34 @@ impl DescentHistory {
     }
 }
 
+// --- constraint problems ---
+
+pub struct ConstraintProblem {
+    pub gram: PartialMatrix,
+    pub frozen: PartialMatrix,
+    pub guess: DMatrix<f64>,
+}
+
+impl ConstraintProblem {
+    pub fn new(element_count: usize) -> ConstraintProblem {
+        const ELEMENT_DIM: usize = 5;
+        ConstraintProblem {
+            gram: PartialMatrix::new(),
+            frozen: PartialMatrix::new(),
+            guess: DMatrix::<f64>::zeros(ELEMENT_DIM, element_count)
+        }
+    }
+    
+    #[cfg(feature = "dev")]
+    pub fn from_guess(guess_columns: &[DVector<f64>]) -> ConstraintProblem {
+        ConstraintProblem {
+            gram: PartialMatrix::new(),
+            frozen: PartialMatrix::new(),
+            guess: DMatrix::from_columns(guess_columns)
+        }
+    }
+}
+
 // --- gram matrix realization ---
 
 // the Lorentz form
@@ -286,12 +378,12 @@ fn seek_better_config(
     None
 }
 
-// seek a matrix `config` for which `config' * Q * config` matches the partial
-// matrix `gram`. use gradient descent starting from `guess`
+// seek a matrix `config` that matches the partial matrix `problem.frozen` and
+// has `config' * Q * config` matching the partial matrix `problem.gram`. start
+// at `problem.guess`, set the frozen entries to their desired values, and then
+// use a regularized Newton's method to seek the desired Gram matrix
 pub fn realize_gram(
-    gram: &PartialMatrix,
-    guess: DMatrix<f64>,
-    frozen: &[(usize, usize)],
+    problem: &ConstraintProblem,
     scaled_tol: f64,
     min_efficiency: f64,
     backoff: f64,
@@ -299,6 +391,11 @@ pub fn realize_gram(
     max_descent_steps: i32,
     max_backoff_steps: i32
 ) -> (DMatrix<f64>, ConfigSubspace, bool, DescentHistory) {
+    // destructure the problem data
+    let ConstraintProblem {
+        gram, guess, frozen
+    } = problem;
+    
     // start the descent history
     let mut history = DescentHistory::new();
     
@@ -313,11 +410,11 @@ pub fn realize_gram(
     
     // convert the frozen indices to stacked format
     let frozen_stacked: Vec<usize> = frozen.into_iter().map(
-        |index| index.1*element_dim + index.0
+        |MatrixEntry { index: (row, col), .. }| col*element_dim + row
     ).collect();
     
-    // use Newton's method with backtracking and gradient descent backup
-    let mut state = SearchState::from_config(gram, guess);
+    // use a regularized Newton's method with backtracking
+    let mut state = SearchState::from_config(gram, frozen.freeze(guess));
     let mut hess = DMatrix::zeros(element_dim, assembly_dim);
     for _ in 0..max_descent_steps {
         // find the negative gradient of the loss function
@@ -415,7 +512,7 @@ pub fn realize_gram(
 
 #[cfg(feature = "dev")]
 pub mod examples {
-    use std::{array, f64::consts::PI};
+    use std::f64::consts::PI;
     
     use super::*;
     
@@ -428,35 +525,7 @@ pub mod examples {
     //   https://www.nippon.com/en/japan-topics/c12801/
     //
     pub fn realize_irisawa_hexlet(scaled_tol: f64) -> (DMatrix<f64>, ConfigSubspace, bool, DescentHistory) {
-        let gram = {
-            let mut gram_to_be = PartialMatrix::new();
-            for s in 0..9 {
-                // each sphere is represented by a spacelike vector
-                gram_to_be.push_sym(s, s, 1.0);
-                
-                // the circumscribing sphere is tangent to all of the other
-                // spheres, with matching orientation
-                if s > 0 {
-                    gram_to_be.push_sym(0, s, 1.0);
-                }
-                
-                if s > 2 {
-                    // each chain sphere is tangent to the "sun" and "moon"
-                    // spheres, with opposing orientation
-                    for n in 1..3 {
-                        gram_to_be.push_sym(s, n, -1.0);
-                    }
-                    
-                    // each chain sphere is tangent to the next chain sphere,
-                    // with opposing orientation
-                    let s_next = 3 + (s-2) % 6;
-                    gram_to_be.push_sym(s, s_next, -1.0);
-                }
-            }
-            gram_to_be
-        };
-        
-        let guess = DMatrix::from_columns(
+        let mut problem = ConstraintProblem::from_guess(
             [
                 sphere(0.0, 0.0, 0.0, 15.0),
                 sphere(0.0, 0.0, -9.0, 5.0),
@@ -471,42 +540,45 @@ pub mod examples {
             ).collect::<Vec<_>>().as_slice()
         );
         
+        for s in 0..9 {
+            // each sphere is represented by a spacelike vector
+            problem.gram.push_sym(s, s, 1.0);
+            
+            // the circumscribing sphere is tangent to all of the other
+            // spheres, with matching orientation
+            if s > 0 {
+                problem.gram.push_sym(0, s, 1.0);
+            }
+            
+            if s > 2 {
+                // each chain sphere is tangent to the "sun" and "moon"
+                // spheres, with opposing orientation
+                for n in 1..3 {
+                    problem.gram.push_sym(s, n, -1.0);
+                }
+                
+                // each chain sphere is tangent to the next chain sphere,
+                // with opposing orientation
+                let s_next = 3 + (s-2) % 6;
+                problem.gram.push_sym(s, s_next, -1.0);
+            }
+        }
+        
         // the frozen entries fix the radii of the circumscribing sphere, the
         // "sun" and "moon" spheres, and one of the chain spheres
-        let frozen: [(usize, usize); 4] = array::from_fn(|k| (3, k));
+        for k in 0..4 {
+            problem.frozen.push(3, k, problem.guess[(3, k)]);
+        }
         
-        realize_gram(
-            &gram, guess, &frozen,
-            scaled_tol, 0.5, 0.9, 1.1, 200, 110
-        )
+        realize_gram(&problem, scaled_tol, 0.5, 0.9, 1.1, 200, 110)
     }
     
     // set up a kaleidocycle, made of points with fixed distances between them,
     // and find its tangent space
     pub fn realize_kaleidocycle(scaled_tol: f64) -> (DMatrix<f64>, ConfigSubspace, bool, DescentHistory) {
-        const N_POINTS: usize = 12;
-        let gram = {
-            let mut gram_to_be = PartialMatrix::new();
-            for block in (0..N_POINTS).step_by(2) {
-                let block_next = (block + 2) % N_POINTS;
-                for j in 0..2 {
-                    // diagonal and hinge edges
-                    for k in j..2 {
-                        gram_to_be.push_sym(block + j, block + k, if j == k { 0.0 } else { -0.5 });
-                    }
-                    
-                    // non-hinge edges
-                    for k in 0..2 {
-                        gram_to_be.push_sym(block + j, block_next + k, -0.625);
-                    }
-                }
-            }
-            gram_to_be
-        };
-        
-        let guess = {
-            const N_HINGES: usize = 6;
-            let guess_elts = (0..N_HINGES).step_by(2).flat_map(
+        const N_HINGES: usize = 6;
+        let mut problem = ConstraintProblem::from_guess(
+            (0..N_HINGES).step_by(2).flat_map(
                 |n| {
                     let ang_hor = (n as f64) * PI/3.0;
                     let ang_vert = ((n + 1) as f64) * PI/3.0;
@@ -519,16 +591,30 @@ pub mod examples {
                         point(x_vert, y_vert, 0.5)
                     ]
                 }
-            ).collect::<Vec<_>>();
-            DMatrix::from_columns(&guess_elts)
-        };
+            ).collect::<Vec<_>>().as_slice()
+        );
         
-        let frozen: [_; N_POINTS] = array::from_fn(|k| (3, k));
+        const N_POINTS: usize = 2 * N_HINGES;
+        for block in (0..N_POINTS).step_by(2) {
+            let block_next = (block + 2) % N_POINTS;
+            for j in 0..2 {
+                // diagonal and hinge edges
+                for k in j..2 {
+                    problem.gram.push_sym(block + j, block + k, if j == k { 0.0 } else { -0.5 });
+                }
+                
+                // non-hinge edges
+                for k in 0..2 {
+                    problem.gram.push_sym(block + j, block_next + k, -0.625);
+                }
+            }
+        }
         
-        realize_gram(
-            &gram, guess, &frozen,
-            scaled_tol, 0.5, 0.9, 1.1, 200, 110
-        )
+        for k in 0..N_POINTS {
+            problem.frozen.push(3, k, problem.guess[(3, k)])
+        }
+        
+        realize_gram(&problem, scaled_tol, 0.5, 0.9, 1.1, 200, 110)
     }
 }
 
@@ -539,6 +625,25 @@ mod tests {
     
     use super::{*, examples::*};
     
+    #[test]
+    fn freeze_test() {
+        let frozen = PartialMatrix(vec![
+            MatrixEntry { index: (0, 0), value: 14.0 },
+            MatrixEntry { index: (0, 2), value: 28.0 },
+            MatrixEntry { index: (1, 1), value: 42.0 },
+            MatrixEntry { index: (1, 2), value: 49.0 }
+        ]);
+        let config = DMatrix::<f64>::from_row_slice(2, 3, &[
+            1.0, 2.0, 3.0,
+            4.0, 5.0, 6.0
+        ]);
+        let expected_result = DMatrix::<f64>::from_row_slice(2, 3, &[
+            14.0, 2.0, 28.0,
+            4.0, 42.0, 49.0
+        ]);
+        assert_eq!(frozen.freeze(&config), expected_result);
+    }
+    
     #[test]
     fn sub_proj_test() {
         let target = PartialMatrix(vec![
@@ -560,18 +665,12 @@ mod tests {
     
     #[test]
     fn zero_loss_test() {
-        let gram = PartialMatrix({
-            let mut entries = Vec::<MatrixEntry>::new();
-            for j in 0..3 {
-                for k in 0..3 {
-                    entries.push(MatrixEntry {
-                        index: (j, k),
-                        value: if j == k { 1.0 } else { -1.0 }
-                    });
-                }
+        let mut gram = PartialMatrix::new();
+        for j in 0..3 {
+            for k in 0..3 {
+                gram.push(j, k, if j == k { 1.0 } else { -1.0 });
             }
-            entries
-        });
+        }
         let config = {
             let a = 0.75_f64.sqrt();
             DMatrix::from_columns(&[
@@ -584,37 +683,33 @@ mod tests {
         assert!(state.loss.abs() < f64::EPSILON);
     }
     
+    /* TO DO */
     // at the frozen indices, the optimization steps should have exact zeros,
-    // and the realized configuration should match the initial guess
+    // and the realized configuration should have the desired values
     #[test]
     fn frozen_entry_test() {
-        let gram = {
-            let mut gram_to_be = PartialMatrix::new();
-            for j in 0..2 {
-                for k in j..2 {
-                    gram_to_be.push_sym(j, k, if (j, k) == (1, 1) { 1.0 } else { 0.0 });
-                }
-            }
-            gram_to_be
-        };
-        let guess = DMatrix::from_columns(&[
+        let mut problem = ConstraintProblem::from_guess(&[
             point(0.0, 0.0, 2.0),
-            sphere(0.0, 0.0, 0.0, 1.0)
+            sphere(0.0, 0.0, 0.0, 0.95)
         ]);
-        let frozen = [(3, 0), (3, 1)];
-        println!();
+        for j in 0..2 {
+            for k in j..2 {
+                problem.gram.push_sym(j, k, if (j, k) == (1, 1) { 1.0 } else { 0.0 });
+            }
+        }
+        problem.frozen.push(3, 0, problem.guess[(3, 0)]);
+        problem.frozen.push(3, 1, 0.5);
         let (config, _, success, history) = realize_gram(
-            &gram, guess.clone(), &frozen,
-            1.0e-12, 0.5, 0.9, 1.1, 200, 110
+            &problem, 1.0e-12, 0.5, 0.9, 1.1, 200, 110
         );
         assert_eq!(success, true);
         for base_step in history.base_step.into_iter() {
-            for index in frozen {
+            for &MatrixEntry { index, .. } in &problem.frozen {
                 assert_eq!(base_step[index], 0.0);
             }
         }
-        for index in frozen {
-            assert_eq!(config[index], guess[index]);
+        for MatrixEntry { index, value } in problem.frozen {
+            assert_eq!(config[index], value);
         }
     }
     
@@ -635,34 +730,32 @@ mod tests {
     #[test]
     fn tangent_test_three_spheres() {
         const SCALED_TOL: f64 = 1.0e-12;
-        let gram = {
-            let mut gram_to_be = PartialMatrix::new();
-            for j in 0..3 {
-                for k in j..3 {
-                    gram_to_be.push_sym(j, k, if j == k { 1.0 } else { -1.0 });
-                }
-            }
-            gram_to_be
-        };
-        let guess = DMatrix::from_columns(&[
+        const ELEMENT_DIM: usize = 5;
+        let mut problem = ConstraintProblem::from_guess(&[
             sphere(0.0, 0.0, 0.0, -2.0),
             sphere(0.0, 0.0, 1.0, 1.0),
             sphere(0.0, 0.0, -1.0, 1.0)
         ]);
-        let frozen: [_; 5] = std::array::from_fn(|k| (k, 0));
+        for j in 0..3 {
+            for k in j..3 {
+                problem.gram.push_sym(j, k, if j == k { 1.0 } else { -1.0 });
+            }
+        }
+        for n in 0..ELEMENT_DIM {
+            problem.frozen.push(n, 0, problem.guess[(n, 0)]);
+        }
         let (config, tangent, success, history) = realize_gram(
-            &gram, guess.clone(), &frozen,
-            SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
+            &problem, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
         );
-        assert_eq!(config, guess);
+        assert_eq!(config, problem.guess);
         assert_eq!(success, true);
         assert_eq!(history.scaled_loss.len(), 1);
         
         // list some motions that should form a basis for the tangent space of
         // the solution variety
         const UNIFORM_DIM: usize = 4;
-        let element_dim = guess.nrows();
-        let assembly_dim = guess.ncols();
+        let element_dim = problem.guess.nrows();
+        let assembly_dim = problem.guess.ncols();
         let tangent_motions_unif = vec![
             basis_matrix((0, 1), UNIFORM_DIM, assembly_dim),
             basis_matrix((1, 1), UNIFORM_DIM, assembly_dim),
@@ -805,22 +898,17 @@ mod tests {
     fn proj_equivar_test() {
         // find a pair of spheres that meet at 120°
         const SCALED_TOL: f64 = 1.0e-12;
-        let gram = {
-            let mut gram_to_be = PartialMatrix::new();
-            gram_to_be.push_sym(0, 0, 1.0);
-            gram_to_be.push_sym(1, 1, 1.0);
-            gram_to_be.push_sym(0, 1, 0.5);
-            gram_to_be
-        };
-        let guess_orig = DMatrix::from_columns(&[
+        let mut problem_orig = ConstraintProblem::from_guess(&[
             sphere(0.0, 0.0, 0.5, 1.0),
             sphere(0.0, 0.0, -0.5, 1.0)
         ]);
+        problem_orig.gram.push_sym(0, 0, 1.0);
+        problem_orig.gram.push_sym(1, 1, 1.0);
+        problem_orig.gram.push_sym(0, 1, 0.5);
         let (config_orig, tangent_orig, success_orig, history_orig) = realize_gram(
-            &gram, guess_orig.clone(), &[],
-            SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
+            &problem_orig, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
         );
-        assert_eq!(config_orig, guess_orig);
+        assert_eq!(config_orig, problem_orig.guess);
         assert_eq!(success_orig, true);
         assert_eq!(history_orig.scaled_loss.len(), 1);
         
@@ -833,11 +921,15 @@ mod tests {
                 sphere(-a, 0.0, 7.0 - a, 1.0)
             ])
         };
+        let problem_tfm = ConstraintProblem {
+            gram: problem_orig.gram,
+            guess: guess_tfm,
+            frozen: problem_orig.frozen
+        };
         let (config_tfm, tangent_tfm, success_tfm, history_tfm) = realize_gram(
-            &gram, guess_tfm.clone(), &[],
-            SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
+            &problem_tfm, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
         );
-        assert_eq!(config_tfm, guess_tfm);
+        assert_eq!(config_tfm, problem_tfm.guess);
         assert_eq!(success_tfm, true);
         assert_eq!(history_tfm.scaled_loss.len(), 1);
         
@@ -869,7 +961,7 @@ mod tests {
         // the comparison tolerance because the transformation seems to
         // introduce some numerical error
         const SCALED_TOL_TFM: f64 = 1.0e-9;
-        let tol_sq = ((guess_orig.nrows() * guess_orig.ncols()) as f64) * SCALED_TOL_TFM * SCALED_TOL_TFM;
+        let tol_sq = ((problem_orig.guess.nrows() * problem_orig.guess.ncols()) as f64) * SCALED_TOL_TFM * SCALED_TOL_TFM;
         assert!((motion_proj_tfm - motion_tfm_proj).norm_squared() < tol_sq);
     }
 }