Vectornaut/dyna3
1
0
Fork 0
forked from StudioInfinity/dyna3
Code Pull requests Activity

Specify the values of the frozen entries

Browse source 
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.
This commit is contained in:
Aaron Fenyes 2025-03-27 13:57:46 -07:00
parent ca9de34427
commit b117c00992
2 changed files with 91 additions and 46 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
app-proto/examples/point-on-sphere.rs
View file

@ -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

129
app-proto/src/engine.rs
View file

@ -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::<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)]
@ -199,8 +229,8 @@ impl DescentHistory {
pub struct ConstraintProblem {
pub gram: PartialMatrix,
pub frozen: PartialMatrix,
pub guess: DMatrix<f64>,
pub frozen: Vec<(usize, usize)>
}
impl ConstraintProblem {
@ -208,8 +238,8 @@ impl ConstraintProblem {
const ELEMENT_DIM: usize = 5;
ConstraintProblem {
gram: PartialMatrix::new(),
guess: DMatrix::<f64>::zeros(ELEMENT_DIM, element_count),
frozen: Vec::new()
frozen: PartialMatrix::new(),
guess: DMatrix::<f64>::zeros(ELEMENT_DIM, element_count)
}
}
@ -217,8 +247,8 @@ impl ConstraintProblem {
pub fn from_guess(guess_columns: &[DVector<f64>]) -> 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<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.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::<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![
@ -580,18 +631,12 @@ mod tests {
#[test]
fn zero_loss_test() {
let gram = PartialMatrix({
let mut entries = Vec::<MatrixEntry>::new();
let mut gram = PartialMatrix::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 }
});
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