Encapsulate the constraint problem data

This will make it easier for elements and regulators to write themselves
into the constraint problem.

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

@@ -1,26 +1,19 @@
-use nalgebra::DMatrix;
+use dyna3::engine::{Q, point, realize_gram, sphere, ConstraintProblem};
-
-use dyna3::engine::{Q, point, realize_gram, sphere, PartialMatrix};
 
 fn main() {
-    let gram = {
+    let mut problem = ConstraintProblem::from_guess(&[
-        let mut gram_to_be = PartialMatrix::new();
-        for j in 0..2 {
-            for k in j..2 {
-                gram_to_be.push_sym(j, k, if (j, k) == (1, 1) { 1.0 } else { 0.0 });
-            }
-        }
-        gram_to_be
-    };
-    let guess = DMatrix::from_columns(&[
         point(0.0, 0.0, 2.0),
         sphere(0.0, 0.0, 0.0, 1.0)
     ]);
-    let frozen = [(3, 0)];
+    for j in 0..2 {
+        for k in j..2 {
+            problem.gram.push_sym(j, k, if (j, k) == (1, 1) { 1.0 } else { 0.0 });
+        }
+    }
+    problem.frozen.push((3, 0));
     println!();
     let (config, _, success, history) = realize_gram(
-        &gram, guess, &frozen,
+        &problem, 1.0e-12, 0.5, 0.9, 1.1, 200, 110
-        1.0e-12, 0.5, 0.9, 1.1, 200, 110
     );
     print!("\nCompleted Gram matrix:{}", config.tr_mul(&*Q) * &config);
     print!("Configuration:{}", config);

app-proto/examples/three-spheres.rs

@@ -1,29 +1,22 @@
-use nalgebra::DMatrix;
+use dyna3::engine::{Q, realize_gram, sphere, ConstraintProblem};
-
-use dyna3::engine::{Q, realize_gram, sphere, PartialMatrix};
 
 fn main() {
-    let gram = {
+    let mut problem = ConstraintProblem::from_guess({
-        let mut gram_to_be = PartialMatrix::new();
-        for j in 0..3 {
-            for k in j..3 {
-                gram_to_be.push_sym(j, k, if j == k { 1.0 } else { -1.0 });
-            }
-        }
-        gram_to_be
-    };
-    let guess = {
         let a: f64 = 0.75_f64.sqrt();
-        DMatrix::from_columns(&[
+        &[
             sphere(1.0, 0.0, 0.0, 1.0),
             sphere(-0.5, a, 0.0, 1.0),
             sphere(-0.5, -a, 0.0, 1.0)
-        ])
+        ]
-    };
+    });
+    for j in 0..3 {
+        for k in j..3 {
+            problem.gram.push_sym(j, k, if j == k { 1.0 } else { -1.0 });
+        }
+    }
     println!();
     let (config, _, success, history) = realize_gram(
-        &gram, guess, &[],
+        &problem, 1.0e-12, 0.5, 0.9, 1.1, 200, 110
-        1.0e-12, 0.5, 0.9, 1.1, 200, 110
     );
     print!("\nCompleted Gram matrix:{}", config.tr_mul(&*Q) * &config);
     if success {

app-proto/src/assembly.rs

@@ -6,7 +6,13 @@ use sycamore::prelude::*;
 use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */
 
 use crate::{
-    engine::{Q, local_unif_to_std, realize_gram, ConfigSubspace, PartialMatrix},
+    engine::{
+        Q,
+        local_unif_to_std,
+        realize_gram,
+        ConfigSubspace,
+        ConstraintProblem
+    },
     specified::SpecifiedValue
 };
 
@@ -268,6 +274,11 @@ impl Assembly {
     // --- realization ---
     
     pub fn realize(&self) {
+        // create a blank constraint problem
+        let mut problem = ConstraintProblem::new(
+            self.elements.with_untracked(|elts| elts.len())
+        );
+        
         // index the elements
         self.elements.update_silent(|elts| {
             for (index, (_, elt)) in elts.into_iter().enumerate() {
@@ -276,9 +287,8 @@ impl Assembly {
         });
         
         // set up the Gram matrix and the initial configuration matrix
-        let (gram, guess) = self.elements.with_untracked(|elts| {
+        self.elements.with_untracked(|elts| {
             // set up the off-diagonal part of the Gram matrix
-            let mut gram_to_be = PartialMatrix::new();
             self.regulators.with_untracked(|regs| {
                 for (_, reg) in regs {
                     reg.set_point.with_untracked(|set_pt| {
@@ -286,7 +296,7 @@ impl Assembly {
                             let subjects = reg.subjects;
                             let row = elts[subjects.0].column_index.unwrap();
                             let col = elts[subjects.1].column_index.unwrap();
-                            gram_to_be.push_sym(row, col, val);
+                            problem.gram.push_sym(row, col, val);
                         }
                     });
                 }
@@ -294,36 +304,32 @@ impl Assembly {
             
             // set up the initial configuration matrix and the diagonal of the
             // Gram matrix
-            let mut guess_to_be = DMatrix::<f64>::zeros(5, elts.len());
             for (_, elt) in elts {
                 let index = elt.column_index.unwrap();
-                gram_to_be.push_sym(index, index, 1.0);
+                problem.gram.push_sym(index, index, 1.0);
-                guess_to_be.set_column(index, &elt.representation.get_clone_untracked());
+                problem.guess.set_column(index, &elt.representation.get_clone_untracked());
             }
-            
-            (gram_to_be, guess_to_be)
         });
         
         /* DEBUG */
         // log the Gram matrix
         console::log_1(&JsValue::from("Gram matrix:"));
-        gram.log_to_console();
+        problem.gram.log_to_console();
         
         /* DEBUG */
         // log the initial configuration matrix
         console::log_1(&JsValue::from("Old configuration:"));
-        for j in 0..guess.nrows() {
+        for j in 0..problem.guess.nrows() {
             let mut row_str = String::new();
-            for k in 0..guess.ncols() {
+            for k in 0..problem.guess.ncols() {
-                row_str.push_str(format!(" {:>8.3}", guess[(j, k)]).as_str());
+                row_str.push_str(format!(" {:>8.3}", problem.guess[(j, k)]).as_str());
             }
             console::log_1(&JsValue::from(row_str));
         }
         
         // look for a configuration with the given Gram matrix
         let (config, tangent, success, history) = realize_gram(
-            &gram, guess, &[],
+            &problem, 1.0e-12, 0.5, 0.9, 1.1, 200, 110
-            1.0e-12, 0.5, 0.9, 1.1, 200, 110
         );
         
         /* DEBUG */

app-proto/src/engine.rs

@@ -195,6 +195,34 @@ impl DescentHistory {
     }
 }
 
+// --- constraint problems ---
+
+pub struct ConstraintProblem {
+    pub gram: PartialMatrix,
+    pub guess: DMatrix<f64>,
+    pub frozen: Vec<(usize, usize)>
+}
+
+impl ConstraintProblem {
+    pub fn new(element_count: usize) -> ConstraintProblem {
+        const ELEMENT_DIM: usize = 5;
+        ConstraintProblem {
+            gram: PartialMatrix::new(),
+            guess: DMatrix::<f64>::zeros(ELEMENT_DIM, element_count),
+            frozen: Vec::new()
+        }
+    }
+    
+    #[cfg(feature = "dev")]
+    pub fn from_guess(guess_columns: &[DVector<f64>]) -> ConstraintProblem {
+        ConstraintProblem {
+            gram: PartialMatrix::new(),
+            guess: DMatrix::from_columns(guess_columns),
+            frozen: Vec::new()
+        }
+    }
+}
+
 // --- gram matrix realization ---
 
 // the Lorentz form
@@ -289,9 +317,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`
 pub fn realize_gram(
-    gram: &PartialMatrix,
+    problem: &ConstraintProblem,
-    guess: DMatrix<f64>,
-    frozen: &[(usize, usize)],
     scaled_tol: f64,
     min_efficiency: f64,
     backoff: f64,
@@ -299,6 +325,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();
     
@@ -317,7 +348,7 @@ pub fn realize_gram(
     ).collect();
     
     // use Newton's method with backtracking and gradient descent backup
-    let mut state = SearchState::from_config(gram, guess);
+    let mut state = SearchState::from_config(gram, guess.clone());
     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 +446,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 +459,7 @@ pub mod examples {
     //   https://www.nippon.com/en/japan-topics/c12801/
     //
     pub fn realize_irisawa_hexlet(scaled_tol: f64) -> (DMatrix<f64>, ConfigSubspace, bool, DescentHistory) {
-        let gram = {
+        let mut problem = ConstraintProblem::from_guess(
-            let mut gram_to_be = PartialMatrix::new();
-            for s in 0..9 {
-                // each sphere is represented by a spacelike vector
-                gram_to_be.push_sym(s, s, 1.0);
-                
-                // the circumscribing sphere is tangent to all of the other
-                // spheres, with matching orientation
-                if s > 0 {
-                    gram_to_be.push_sym(0, s, 1.0);
-                }
-                
-                if s > 2 {
-                    // each chain sphere is tangent to the "sun" and "moon"
-                    // spheres, with opposing orientation
-                    for n in 1..3 {
-                        gram_to_be.push_sym(s, n, -1.0);
-                    }
-                    
-                    // each chain sphere is tangent to the next chain sphere,
-                    // with opposing orientation
-                    let s_next = 3 + (s-2) % 6;
-                    gram_to_be.push_sym(s, s_next, -1.0);
-                }
-            }
-            gram_to_be
-        };
-        
-        let guess = DMatrix::from_columns(
             [
                 sphere(0.0, 0.0, 0.0, 15.0),
                 sphere(0.0, 0.0, -9.0, 5.0),
@@ -471,42 +474,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))
+        }
         
-        realize_gram(
+        realize_gram(&problem, scaled_tol, 0.5, 0.9, 1.1, 200, 110)
-            &gram, guess, &frozen,
-            scaled_tol, 0.5, 0.9, 1.1, 200, 110
-        )
     }
     
     // set up a kaleidocycle, made of points with fixed distances between them,
     // and find its tangent space
     pub fn realize_kaleidocycle(scaled_tol: f64) -> (DMatrix<f64>, ConfigSubspace, bool, DescentHistory) {
-        const N_POINTS: usize = 12;
+        const N_HINGES: usize = 6;
-        let gram = {
+        let mut problem = ConstraintProblem::from_guess(
-            let mut gram_to_be = PartialMatrix::new();
+            (0..N_HINGES).step_by(2).flat_map(
-            for block in (0..N_POINTS).step_by(2) {
-                let block_next = (block + 2) % N_POINTS;
-                for j in 0..2 {
-                    // diagonal and hinge edges
-                    for k in j..2 {
-                        gram_to_be.push_sym(block + j, block + k, if j == k { 0.0 } else { -0.5 });
-                    }
-                    
-                    // non-hinge edges
-                    for k in 0..2 {
-                        gram_to_be.push_sym(block + j, block_next + k, -0.625);
-                    }
-                }
-            }
-            gram_to_be
-        };
-        
-        let guess = {
-            const N_HINGES: usize = 6;
-            let guess_elts = (0..N_HINGES).step_by(2).flat_map(
                 |n| {
                     let ang_hor = (n as f64) * PI/3.0;
                     let ang_vert = ((n + 1) as f64) * PI/3.0;
@@ -519,16 +525,30 @@ pub mod examples {
                         point(x_vert, y_vert, 0.5)
                     ]
                 }
-            ).collect::<Vec<_>>();
+            ).collect::<Vec<_>>().as_slice()
-            DMatrix::from_columns(&guess_elts)
+        );
-        };
         
-        let frozen: [_; N_POINTS] = array::from_fn(|k| (3, k));
+        const N_POINTS: usize = 2 * N_HINGES;
+        for block in (0..N_POINTS).step_by(2) {
+            let block_next = (block + 2) % N_POINTS;
+            for j in 0..2 {
+                // diagonal and hinge edges
+                for k in j..2 {
+                    problem.gram.push_sym(block + j, block + k, if j == k { 0.0 } else { -0.5 });
+                }
+                
+                // non-hinge edges
+                for k in 0..2 {
+                    problem.gram.push_sym(block + j, block_next + k, -0.625);
+                }
+            }
+        }
         
-        realize_gram(
+        for k in 0..N_POINTS {
-            &gram, guess, &frozen,
+            problem.frozen.push((3, k))
-            scaled_tol, 0.5, 0.9, 1.1, 200, 110
+        }
-        )
+        
+        realize_gram(&problem, scaled_tol, 0.5, 0.9, 1.1, 200, 110)
     }
 }
 
@@ -588,33 +608,29 @@ mod tests {
     // and the realized configuration should match the initial guess
     #[test]
     fn frozen_entry_test() {
-        let gram = {
+        let mut problem = ConstraintProblem::from_guess(&[
-            let mut gram_to_be = PartialMatrix::new();
-            for j in 0..2 {
-                for k in j..2 {
-                    gram_to_be.push_sym(j, k, if (j, k) == (1, 1) { 1.0 } else { 0.0 });
-                }
-            }
-            gram_to_be
-        };
-        let guess = DMatrix::from_columns(&[
             point(0.0, 0.0, 2.0),
             sphere(0.0, 0.0, 0.0, 1.0)
         ]);
-        let frozen = [(3, 0), (3, 1)];
+        for j in 0..2 {
-        println!();
+            for k in j..2 {
+                problem.gram.push_sym(j, k, if (j, k) == (1, 1) { 1.0 } else { 0.0 });
+            }
+        }
+        for k in 0..2 {
+            problem.frozen.push((3, k));
+        }
         let (config, _, success, history) = realize_gram(
-            &gram, guess.clone(), &frozen,
+            &problem, 1.0e-12, 0.5, 0.9, 1.1, 200, 110
-            1.0e-12, 0.5, 0.9, 1.1, 200, 110
         );
         assert_eq!(success, true);
         for base_step in history.base_step.into_iter() {
-            for index in frozen {
+            for &index in &problem.frozen {
                 assert_eq!(base_step[index], 0.0);
             }
         }
-        for index in frozen {
+        for index in problem.frozen {
-            assert_eq!(config[index], guess[index]);
+            assert_eq!(config[index], problem.guess[index]);
         }
     }
     
@@ -635,34 +651,32 @@ mod tests {
     #[test]
     fn tangent_test_three_spheres() {
         const SCALED_TOL: f64 = 1.0e-12;
-        let gram = {
+        const ELEMENT_DIM: usize = 5;
-            let mut gram_to_be = PartialMatrix::new();
+        let mut problem = ConstraintProblem::from_guess(&[
-            for j in 0..3 {
-                for k in j..3 {
-                    gram_to_be.push_sym(j, k, if j == k { 1.0 } else { -1.0 });
-                }
-            }
-            gram_to_be
-        };
-        let guess = DMatrix::from_columns(&[
             sphere(0.0, 0.0, 0.0, -2.0),
             sphere(0.0, 0.0, 1.0, 1.0),
             sphere(0.0, 0.0, -1.0, 1.0)
         ]);
-        let frozen: [_; 5] = std::array::from_fn(|k| (k, 0));
+        for j in 0..3 {
+            for k in j..3 {
+                problem.gram.push_sym(j, k, if j == k { 1.0 } else { -1.0 });
+            }
+        }
+        for n in 0..ELEMENT_DIM {
+            problem.frozen.push((n, 0));
+        }
         let (config, tangent, success, history) = realize_gram(
-            &gram, guess.clone(), &frozen,
+            &problem, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
-            SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
         );
-        assert_eq!(config, guess);
+        assert_eq!(config, problem.guess);
         assert_eq!(success, true);
         assert_eq!(history.scaled_loss.len(), 1);
         
         // list some motions that should form a basis for the tangent space of
         // the solution variety
         const UNIFORM_DIM: usize = 4;
-        let element_dim = guess.nrows();
+        let element_dim = problem.guess.nrows();
-        let assembly_dim = guess.ncols();
+        let assembly_dim = problem.guess.ncols();
         let tangent_motions_unif = vec![
             basis_matrix((0, 1), UNIFORM_DIM, assembly_dim),
             basis_matrix((1, 1), UNIFORM_DIM, assembly_dim),
@@ -805,22 +819,17 @@ mod tests {
     fn proj_equivar_test() {
         // find a pair of spheres that meet at 120°
         const SCALED_TOL: f64 = 1.0e-12;
-        let gram = {
+        let mut problem_orig = ConstraintProblem::from_guess(&[
-            let mut gram_to_be = PartialMatrix::new();
-            gram_to_be.push_sym(0, 0, 1.0);
-            gram_to_be.push_sym(1, 1, 1.0);
-            gram_to_be.push_sym(0, 1, 0.5);
-            gram_to_be
-        };
-        let guess_orig = DMatrix::from_columns(&[
             sphere(0.0, 0.0, 0.5, 1.0),
             sphere(0.0, 0.0, -0.5, 1.0)
         ]);
+        problem_orig.gram.push_sym(0, 0, 1.0);
+        problem_orig.gram.push_sym(1, 1, 1.0);
+        problem_orig.gram.push_sym(0, 1, 0.5);
         let (config_orig, tangent_orig, success_orig, history_orig) = realize_gram(
-            &gram, guess_orig.clone(), &[],
+            &problem_orig, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
-            SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
         );
-        assert_eq!(config_orig, guess_orig);
+        assert_eq!(config_orig, problem_orig.guess);
         assert_eq!(success_orig, true);
         assert_eq!(history_orig.scaled_loss.len(), 1);
         
@@ -833,11 +842,15 @@ mod tests {
                 sphere(-a, 0.0, 7.0 - a, 1.0)
             ])
         };
+        let problem_tfm = ConstraintProblem {
+            gram: problem_orig.gram,
+            guess: guess_tfm,
+            frozen: problem_orig.frozen
+        };
         let (config_tfm, tangent_tfm, success_tfm, history_tfm) = realize_gram(
-            &gram, guess_tfm.clone(), &[],
+            &problem_tfm, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
-            SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
         );
-        assert_eq!(config_tfm, guess_tfm);
+        assert_eq!(config_tfm, problem_tfm.guess);
         assert_eq!(success_tfm, true);
         assert_eq!(history_tfm.scaled_loss.len(), 1);
         
@@ -869,7 +882,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);
     }
 }