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

chore: remove trailing whitespace, add CR at end of file
All checks were successful
/ test (pull_request) Successful in 3m38s
Details

Browse source
This commit is contained in:
Glen Whitney 2025-10-10 10:20:38 -07:00
parent 6dbbe2ce2d
commit a4101ab81c
11 changed files with 320 additions and 320 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -62,19 +62,19 @@ impl PartialMatrix {
pub fn new() -> Self {
Self(Vec::<MatrixEntry>::new())
}
pub fn push(&mut self, row: usize, col: usize, value: f64) {
let Self(entries) = self;
entries.push(MatrixEntry { index: (row, col), value });
}
pub fn push_sym(&mut self, row: usize, col: usize, value: f64) {
self.push(row, col, value);
if row != col {
self.push(col, row, value);
}
}
fn freeze(&self, a: &DMatrix<f64>) -> DMatrix<f64> {
let mut result = a.clone();
for &MatrixEntry { index, value } in self {
@ -82,7 +82,7 @@ impl PartialMatrix {
}
result
}
fn proj(&self, a: &DMatrix<f64>) -> DMatrix<f64> {
let mut result = DMatrix::<f64>::zeros(a.nrows(), a.ncols());
for &MatrixEntry { index, .. } in self {
@ -90,7 +90,7 @@ impl PartialMatrix {
}
result
}
fn sub_proj(&self, rhs: &DMatrix<f64>) -> DMatrix<f64> {
let mut result = DMatrix::<f64>::zeros(rhs.nrows(), rhs.ncols());
for &MatrixEntry { index, value } in self {
@ -112,7 +112,7 @@ impl Display for PartialMatrix {
impl IntoIterator for PartialMatrix {
type Item = MatrixEntry;
type IntoIter = std::vec::IntoIter<Self::Item>;
fn into_iter(self) -> Self::IntoIter {
let Self(entries) = self;
entries.into_iter()
@ -122,7 +122,7 @@ impl IntoIterator for PartialMatrix {
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()
@ -146,7 +146,7 @@ impl ConfigSubspace {
basis_std: Vec::new(),
}
}
// approximate the kernel of a symmetric endomorphism of the configuration
// space for `assembly_dim` elements. we consider an eigenvector to be part
// of the kernel if its eigenvalue is smaller than the constant `THRESHOLD`
@ -167,10 +167,10 @@ impl ConfigSubspace {
|(λ, v)| (λ.abs() < THRESHOLD).then_some(v)
).collect::<Vec<_>>().as_slice()
);
// express the basis in the standard coordinates
let basis_std = proj_to_std * &basis_proj;
const ELEMENT_DIM: usize = 5;
const UNIFORM_DIM: usize = 4;
Self {
@ -187,15 +187,15 @@ impl ConfigSubspace {
).collect(),
}
}
pub fn dim(&self) -> usize {
self.basis_std.len()
}
pub fn assembly_dim(&self) -> usize {
self.assembly_dim
}
// find the projection onto this subspace of the motion where the element
// with the given column index has velocity `v`. the velocity is given in
// projection coordinates, and the projection is done with respect to the
@ -253,7 +253,7 @@ impl ConstraintProblem {
guess: DMatrix::<f64>::zeros(ELEMENT_DIM, element_count),
}
}
#[cfg(feature = "dev")]
pub fn from_guess(guess_columns: &[DVector<f64>]) -> Self {
Self {
@ -377,10 +377,10 @@ pub fn realize_gram(
) -> Realization {
// destructure the problem data
let ConstraintProblem { gram, guess, frozen } = problem;
// start the descent history
let mut history = DescentHistory::new();
// handle the case where the assembly is empty. our general realization
// routine can't handle this case because it builds the Hessian using
// `DMatrix::from_columns`, which panics when the list of columns is empty
@ -394,20 +394,20 @@ pub fn realize_gram(
);
return Realization { result, history };
}
// find the dimension of the search space
let element_dim = guess.nrows();
let total_dim = element_dim * assembly_dim;
// scale the tolerance
let scale_adjustment = (gram.0.len() as f64).sqrt();
let tol = scale_adjustment * scaled_tol;
// convert the frozen indices to stacked format
let frozen_stacked: Vec<usize> = frozen.into_iter().map(
|MatrixEntry { index: (row, col), .. }| col*element_dim + row
).collect();
// 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);
@ -416,7 +416,7 @@ pub fn realize_gram(
let neg_grad = 4.0 * &*Q * &state.config * &state.err_proj;
let mut neg_grad_stacked = neg_grad.clone().reshape_generic(Dyn(total_dim), Const::<1>);
history.neg_grad.push(neg_grad.clone());
// find the negative Hessian of the loss function
let mut hess_cols = Vec::<DVector<f64>>::with_capacity(total_dim);
for col in 0..assembly_dim {
@ -435,7 +435,7 @@ pub fn realize_gram(
}
}
hess = DMatrix::from_columns(hess_cols.as_slice());
// regularize the Hessian
let hess_eigvals = hess.symmetric_eigenvalues();
let min_eigval = hess_eigvals.min();
@ -443,7 +443,7 @@ pub fn realize_gram(
hess -= reg_scale * min_eigval * DMatrix::identity(total_dim, total_dim);
}
history.hess_eigvals.push(hess_eigvals);
// project the negative gradient and negative Hessian onto the
// orthogonal complement of the frozen subspace
let zero_col = DVector::zeros(total_dim);
@ -454,12 +454,12 @@ pub fn realize_gram(
hess.set_column(k, &zero_col);
hess[(k, k)] = 1.0;
}
// stop if the loss is tolerably low
history.config.push(state.config.clone());
history.scaled_loss.push(state.loss / scale_adjustment);
if state.loss < tol { break; }
// compute the Newton step
/* TO DO */
/*
@ -479,7 +479,7 @@ pub fn realize_gram(
let base_step_stacked = hess_cholesky.solve(&neg_grad_stacked);
let base_step = base_step_stacked.reshape_generic(Dyn(element_dim), Dyn(assembly_dim));
history.base_step.push(base_step.clone());
// use backtracking line search to find a better configuration
if let Some((better_state, backoff_steps)) = seek_better_config(
gram, &state, &base_step, neg_grad.dot(&base_step),
@ -505,10 +505,10 @@ pub fn realize_gram(
.view_mut(block_start, (element_dim, UNIFORM_DIM))
.copy_from(&local_unif_to_std(state.config.column(n)));
}
// find the kernel of the Hessian. give it the uniform inner product
let tangent = ConfigSubspace::symmetric_kernel(hess, unif_to_std, assembly_dim);
Ok(ConfigNeighborhood { #[cfg(feature = "dev")] config: state.config, nbhd: tangent })
} else {
Err("Failed to reach target accuracy".to_string())
@ -521,9 +521,9 @@ pub fn realize_gram(
#[cfg(feature = "dev")]
pub mod examples {
use std::f64::consts::PI;
use super::*;
// this problem is from a sangaku by Irisawa Shintarō Hiroatsu. the article
// below includes a nice translation of the problem statement, which was
// recorded in Uchida Itsumi's book _Kokon sankan_ (_Mathematics, Past and
@ -547,40 +547,40 @@ 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
for k in 0..4 {
problem.frozen.push(3, k, problem.guess[(3, k)]);
}
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) -> Realization {
@ -601,7 +601,7 @@ pub mod examples {
}
).collect::<Vec<_>>().as_slice()
);
const N_POINTS: usize = 2 * N_HINGES;
for block in (0..N_POINTS).step_by(2) {
let block_next = (block + 2) % N_POINTS;
@ -610,18 +610,18 @@ pub mod examples {
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);
}
}
}
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)
}
}
@ -630,9 +630,9 @@ pub mod examples {
mod tests {
use nalgebra::Vector3;
use std::{f64::consts::{FRAC_1_SQRT_2, PI}, iter};
use super::{*, examples::*};
#[test]
fn freeze_test() {
let frozen = PartialMatrix(vec![
@ -651,7 +651,7 @@ mod tests {
]);
assert_eq!(frozen.freeze(&config), expected_result);
}
#[test]
fn sub_proj_test() {
let target = PartialMatrix(vec![
@ -670,7 +670,7 @@ mod tests {
]);
assert_eq!(target.sub_proj(&attempt), expected_result);
}
#[test]
fn zero_loss_test() {
let mut gram = PartialMatrix::new();
@ -690,7 +690,7 @@ mod tests {
let state = SearchState::from_config(&gram, config);
assert!(state.loss.abs() < f64::EPSILON);
}
/* TO DO */
// at the frozen indices, the optimization steps should have exact zeros,
// and the realized configuration should have the desired values
@ -720,13 +720,13 @@ mod tests {
assert_eq!(config[index], value);
}
}
#[test]
fn irisawa_hexlet_test() {
// solve Irisawa's problem
const SCALED_TOL: f64 = 1.0e-12;
let config = realize_irisawa_hexlet(SCALED_TOL).result.unwrap().config;
// check against Irisawa's solution
let entry_tol = SCALED_TOL.sqrt();
let solution_diams = [30.0, 10.0, 6.0, 5.0, 15.0, 10.0, 3.75, 2.5, 2.0 + 8.0/11.0];
@ -734,7 +734,7 @@ mod tests {
assert!((config[(3, k)] - 1.0 / diam).abs() < entry_tol);
}
}
#[test]
fn tangent_test_three_spheres() {
const SCALED_TOL: f64 = 1.0e-12;
@ -758,7 +758,7 @@ mod tests {
let ConfigNeighborhood { config, nbhd: tangent } = result.unwrap();
assert_eq!(config, problem.guess);
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;
@ -786,11 +786,11 @@ mod tests {
0.0, 0.0, -1.0, 0.25, 1.0,
]),
];
// confirm that the dimension of the tangent space is no greater than
// expected
assert_eq!(tangent.basis_std.len(), tangent_motions_std.len());
// confirm that the tangent space contains all the motions we expect it
// to. since we've already bounded the dimension of the tangent space,
// this confirms that the tangent space is what we expect it to be
@ -802,13 +802,13 @@ mod tests {
assert!((motion_std - motion_proj).norm_squared() < tol_sq);
}
}
fn translation_motion_unif(vel: &Vector3<f64>, assembly_dim: usize) -> Vec<DVector<f64>> {
let mut elt_motion = DVector::zeros(4);
elt_motion.fixed_rows_mut::<3>(0).copy_from(vel);
iter::repeat(elt_motion).take(assembly_dim).collect()
}
fn rotation_motion_unif(ang_vel: &Vector3<f64>, points: Vec<DVectorView<f64>>) -> Vec<DVector<f64>> {
points.into_iter().map(
|pt| {
@ -819,7 +819,7 @@ mod tests {
}
).collect()
}
#[test]
fn tangent_test_kaleidocycle() {
// set up a kaleidocycle and find its tangent space
@ -827,7 +827,7 @@ mod tests {
let Realization { result, history } = realize_kaleidocycle(SCALED_TOL);
let ConfigNeighborhood { config, nbhd: tangent } = result.unwrap();
assert_eq!(history.scaled_loss.len(), 1);
// list some motions that should form a basis for the tangent space of
// the solution variety
const N_HINGES: usize = 6;
@ -838,12 +838,12 @@ mod tests {
translation_motion_unif(&Vector3::new(1.0, 0.0, 0.0), assembly_dim),
translation_motion_unif(&Vector3::new(0.0, 1.0, 0.0), assembly_dim),
translation_motion_unif(&Vector3::new(0.0, 0.0, 1.0), assembly_dim),
// the rotations about the coordinate axes
rotation_motion_unif(&Vector3::new(1.0, 0.0, 0.0), config.column_iter().collect()),
rotation_motion_unif(&Vector3::new(0.0, 1.0, 0.0), config.column_iter().collect()),
rotation_motion_unif(&Vector3::new(0.0, 0.0, 1.0), config.column_iter().collect()),
// the twist motion. more precisely: a motion that keeps the center
// of mass stationary and preserves the distances between the
// vertices to first order. this has to be the twist as long as:
@ -872,11 +872,11 @@ mod tests {
).collect::<Vec<_>>()
)
).collect::<Vec<_>>();
// confirm that the dimension of the tangent space is no greater than
// expected
assert_eq!(tangent.basis_std.len(), tangent_motions_unif.len());
// confirm that the tangent space contains all the motions we expect it
// to. since we've already bounded the dimension of the tangent space,
// this confirms that the tangent space is what we expect it to be
@ -888,7 +888,7 @@ mod tests {
assert!((motion_std - motion_proj).norm_squared() < tol_sq);
}
}
fn translation(dis: Vector3<f64>) -> DMatrix<f64> {
const ELEMENT_DIM: usize = 5;
DMatrix::from_column_slice(ELEMENT_DIM, ELEMENT_DIM, &[
@ -899,7 +899,7 @@ mod tests {
0.0, 0.0, 0.0, 0.0, 1.0,
])
}
// confirm that projection onto a configuration subspace is equivariant with
// respect to Euclidean motions
#[test]
@ -919,7 +919,7 @@ mod tests {
let ConfigNeighborhood { config: config_orig, nbhd: tangent_orig } = result_orig.unwrap();
assert_eq!(config_orig, problem_orig.guess);
assert_eq!(history_orig.scaled_loss.len(), 1);
// find another pair of spheres that meet at 120°. we'll think of this
// solution as a transformed version of the original one
let guess_tfm = {
@ -940,17 +940,17 @@ mod tests {
let ConfigNeighborhood { config: config_tfm, nbhd: tangent_tfm } = result_tfm.unwrap();
assert_eq!(config_tfm, problem_tfm.guess);
assert_eq!(history_tfm.scaled_loss.len(), 1);
// project a nudge to the tangent space of the solution variety at the
// original solution
let motion_orig = DVector::from_column_slice(&[0.0, 0.0, 1.0, 0.0]);
let motion_orig_proj = tangent_orig.proj(&motion_orig.as_view(), 0);
// project the equivalent nudge to the tangent space of the solution
// variety at the transformed solution
let motion_tfm = DVector::from_column_slice(&[FRAC_1_SQRT_2, 0.0, FRAC_1_SQRT_2, 0.0]);
let motion_tfm_proj = tangent_tfm.proj(&motion_tfm.as_view(), 0);
// take the transformation that sends the original solution to the
// transformed solution and apply it to the motion that the original
// solution makes in response to the nudge
@ -964,7 +964,7 @@ mod tests {
]);
let transl = translation(Vector3::new(0.0, 0.0, 7.0));
let motion_proj_tfm = transl * rot * motion_orig_proj;
// confirm that the projection of the nudge is equivariant. we loosen
// the comparison tolerance because the transformation seems to
// introduce some numerical error
@ -972,4 +972,4 @@ mod tests {
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);
}
}
}