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

Compare commits

base: glen:51df7defc52a0efc93f6a805187859aa28c1e51d
Branches Tags
glen:main
glen:study554
glen:pointCoordRegulators
glen:readme-toplevel
glen:kaleidocycle-example
glen:outline-update-fix
glen:lang-trials
glen:seed-problem
StudioInfinity:main
StudioInfinity:readme-toplevel
StudioInfinity:kaleidocycle-example
StudioInfinity:outline-update-fix
StudioInfinity:lang-trials
StudioInfinity:seed-problem
glen:0.3pre
glen:0.2
glen:v0.1.1
glen:v0.1.0
StudioInfinity:0.2
StudioInfinity:0.3pre
StudioInfinity:v0.1.0
StudioInfinity:v0.1.1
..
compare: glen:21cefa9f8a2f297b6cf36fbc31ae5249cca3a88f
Branches Tags
glen:study554
glen:pointCoordRegulators
glen:main
glen:readme-toplevel
glen:kaleidocycle-example
glen:outline-update-fix
glen:lang-trials
glen:seed-problem
StudioInfinity:main
StudioInfinity:readme-toplevel
StudioInfinity:kaleidocycle-example
StudioInfinity:outline-update-fix
StudioInfinity:lang-trials
StudioInfinity:seed-problem
glen:0.3pre
glen:0.2
glen:v0.1.1
glen:v0.1.0
StudioInfinity:0.2
StudioInfinity:0.3pre
StudioInfinity:v0.1.0
StudioInfinity:v0.1.1

No commits in common. "51df7defc52a0efc93f6a805187859aa28c1e51d" and "21cefa9f8a2f297b6cf36fbc31ae5249cca3a88f" have entirely different histories.
51df7defc5 ... 21cefa9f8a

