StudioInfinity/dyna3
2
0
Fork 1
Code Issues 56 Pull requests 1 Projects 3 Releases 4 Wiki Activity Actions

Curvature regulators #80

Merged
glen merged 21 commits from Vectornaut/dyna3:curvature-regulators into main 2025-04-21 23:40:43 +00:00
Conversation 22 Commits 21 Files changed 6 +91 -46
2 changed files with 91 additions and 46 deletions
Download patch file Download diff file Expand all files Collapse all files
Showing only changes of commit d57ff59730 - Show all commits

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.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!(); println!();
let (config, _, success, history) = realize_gram( let (config, _, success, history) = realize_gram(
&problem, 1.0e-12, 0.5, 0.9, 1.1, 200, 110 &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 --- // --- partial matrices ---
struct MatrixEntry { pub struct MatrixEntry {
index: (usize, usize), index: (usize, usize),
value: f64 value: f64
} }
@ -49,42 +49,72 @@ impl PartialMatrix {
PartialMatrix(Vec::<MatrixEntry>::new()) 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; let PartialMatrix(entries) = self;
entries.push(MatrixEntry { index: (row, col), value: value }); 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 { if row != col {
entries.push(MatrixEntry { index: (col, row), value: value }); self.push(col, row, value);
} }
} }
/* DEBUG */ /* DEBUG */
pub fn log_to_console(&self) { pub fn log_to_console(&self) {
let PartialMatrix(entries) = self; for &MatrixEntry { index: (row, col), value } in self {
for ent in entries { console::log_1(&JsValue::from(
let ent_str = format!(" {} {} {}", ent.index.0, ent.index.1, ent.value); format!(" {} {} {}", row, col, value)
console::log_1(&JsValue::from(ent_str.as_str())); ));
} }
} }
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> { fn proj(&self, a: &DMatrix<f64>) -> DMatrix<f64> {
let mut result = DMatrix::<f64>::zeros(a.nrows(), a.ncols()); let mut result = DMatrix::<f64>::zeros(a.nrows(), a.ncols());
let PartialMatrix(entries) = self; for &MatrixEntry { index, .. } in self {
for ent in entries { result[index] = a[index];
result[ent.index] = a[ent.index];
} }
result result
} }
fn sub_proj(&self, rhs: &DMatrix<f64>) -> DMatrix<f64> { fn sub_proj(&self, rhs: &DMatrix<f64>) -> DMatrix<f64> {
let mut result = DMatrix::<f64>::zeros(rhs.nrows(), rhs.ncols()); let mut result = DMatrix::<f64>::zeros(rhs.nrows(), rhs.ncols());
let PartialMatrix(entries) = self; for &MatrixEntry { index, value } in self {
for ent in entries { result[index] = value - rhs[index];
result[ent.index] = ent.value - rhs[ent.index];
} }
result 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 --- // --- configuration subspaces ---
#[derive(Clone)] #[derive(Clone)]
@ -199,8 +229,8 @@ impl DescentHistory {
pub struct ConstraintProblem { pub struct ConstraintProblem {
pub gram: PartialMatrix, pub gram: PartialMatrix,
pub frozen: PartialMatrix,
pub guess: DMatrix<f64>, pub guess: DMatrix<f64>,
pub frozen: Vec<(usize, usize)>
} }
impl ConstraintProblem { impl ConstraintProblem {
@ -208,8 +238,8 @@ impl ConstraintProblem {
const ELEMENT_DIM: usize = 5; const ELEMENT_DIM: usize = 5;
ConstraintProblem { ConstraintProblem {
gram: PartialMatrix::new(), gram: PartialMatrix::new(),
guess: DMatrix::<f64>::zeros(ELEMENT_DIM, element_count), frozen: PartialMatrix::new(),
frozen: Vec::new() guess: DMatrix::<f64>::zeros(ELEMENT_DIM, element_count)
} }
} }
@ -217,8 +247,8 @@ impl ConstraintProblem {
pub fn from_guess(guess_columns: &[DVector<f64>]) -> ConstraintProblem { pub fn from_guess(guess_columns: &[DVector<f64>]) -> ConstraintProblem {
ConstraintProblem { ConstraintProblem {
gram: PartialMatrix::new(), gram: PartialMatrix::new(),
guess: DMatrix::from_columns(guess_columns), frozen: PartialMatrix::new(),
frozen: Vec::new() guess: DMatrix::from_columns(guess_columns)
} }
} }
} }
@ -314,8 +344,10 @@ fn seek_better_config(
None None
} }
// seek a matrix `config` for which `config' * Q * config` matches the partial // seek a matrix `config` that matches the partial matrix `problem.frozen` and
// matrix `gram`. use gradient descent starting from `guess` // 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( pub fn realize_gram(
problem: &ConstraintProblem, problem: &ConstraintProblem,
scaled_tol: f64, scaled_tol: f64,
@ -344,11 +376,11 @@ pub fn realize_gram(
// convert the frozen indices to stacked format // convert the frozen indices to stacked format
let frozen_stacked: Vec<usize> = frozen.into_iter().map( 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(); ).collect();
// use Newton's method with backtracking and gradient descent backup // use a regularized Newton's method with backtracking
let mut state = SearchState::from_config(gram, guess.clone()); let mut state = SearchState::from_config(gram, frozen.freeze(guess));
let mut hess = DMatrix::zeros(element_dim, assembly_dim); let mut hess = DMatrix::zeros(element_dim, assembly_dim);
for _ in 0..max_descent_steps { for _ in 0..max_descent_steps {
// find the negative gradient of the loss function // 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 // the frozen entries fix the radii of the circumscribing sphere, the
// "sun" and "moon" spheres, and one of the chain spheres // "sun" and "moon" spheres, and one of the chain spheres
for k in 0..4 { 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) 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 { 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) realize_gram(&problem, scaled_tol, 0.5, 0.9, 1.1, 200, 110)
@ -559,6 +591,25 @@ mod tests {
use super::{*, examples::*}; 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] #[test]
fn sub_proj_test() { fn sub_proj_test() {
let target = PartialMatrix(vec![ let target = PartialMatrix(vec![
@ -580,18 +631,12 @@ mod tests {
#[test] #[test]
fn zero_loss_test() { fn zero_loss_test() {
let gram = PartialMatrix({ let mut gram = PartialMatrix::new();
let mut entries = Vec::<MatrixEntry>::new();
for j in 0..3 { for j in 0..3 {
for k in 0..3 { for k in 0..3 {
entries.push(MatrixEntry { gram.push(j, k, if j == k { 1.0 } else { -1.0 });
index: (j, k),
value: if j == k { 1.0 } else { -1.0 }
});
} }
} }
entries
});
let config = { let config = {
let a = 0.75_f64.sqrt(); let a = 0.75_f64.sqrt();
DMatrix::from_columns(&[ DMatrix::from_columns(&[
@ -604,33 +649,33 @@ mod tests {
assert!(state.loss.abs() < f64::EPSILON); assert!(state.loss.abs() < f64::EPSILON);
} }
/* TO DO */
// at the frozen indices, the optimization steps should have exact zeros, // 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] #[test]
fn frozen_entry_test() { fn frozen_entry_test() {
let mut problem = ConstraintProblem::from_guess(&[ let mut problem = ConstraintProblem::from_guess(&[
point(0.0, 0.0, 2.0), 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 j in 0..2 {
for k in j..2 { for k in j..2 {
problem.gram.push_sym(j, k, if (j, k) == (1, 1) { 1.0 } else { 0.0 }); 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, 0, problem.guess[(3, 0)]);
problem.frozen.push((3, k)); problem.frozen.push(3, 1, 0.5);
}
let (config, _, success, history) = realize_gram( let (config, _, success, history) = realize_gram(
&problem, 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); assert_eq!(success, true);
for base_step in history.base_step.into_iter() { 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); assert_eq!(base_step[index], 0.0);
} }
} }
for index in problem.frozen { for MatrixEntry { index, value } in problem.frozen {
assert_eq!(config[index], problem.guess[index]); assert_eq!(config[index], value);
} }
} }
@ -663,7 +708,7 @@ mod tests {
} }
} }
for n in 0..ELEMENT_DIM { 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( let (config, tangent, success, history) = realize_gram(
&problem, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110 &problem, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110