From b117c00992209bad400e0b764ef09660c9ae1574 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 27 Mar 2025 13:57:46 -0700 Subject: [PATCH] Specify the values of the frozen entries Before, a `ConstraintProblem` only specified the indices of the frozen entries. During realization, the frozen entries kept whatever values they had in the initial guess. This commit adds the values of the frozen entries to the `frozen` field of `ConstraintProblem`. The frozen entries of the guess are set to the desired values at the beginning of realization. This commit also improves the `PartialMatrix` structure, which is used to specify the indices and values of the frozen entries. --- app-proto/examples/point-on-sphere.rs | 2 +- app-proto/src/engine.rs | 135 +++++++++++++++++--------- 2 files changed, 91 insertions(+), 46 deletions(-) diff --git a/app-proto/examples/point-on-sphere.rs b/app-proto/examples/point-on-sphere.rs index 2820793..880d7b0 100644 --- a/app-proto/examples/point-on-sphere.rs +++ b/app-proto/examples/point-on-sphere.rs @@ -10,7 +10,7 @@ fn main() { problem.gram.push_sym(j, k, if (j, k) == (1, 1) { 1.0 } else { 0.0 }); } } - problem.frozen.push((3, 0)); + problem.frozen.push(3, 0, problem.guess[(3, 0)]); println!(); let (config, _, success, history) = realize_gram( &problem, 1.0e-12, 0.5, 0.9, 1.1, 200, 110 diff --git a/app-proto/src/engine.rs b/app-proto/src/engine.rs index deb88bd..44f44e0 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -37,7 +37,7 @@ pub fn sphere_with_offset(dir_x: f64, dir_y: f64, dir_z: f64, off: f64, curv: f6 // --- partial matrices --- -struct MatrixEntry { +pub struct MatrixEntry { index: (usize, usize), value: f64 } @@ -49,42 +49,72 @@ impl PartialMatrix { PartialMatrix(Vec::::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) -> DMatrix { + let mut result = a.clone(); + for &MatrixEntry { index, value } in self { + result[index] = value; + } + result + } + fn proj(&self, a: &DMatrix) -> DMatrix { let mut result = DMatrix::::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) -> DMatrix { let mut result = DMatrix::::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; + + 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)] @@ -199,8 +229,8 @@ impl DescentHistory { pub struct ConstraintProblem { pub gram: PartialMatrix, + pub frozen: PartialMatrix, pub guess: DMatrix, - pub frozen: Vec<(usize, usize)> } impl ConstraintProblem { @@ -208,8 +238,8 @@ impl ConstraintProblem { const ELEMENT_DIM: usize = 5; ConstraintProblem { gram: PartialMatrix::new(), - guess: DMatrix::::zeros(ELEMENT_DIM, element_count), - frozen: Vec::new() + frozen: PartialMatrix::new(), + guess: DMatrix::::zeros(ELEMENT_DIM, element_count) } } @@ -217,8 +247,8 @@ impl ConstraintProblem { pub fn from_guess(guess_columns: &[DVector]) -> ConstraintProblem { ConstraintProblem { gram: PartialMatrix::new(), - guess: DMatrix::from_columns(guess_columns), - frozen: Vec::new() + frozen: PartialMatrix::new(), + guess: DMatrix::from_columns(guess_columns) } } } @@ -314,8 +344,10 @@ 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( problem: &ConstraintProblem, scaled_tol: f64, @@ -344,11 +376,11 @@ pub fn realize_gram( // convert the frozen indices to stacked format let frozen_stacked: Vec = 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.clone()); + // 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 @@ -501,7 +533,7 @@ pub mod examples { // the frozen entries fix the radii of the circumscribing sphere, the // "sun" and "moon" spheres, and one of the chain spheres for k in 0..4 { - problem.frozen.push((3, k)) + problem.frozen.push(3, k, problem.guess[(3, k)]); } realize_gram(&problem, scaled_tol, 0.5, 0.9, 1.1, 200, 110) @@ -545,7 +577,7 @@ pub mod examples { } for k in 0..N_POINTS { - problem.frozen.push((3, k)) + problem.frozen.push(3, k, problem.guess[(3, k)]) } realize_gram(&problem, scaled_tol, 0.5, 0.9, 1.1, 200, 110) @@ -559,6 +591,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::::from_row_slice(2, 3, &[ + 1.0, 2.0, 3.0, + 4.0, 5.0, 6.0 + ]); + let expected_result = DMatrix::::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![ @@ -580,18 +631,12 @@ mod tests { #[test] fn zero_loss_test() { - let gram = PartialMatrix({ - let mut entries = Vec::::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(&[ @@ -604,33 +649,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 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) ]); 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 }); } } - for k in 0..2 { - problem.frozen.push((3, k)); - } + problem.frozen.push(3, 0, problem.guess[(3, 0)]); + problem.frozen.push(3, 1, 0.5); let (config, _, success, history) = realize_gram( &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 &problem.frozen { + for &MatrixEntry { index, .. } in &problem.frozen { assert_eq!(base_step[index], 0.0); } } - for index in problem.frozen { - assert_eq!(config[index], problem.guess[index]); + for MatrixEntry { index, value } in problem.frozen { + assert_eq!(config[index], value); } } @@ -663,7 +708,7 @@ mod tests { } } for n in 0..ELEMENT_DIM { - problem.frozen.push((n, 0)); + problem.frozen.push(n, 0, problem.guess[(n, 0)]); } let (config, tangent, success, history) = realize_gram( &problem, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110