1 changed files with 17 additions and 244 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -490,9 +490,6 @@ pub mod irisawa {
#[cfg(test)]
mod tests {
use nalgebra::Vector3;
use std::{array, f64::consts::{FRAC_1_SQRT_2, PI}, iter};
use super::{*, irisawa::realize_irisawa_hexlet};
#[test]
@ -555,8 +552,10 @@ mod tests {
}
#[test]
fn tangent_test_three_spheres() {
fn tangent_test() {
const SCALED_TOL: f64 = 1.0e-12;
const ELEMENT_DIM: usize = 5;
const ASSEMBLY_DIM: usize = 3;
let gram = {
let mut gram_to_be = PartialMatrix::new();
for j in 0..3 {
@ -581,258 +580,32 @@ mod tests {
assert_eq!(history.scaled_loss.len(), 1);
// confirm that the tangent space has dimension five or less
assert_eq!(tangent.basis_std.len(), 5);
let ConfigSubspace(ref tangent_basis) = tangent;
assert_eq!(tangent_basis.len(), 5);
// 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
const UNIFORM_DIM: usize = 4;
let element_dim = guess.nrows();
let assembly_dim = guess.ncols();
let tangent_motions_unif = vec![
basis_matrix((0, 1), UNIFORM_DIM, assembly_dim),
basis_matrix((1, 1), UNIFORM_DIM, assembly_dim),
basis_matrix((0, 2), UNIFORM_DIM, assembly_dim),
basis_matrix((1, 2), UNIFORM_DIM, assembly_dim),
DMatrix::<f64>::from_column_slice(UNIFORM_DIM, assembly_dim, &[
0.0, 0.0, 0.0, 0.0,
0.0, 0.0, -0.5, -0.5,
0.0, 0.0, -0.5, 0.5
])
];
let tangent_motions_std = vec![
basis_matrix((0, 1), element_dim, assembly_dim),
basis_matrix((1, 1), element_dim, assembly_dim),
basis_matrix((0, 2), element_dim, assembly_dim),
basis_matrix((1, 2), element_dim, assembly_dim),
DMatrix::<f64>::from_column_slice(element_dim, assembly_dim, &[
0.0, 0.0, 0.0, 0.00, 0.0,
let tangent_motions = vec![
basis_matrix((0, 1), ELEMENT_DIM, ASSEMBLY_DIM),
basis_matrix((1, 1), ELEMENT_DIM, ASSEMBLY_DIM),
basis_matrix((0, 2), ELEMENT_DIM, ASSEMBLY_DIM),
basis_matrix((1, 2), ELEMENT_DIM, ASSEMBLY_DIM),
DMatrix::<f64>::from_column_slice(ELEMENT_DIM, 3, &[
0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, -1.0, -0.25, -1.0,
0.0, 0.0, -1.0, 0.25, 1.0
0.0, 0.0, -1.0, 0.25, 1.0
])
];
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(
let tol_sq = ((ELEMENT_DIM * ASSEMBLY_DIM) as f64) * SCALED_TOL * SCALED_TOL;
for motion in tangent_motions {
let motion_proj: DMatrix<_> = motion.column_iter().enumerate().map(
|(k, v)| tangent.proj(&v, k)
).sum();
assert!((motion_std - motion_proj).norm_squared() < tol_sq);
assert!((motion - 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| {
let vel = ang_vel.cross(&pt.fixed_rows::<3>(0));
let mut elt_motion = DVector::zeros(4);
elt_motion.fixed_rows_mut::<3>(0).copy_from(&vel);
elt_motion
}
).collect()
}
#[test]
fn tangent_test_kaleidocycle() {
const SCALED_TOL: f64 = 1.0e-12;
// set up a kaleidocycle, made of points with fixed distances between
// them, and find its tangent space
const N_POINTS: usize = 12;
const N_HINGES: usize = 6;
let gram = {
let mut gram_to_be = PartialMatrix::new();
for block in (0..N_POINTS).step_by(2) {
let block_next = (block + 2) % N_POINTS;
for j in 0..2 {
// diagonal and hinge edges
for k in j..2 {
gram_to_be.push_sym(block + j, block + k, if j == k { 0.0 } else { -0.5 });
}
// non-hinge edges
for k in 0..2 {
gram_to_be.push_sym(block + j, block_next + k, -0.625);
}
}
}
gram_to_be
};
let guess = {
let guess_elts = (0..N_HINGES).step_by(2).flat_map(
|n| {
let ang_hor = (n as f64) * PI/3.0;
let ang_vert = ((n + 1) as f64) * PI/3.0;
let x_vert = ang_vert.cos();
let y_vert = ang_vert.sin();
[
point(0.0, 0.0, 0.0),
point(ang_hor.cos(), ang_hor.sin(), 0.0),
point(x_vert, y_vert, -0.5),
point(x_vert, y_vert, 0.5)
]
}
).collect::<Vec<_>>();
DMatrix::from_columns(&guess_elts)
};
let frozen: [_; N_POINTS] = array::from_fn(|k| (3, k));
let (config, tangent, success, history) = realize_gram(
&gram, guess.clone(), &frozen,
SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
);
assert_eq!(config, guess);
assert_eq!(success, true);
assert_eq!(history.scaled_loss.len(), 1);
// list some motions that should form a basis for the tangent space of
// the solution variety
let element_dim = guess.nrows();
let assembly_dim = guess.ncols();
let tangent_motions_unif = vec![
// the translations along the coordinate axes
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), guess.column_iter().collect()),
rotation_motion_unif(&Vector3::new(0.0, 1.0, 0.0), guess.column_iter().collect()),
rotation_motion_unif(&Vector3::new(0.0, 0.0, 1.0), guess.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:
// - twisting is the kaleidocycle's only internal degree of
// freedom
// - every first-order motion of the kaleidocycle comes from an
// actual motion
(0..N_HINGES).step_by(2).flat_map(
|n| {
let ang_vert = ((n + 1) as f64) * PI/3.0;
let vel_vert_x = 4.0 * ang_vert.cos();
let vel_vert_y = 4.0 * ang_vert.sin();
[
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])
]
}
).collect::<Vec<_>>()
];
let tangent_motions_std = tangent_motions_unif.iter().map(
|motion| DMatrix::from_columns(
&guess.column_iter().zip(motion).map(
|(v, elt_motion)| local_unif_to_std(v) * elt_motion
).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
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);
}
}
fn translation(dis: Vector3<f64>) -> DMatrix<f64> {
const ELEMENT_DIM: usize = 5;
DMatrix::from_column_slice(ELEMENT_DIM, ELEMENT_DIM, &[
1.0, 0.0, 0.0, 0.0, dis[0],
0.0, 1.0, 0.0, 0.0, dis[1],
0.0, 0.0, 1.0, 0.0, dis[2],
2.0*dis[0], 2.0*dis[1], 2.0*dis[2], 1.0, dis.norm_squared(),
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]
fn proj_equivar_test() {
// find a pair of spheres that meet at 120°
const SCALED_TOL: f64 = 1.0e-12;
let gram = {
let mut gram_to_be = PartialMatrix::new();
gram_to_be.push_sym(0, 0, 1.0);
gram_to_be.push_sym(1, 1, 1.0);
gram_to_be.push_sym(0, 1, 0.5);
gram_to_be
};
let guess_orig = DMatrix::from_columns(&[
sphere(0.0, 0.0, 0.5, 1.0),
sphere(0.0, 0.0, -0.5, 1.0)
]);
let (config_orig, tangent_orig, success_orig, history_orig) = realize_gram(
&gram, guess_orig.clone(), &[],
SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
);
assert_eq!(config_orig, guess_orig);
assert_eq!(success_orig, true);
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 = {
let a = 0.5 * FRAC_1_SQRT_2;
DMatrix::from_columns(&[
sphere(a, 0.0, 7.0 + a, 1.0),
sphere(-a, 0.0, 7.0 - a, 1.0)
])
};
let (config_tfm, tangent_tfm, success_tfm, history_tfm) = realize_gram(
&gram, guess_tfm.clone(), &[],
SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
);
assert_eq!(config_tfm, guess_tfm);
assert_eq!(success_tfm, true);
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
const ELEMENT_DIM: usize = 5;
let rot = DMatrix::from_column_slice(ELEMENT_DIM, ELEMENT_DIM, &[
FRAC_1_SQRT_2, 0.0, -FRAC_1_SQRT_2, 0.0, 0.0,
0.0, 1.0, 0.0, 0.0, 0.0,
FRAC_1_SQRT_2, 0.0, FRAC_1_SQRT_2, 0.0, 0.0,
0.0, 0.0, 0.0, 1.0, 0.0,
0.0, 0.0, 0.0, 0.0, 1.0
]);
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
const SCALED_TOL_TFM: f64 = 1.0e-9;
let tol_sq = ((guess_orig.nrows() * guess_orig.ncols()) as f64) * SCALED_TOL_TFM * SCALED_TOL_TFM;
assert!((motion_proj_tfm - motion_tfm_proj).norm_squared() < tol_sq);
}
// at the frozen indices, the optimization steps should have exact zeros,
// and the realized configuration should match the initial guess
#[test]