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

chore: wrap lines at 80 characters
All checks were successful
/ test (pull_request) Successful in 3m36s
Details

Browse source
This commit is contained in:
Glen Whitney 2025-10-10 22:34:18 -07:00
parent a4101ab81c
commit 19e5458e1b
7 changed files with 310 additions and 150 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -9,8 +9,11 @@ pub fn point(x: f64, y: f64, z: f64) -> DVector<f64> {
}
// the sphere with the given center and radius, with inward-pointing normals
pub fn sphere(center_x: f64, center_y: f64, center_z: f64, radius: f64) -> DVector<f64> {
let center_norm_sq = center_x * center_x + center_y * center_y + center_z * center_z;
pub fn sphere(center_x: f64, center_y: f64, center_z: f64, radius: f64)
-> DVector<f64>
{
let center_norm_sq =
center_x * center_x + center_y * center_y + center_z * center_z;
DVector::from_column_slice(&[
center_x / radius,
center_y / radius,
@ -23,7 +26,9 @@ pub fn sphere(center_x: f64, center_y: f64, center_z: f64, radius: f64) -> DVect
// the sphere of curvature `curv` whose closest point to the origin has position
// `off * dir` and normal `dir`, where `dir` is a unit vector. setting the
// curvature to zero gives a plane
pub fn sphere_with_offset(dir_x: f64, dir_y: f64, dir_z: f64, off: f64, curv: f64) -> DVector<f64> {
pub fn sphere_with_offset(
dir_x: f64, dir_y: f64, dir_z: f64, off: f64, curv: f64) -> DVector<f64>
{
let norm_sp = 1.0 + off * curv;
DVector::from_column_slice(&[
norm_sp * dir_x,
@ -200,7 +205,9 @@ impl ConfigSubspace {
// with the given column index has velocity `v`. the velocity is given in
// projection coordinates, and the projection is done with respect to the
// projection inner product
pub fn proj(&self, v: &DVectorView<f64>, column_index: usize) -> DMatrix<f64> {
pub fn proj(&self, v: &DVectorView<f64>, column_index: usize)
-> DMatrix<f64>
{
if self.dim() == 0 {
const ELEMENT_DIM: usize = 5;
DMatrix::zeros(ELEMENT_DIM, self.assembly_dim)
@ -291,7 +298,9 @@ impl SearchState {
}
}
fn basis_matrix(index: (usize, usize), nrows: usize, ncols: usize) -> DMatrix<f64> {
fn basis_matrix(index: (usize, usize), nrows: usize, ncols: usize)
-> DMatrix<f64>
{
let mut result = DMatrix::<f64>::zeros(nrows, ncols);
result[index] = 1.0;
result
@ -414,7 +423,8 @@ pub fn realize_gram(
for _ in 0..max_descent_steps {
// find the negative gradient of the loss function
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>);
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
@ -431,7 +441,8 @@ pub fn realize_gram(
-&basis_mat * &state.err_proj
+ &state.config * &neg_d_err_proj
);
hess_cols.push(deriv_grad.reshape_generic(Dyn(total_dim), Const::<1>));
hess_cols.push(
deriv_grad.reshape_generic(Dyn(total_dim), Const::<1>));
}
}
hess = DMatrix::from_columns(hess_cols.as_slice());
@ -440,7 +451,8 @@ pub fn realize_gram(
let hess_eigvals = hess.symmetric_eigenvalues();
let min_eigval = hess_eigvals.min();
if min_eigval <= 0.0 {
hess -= reg_scale * min_eigval * DMatrix::identity(total_dim, total_dim);
hess -= reg_scale * min_eigval
* DMatrix::identity(total_dim, total_dim);
}
history.hess_eigvals.push(hess_eigvals);
@ -477,7 +489,8 @@ 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));
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
@ -507,9 +520,12 @@ pub fn realize_gram(
}
// find the kernel of the Hessian. give it the uniform inner product
let tangent = ConfigSubspace::symmetric_kernel(hess, unif_to_std, assembly_dim);
let tangent =
ConfigSubspace::symmetric_kernel(hess, unif_to_std, assembly_dim);
Ok(ConfigNeighborhood { #[cfg(feature = "dev")] config: state.config, nbhd: tangent })
Ok(ConfigNeighborhood {
#[cfg(feature = "dev")] config: state.config, nbhd: tangent
})
} else {
Err("Failed to reach target accuracy".to_string())
};
@ -608,7 +624,8 @@ pub mod examples {
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 });
problem.gram.push_sym(
block + j, block + k, if j == k { 0.0 } else { -0.5 });
}
// non-hinge edges
@ -702,7 +719,8 @@ mod tests {
]);
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.gram.push_sym(
j, k, if (j, k) == (1, 1) { 1.0 } else { 0.0 });
}
}
problem.frozen.push(3, 0, problem.guess[(3, 0)]);
@ -729,7 +747,8 @@ mod tests {
// 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];
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];
for (k, diam) in solution_diams.into_iter().enumerate() {
assert!((config[(3, k)] - 1.0 / diam).abs() < entry_tol);
}
@ -794,22 +813,29 @@ mod tests {
// 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
let tol_sq = ((element_dim * assembly_dim) as f64) * SCALED_TOL * SCALED_TOL;
for (motion_unif, motion_std) in tangent_motions_unif.into_iter().zip(tangent_motions_std) {
let motion_proj: DMatrix<_> = motion_unif.column_iter().enumerate().map(
|(k, v)| tangent.proj(&v, k)
).sum();
let tol_sq = ((element_dim * assembly_dim) as f64)
* SCALED_TOL * SCALED_TOL;
for (motion_unif, motion_std)
in tangent_motions_unif.into_iter().zip(tangent_motions_std) {
let motion_proj: DMatrix<_> =
motion_unif.column_iter().enumerate().map(
|(k, v)| tangent.proj(&v, k)
).sum();
assert!((motion_std - motion_proj).norm_squared() < tol_sq);
}
}
fn translation_motion_unif(vel: &Vector3<f64>, assembly_dim: usize) -> Vec<DVector<f64>> {
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>> {
fn rotation_motion_unif(
ang_vel: &Vector3<f64>, points: Vec<DVectorView<f64>>
) -> Vec<DVector<f64>> {
points.into_iter().map(
|pt| {
let vel = ang_vel.cross(&pt.fixed_rows::<3>(0));
@ -840,9 +866,12 @@ mod tests {
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()),
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
@ -859,8 +888,10 @@ mod tests {
[
DVector::from_column_slice(&[0.0, 0.0, 5.0, 0.0]),
DVector::from_column_slice(&[0.0, 0.0, 1.0, 0.0]),
DVector::from_column_slice(&[-vel_vert_x, -vel_vert_y, -3.0, 0.0]),
DVector::from_column_slice(&[vel_vert_x, vel_vert_y, -3.0, 0.0]),
DVector::from_column_slice(
&[-vel_vert_x, -vel_vert_y, -3.0, 0.0]),
DVector::from_column_slice(
&[vel_vert_x, vel_vert_y, -3.0, 0.0]),
]
}
).collect::<Vec<_>>(),
@ -880,11 +911,14 @@ mod tests {
// 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
let tol_sq = ((element_dim * assembly_dim) as f64) * SCALED_TOL * SCALED_TOL;
for (motion_unif, motion_std) in tangent_motions_unif.into_iter().zip(tangent_motions_std) {
let motion_proj: DMatrix<_> = motion_unif.into_iter().enumerate().map(
|(k, v)| tangent.proj(&v.as_view(), k)
).sum();
let tol_sq = ((element_dim * assembly_dim) as f64)
* SCALED_TOL * SCALED_TOL;
for (motion_unif, motion_std)
in tangent_motions_unif.into_iter().zip(tangent_motions_std) {
let motion_proj: DMatrix<_> =
motion_unif.into_iter().enumerate().map(
|(k, v)| tangent.proj(&v.as_view(), k)
).sum();
assert!((motion_std - motion_proj).norm_squared() < tol_sq);
}
}
@ -913,10 +947,10 @@ mod tests {
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 Realization { result: result_orig, history: history_orig } = realize_gram(
&problem_orig, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
);
let ConfigNeighborhood { config: config_orig, nbhd: tangent_orig } = result_orig.unwrap();
let Realization { result: result_orig, history: history_orig } =
realize_gram(&problem_orig, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110);
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);
@ -934,10 +968,10 @@ mod tests {
frozen: problem_orig.frozen,
guess: guess_tfm,
};
let Realization { result: result_tfm, history: history_tfm } = realize_gram(
&problem_tfm, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
);
let ConfigNeighborhood { config: config_tfm, nbhd: tangent_tfm } = result_tfm.unwrap();
let Realization { result: result_tfm, history: history_tfm } =
realize_gram(&problem_tfm, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110);
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);
@ -948,7 +982,8 @@ mod tests {
// 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 = 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
@ -969,7 +1004,9 @@ 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 = ((problem_orig.guess.nrows() * problem_orig.guess.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);
}
}