Switch to Euclidean-invariant projection onto tangent space of solution variety (#34)
This pull request addresses issues #32 and #33 by projecting nudges onto the tangent space of the solution variety using a Euclidean-invariant inner product, which I'm calling the *uniform* inner product. ### Definition of the uniform inner product For spheres and planes, the uniform inner product is defined on the tangent space of the hyperboloid $\langle v, v \rangle = 1$. For points, it's defined on the tangent space of the paraboloid $\langle v, v \rangle = 0,\; \langle v, I_\infty \rangle = 1$. The tangent space of an assembly can be expressed as the direct sum of the tangent spaces of the elements. We extend the uniform inner product to assemblies by declaring the tangent spaces of different elements to be orthogonal. #### For spheres and planes If $v = [x, y, z, b, c]^\top$ is on the hyperboloid $\langle v, v \rangle = 1$, the vectors $$\left[ \begin{array}{c} 2b \\ \cdot \\ \cdot \\ \cdot \\ x \end{array} \right],\;\left[ \begin{array}{c} \cdot \\ 2b \\ \cdot \\ \cdot \\ y \end{array} \right],\;\left[ \begin{array}{c} \cdot \\ \cdot \\ 2b \\ \cdot \\ z \end{array} \right],\;\left[ \begin{array}{l} 2bx \\ 2by \\ 2bz \\ 2b^2 \\ 2bc + 1 \end{array} \right]$$ form a basis for the tangent space of hyperboloid at $v$. We declare this basis to be orthonormal with respect to the uniform inner product. The first three vectors in the basis are unit-speed translations along the coordinate axes. The last vector moves the surface at unit speed along its normal field. For spheres, this increases the radius at unit rate. For planes, this translates the plane parallel to itself at unit speed. This description makes it clear that the uniform inner product is invariant under Euclidean motions. #### For points If $v = [x, y, z, b, c]^\top$ is on the paraboloid $\langle v, v \rangle = 0,\; \langle v, I_\infty \rangle = 1$, the vectors $$\left[ \begin{array}{c} 2b \\ \cdot \\ \cdot \\ \cdot \\ x \end{array} \right],\;\left[ \begin{array}{c} \cdot \\ 2b \\ \cdot \\ \cdot \\ y \end{array} \right],\;\left[ \begin{array}{c} \cdot \\ \cdot \\ 2b \\ \cdot \\ z \end{array} \right]$$ form a basis for the tangent space of paraboloid at $v$. We declare this basis to be orthonormal with respect to the uniform inner product. The meanings of the basis vectors, and the argument that the uniform inner product is Euclidean-invariant, are the same as for spheres and planes. In the engine, we pad the basis with $[0, 0, 0, 0, 1]^\top$ to keep the number of uniform coordinates consistent across element types. ### Confirmation of intended behavior Two new tests confirm that we've corrected the misbehaviors described in issues #32 and #33. Issue | Test ---|--- #32 | `proj_equivar_test` #33 | `tangent_test_kaleidocycle` Co-authored-by: Aaron Fenyes <aaron.fenyes@fareycircles.ooo> Reviewed-on: #34 Co-authored-by: Vectornaut <vectornaut@nobody@nowhere.net> Co-committed-by: Vectornaut <vectornaut@nobody@nowhere.net>
This commit is contained in:
parent
22870342f3
commit
817a446fad
72
app-proto/examples/kaleidocycle.rs
Normal file
72
app-proto/examples/kaleidocycle.rs
Normal file
@ -0,0 +1,72 @@
|
||||
use nalgebra::{DMatrix, DVector};
|
||||
use std::{array, f64::consts::PI};
|
||||
|
||||
use dyna3::engine::{Q, point, realize_gram, PartialMatrix};
|
||||
|
||||
fn main() {
|
||||
// set up a kaleidocycle, made of points with fixed distances between them,
|
||||
// and find its tangent space
|
||||
const N_POINTS: usize = 12;
|
||||
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 = {
|
||||
const N_HINGES: usize = 6;
|
||||
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, &frozen,
|
||||
1.0e-12, 0.5, 0.9, 1.1, 200, 110
|
||||
);
|
||||
print!("Completed Gram matrix:{}", config.tr_mul(&*Q) * &config);
|
||||
print!("Configuration:{}", config);
|
||||
if success {
|
||||
println!("Target accuracy achieved!");
|
||||
} else {
|
||||
println!("Failed to reach target accuracy");
|
||||
}
|
||||
println!("Steps: {}", history.scaled_loss.len() - 1);
|
||||
println!("Loss: {}\n", history.scaled_loss.last().unwrap());
|
||||
|
||||
// find the kaleidocycle's twist motion
|
||||
let up = DVector::from_column_slice(&[0.0, 0.0, 1.0, 0.0]);
|
||||
let down = -&up;
|
||||
let twist_motion: DMatrix<_> = (0..N_POINTS).step_by(4).flat_map(
|
||||
|n| [
|
||||
tangent.proj(&up.as_view(), n),
|
||||
tangent.proj(&down.as_view(), n+1)
|
||||
]
|
||||
).sum();
|
||||
let normalization = 5.0 / twist_motion[(2, 0)];
|
||||
print!("Twist motion:{}", normalization * twist_motion);
|
||||
}
|
||||
@ -9,3 +9,4 @@
|
||||
cargo run --example irisawa-hexlet
|
||||
cargo run --example three-spheres
|
||||
cargo run --example point-on-sphere
|
||||
cargo run --example kaleidocycle
|
||||
@ -5,7 +5,7 @@ use std::{collections::BTreeSet, sync::atomic::{AtomicU64, Ordering}};
|
||||
use sycamore::prelude::*;
|
||||
use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */
|
||||
|
||||
use crate::engine::{realize_gram, ConfigSubspace, PartialMatrix, Q};
|
||||
use crate::engine::{realize_gram, local_unif_to_std, ConfigSubspace, PartialMatrix, Q};
|
||||
|
||||
// the types of the keys we use to access an assembly's elements and constraints
|
||||
pub type ElementKey = usize;
|
||||
@ -120,6 +120,7 @@ pub struct Constraint {
|
||||
pub active: Signal<bool>
|
||||
}
|
||||
|
||||
// the velocity is expressed in uniform coordinates
|
||||
pub struct ElementMotion<'a> {
|
||||
pub key: ElementKey,
|
||||
pub velocity: DVectorView<'a, f64>
|
||||
@ -359,7 +360,14 @@ impl Assembly {
|
||||
// this element didn't have a column index when we started, so
|
||||
// by invariant (2), it's unconstrained
|
||||
let mut target_column = motion_proj.column_mut(column_index);
|
||||
target_column += elt_motion.velocity;
|
||||
let unif_to_std = self.elements.with_untracked(
|
||||
|elts| {
|
||||
elts[elt_motion.key].representation.with_untracked(
|
||||
|rep| local_unif_to_std(rep.as_view())
|
||||
)
|
||||
}
|
||||
);
|
||||
target_column += unif_to_std * elt_motion.velocity;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -130,6 +130,8 @@ pub fn Display() -> View {
|
||||
let translate_pos_y = create_signal(0.0);
|
||||
let translate_neg_z = create_signal(0.0);
|
||||
let translate_pos_z = create_signal(0.0);
|
||||
let shrink_neg = create_signal(0.0);
|
||||
let shrink_pos = create_signal(0.0);
|
||||
|
||||
// change listener
|
||||
let scene_changed = create_signal(true);
|
||||
@ -164,6 +166,7 @@ pub fn Display() -> View {
|
||||
|
||||
// manipulation
|
||||
const TRANSLATION_SPEED: f64 = 0.15; // in length units per second
|
||||
const SHRINKING_SPEED: f64 = 0.15; // in length units per second
|
||||
|
||||
// display parameters
|
||||
const OPACITY: f32 = 0.5; /* SCAFFOLDING */
|
||||
@ -292,6 +295,8 @@ pub fn Display() -> View {
|
||||
let translate_pos_y_val = translate_pos_y.get();
|
||||
let translate_neg_z_val = translate_neg_z.get();
|
||||
let translate_pos_z_val = translate_pos_z.get();
|
||||
let shrink_neg_val = shrink_neg.get();
|
||||
let shrink_pos_val = shrink_pos.get();
|
||||
|
||||
// update the assembly's orientation
|
||||
let ang_vel = {
|
||||
@ -323,24 +328,27 @@ pub fn Display() -> View {
|
||||
let sel = state.selection.with(
|
||||
|sel| *sel.into_iter().next().unwrap()
|
||||
);
|
||||
let rep = state.assembly.elements.with_untracked(
|
||||
|elts| elts[sel].representation.get_clone_untracked()
|
||||
);
|
||||
let translate_x = translate_pos_x_val - translate_neg_x_val;
|
||||
let translate_y = translate_pos_y_val - translate_neg_y_val;
|
||||
let translate_z = translate_pos_z_val - translate_neg_z_val;
|
||||
if translate_x != 0.0 || translate_y != 0.0 || translate_z != 0.0 {
|
||||
let vel_field = {
|
||||
let u = Vector3::new(translate_x, translate_y, translate_z).normalize();
|
||||
DMatrix::from_column_slice(5, 5, &[
|
||||
0.0, 0.0, 0.0, 0.0, u[0],
|
||||
0.0, 0.0, 0.0, 0.0, u[1],
|
||||
0.0, 0.0, 0.0, 0.0, u[2],
|
||||
2.0*u[0], 2.0*u[1], 2.0*u[2], 0.0, 0.0,
|
||||
0.0, 0.0, 0.0, 0.0, 0.0
|
||||
])
|
||||
let shrink = shrink_pos_val - shrink_neg_val;
|
||||
let translating =
|
||||
translate_x != 0.0
|
||||
|| translate_y != 0.0
|
||||
|| translate_z != 0.0;
|
||||
if translating || shrink != 0.0 {
|
||||
let elt_motion = {
|
||||
let u = if translating {
|
||||
TRANSLATION_SPEED * Vector3::new(
|
||||
translate_x, translate_y, translate_z
|
||||
).normalize()
|
||||
} else {
|
||||
Vector3::zeros()
|
||||
};
|
||||
time_step * DVector::from_column_slice(
|
||||
&[u[0], u[1], u[2], SHRINKING_SPEED * shrink]
|
||||
)
|
||||
};
|
||||
let elt_motion: DVector<f64> = time_step * TRANSLATION_SPEED * vel_field * rep;
|
||||
assembly_for_raf.deform(
|
||||
vec![
|
||||
ElementMotion {
|
||||
@ -501,6 +509,8 @@ pub fn Display() -> View {
|
||||
"s" | "S" if shift => translate_pos_z.set(value),
|
||||
"w" | "W" => translate_pos_y.set(value),
|
||||
"s" | "S" => translate_neg_y.set(value),
|
||||
"]" | "}" => shrink_neg.set(value),
|
||||
"[" | "{" => shrink_pos.set(value),
|
||||
_ => manipulating = false
|
||||
};
|
||||
if manipulating {
|
||||
|
||||
@ -90,32 +90,34 @@ impl PartialMatrix {
|
||||
#[derive(Clone)]
|
||||
pub struct ConfigSubspace {
|
||||
assembly_dim: usize,
|
||||
basis: Vec<DMatrix<f64>>
|
||||
basis_std: Vec<DMatrix<f64>>,
|
||||
basis_proj: Vec<DMatrix<f64>>
|
||||
}
|
||||
|
||||
impl ConfigSubspace {
|
||||
pub fn zero(assembly_dim: usize) -> ConfigSubspace {
|
||||
ConfigSubspace {
|
||||
assembly_dim: assembly_dim,
|
||||
basis: Vec::new()
|
||||
basis_proj: Vec::new(),
|
||||
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`
|
||||
fn symmetric_kernel(a: DMatrix<f64>, assembly_dim: usize) -> ConfigSubspace {
|
||||
const ELEMENT_DIM: usize = 5;
|
||||
const THRESHOLD: f64 = 1.0e-4;
|
||||
let eig = SymmetricEigen::new(a);
|
||||
fn symmetric_kernel(a: DMatrix<f64>, proj_to_std: DMatrix<f64>, assembly_dim: usize) -> ConfigSubspace {
|
||||
// find a basis for the kernel. the basis is expressed in the projection
|
||||
// coordinates, and it's orthonormal with respect to the projection
|
||||
// inner product
|
||||
const THRESHOLD: f64 = 0.1;
|
||||
let eig = SymmetricEigen::new(proj_to_std.tr_mul(&a) * &proj_to_std);
|
||||
let eig_vecs = eig.eigenvectors.column_iter();
|
||||
let eig_pairs = eig.eigenvalues.iter().zip(eig_vecs);
|
||||
let basis = eig_pairs.filter_map(
|
||||
|(λ, v)| (λ.abs() < THRESHOLD).then_some(
|
||||
Into::<DMatrix<f64>>::into(
|
||||
v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim))
|
||||
)
|
||||
)
|
||||
let basis_proj = DMatrix::from_columns(
|
||||
eig_pairs.filter_map(
|
||||
|(λ, v)| (λ.abs() < THRESHOLD).then_some(v)
|
||||
).collect::<Vec<_>>().as_slice()
|
||||
);
|
||||
|
||||
/* DEBUG */
|
||||
@ -125,30 +127,45 @@ impl ConfigSubspace {
|
||||
format!("Eigenvalues used to find kernel:{}", eig.eigenvalues)
|
||||
));
|
||||
|
||||
// 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;
|
||||
ConfigSubspace {
|
||||
assembly_dim: assembly_dim,
|
||||
basis: basis.collect()
|
||||
basis_std: basis_std.column_iter().map(
|
||||
|v| Into::<DMatrix<f64>>::into(
|
||||
v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim))
|
||||
)
|
||||
).collect(),
|
||||
basis_proj: basis_proj.column_iter().map(
|
||||
|v| Into::<DMatrix<f64>>::into(
|
||||
v.reshape_generic(Dyn(UNIFORM_DIM), Dyn(assembly_dim))
|
||||
)
|
||||
).collect()
|
||||
}
|
||||
}
|
||||
|
||||
pub fn dim(&self) -> usize {
|
||||
self.basis.len()
|
||||
self.basis_std.len()
|
||||
}
|
||||
|
||||
pub fn assembly_dim(&self) -> usize {
|
||||
self.assembly_dim
|
||||
}
|
||||
|
||||
// find the projection onto this subspace, with respect to the Euclidean
|
||||
// inner product, of the motion where the element with the given column
|
||||
// index has velocity `v`
|
||||
// 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
|
||||
// projection inner product
|
||||
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)
|
||||
} else {
|
||||
self.basis.iter().map(
|
||||
|b| b.column(column_index).dot(&v) * b
|
||||
self.basis_proj.iter().zip(self.basis_std.iter()).map(
|
||||
|(b_proj, b_std)| b_proj.column(column_index).dot(&v) * b_std
|
||||
).sum()
|
||||
}
|
||||
}
|
||||
@ -215,6 +232,37 @@ fn basis_matrix(index: (usize, usize), nrows: usize, ncols: usize) -> DMatrix<f6
|
||||
result
|
||||
}
|
||||
|
||||
// given a normalized vector `v` representing an element, build a basis for the
|
||||
// element's linear configuration space consisting of:
|
||||
// - the unit translation motions of the element
|
||||
// - the unit shrinking motion of the element, if it's a sphere
|
||||
// - one or two vectors whose coefficients vanish on the tangent space of the
|
||||
// normalization variety
|
||||
pub fn local_unif_to_std(v: DVectorView<f64>) -> DMatrix<f64> {
|
||||
const ELEMENT_DIM: usize = 5;
|
||||
const UNIFORM_DIM: usize = 4;
|
||||
let curv = 2.0*v[3];
|
||||
if v.dot(&(&*Q * v)) < 0.5 {
|
||||
// `v` represents a point. the normalization condition says that the
|
||||
// curvature component of `v` is 1/2
|
||||
DMatrix::from_column_slice(ELEMENT_DIM, UNIFORM_DIM, &[
|
||||
curv, 0.0, 0.0, 0.0, v[0],
|
||||
0.0, curv, 0.0, 0.0, v[1],
|
||||
0.0, 0.0, curv, 0.0, v[2],
|
||||
0.0, 0.0, 0.0, 0.0, 1.0
|
||||
])
|
||||
} else {
|
||||
// `v` represents a sphere. the normalization condition says that the
|
||||
// Lorentz product of `v` with itself is 1
|
||||
DMatrix::from_column_slice(ELEMENT_DIM, UNIFORM_DIM, &[
|
||||
curv, 0.0, 0.0, 0.0, v[0],
|
||||
0.0, curv, 0.0, 0.0, v[1],
|
||||
0.0, 0.0, curv, 0.0, v[2],
|
||||
curv*v[0], curv*v[1], curv*v[2], curv*v[3], curv*v[4] + 1.0
|
||||
])
|
||||
}
|
||||
}
|
||||
|
||||
// use backtracking line search to find a better configuration
|
||||
fn seek_better_config(
|
||||
gram: &PartialMatrix,
|
||||
@ -344,7 +392,19 @@ pub fn realize_gram(
|
||||
}
|
||||
let success = state.loss < tol;
|
||||
let tangent = if success {
|
||||
ConfigSubspace::symmetric_kernel(hess, assembly_dim)
|
||||
// express the uniform basis in the standard basis
|
||||
const UNIFORM_DIM: usize = 4;
|
||||
let total_dim_unif = UNIFORM_DIM * assembly_dim;
|
||||
let mut unif_to_std = DMatrix::<f64>::zeros(total_dim, total_dim_unif);
|
||||
for n in 0..assembly_dim {
|
||||
let block_start = (element_dim * n, UNIFORM_DIM * n);
|
||||
unif_to_std
|
||||
.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
|
||||
ConfigSubspace::symmetric_kernel(hess, unif_to_std, assembly_dim)
|
||||
} else {
|
||||
ConfigSubspace::zero(assembly_dim)
|
||||
};
|
||||
@ -424,6 +484,9 @@ 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]
|
||||
@ -486,10 +549,8 @@ mod tests {
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn tangent_test() {
|
||||
fn tangent_test_three_spheres() {
|
||||
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 {
|
||||
@ -514,32 +575,258 @@ mod tests {
|
||||
assert_eq!(history.scaled_loss.len(), 1);
|
||||
|
||||
// confirm that the tangent space has dimension five or less
|
||||
let ConfigSubspace(ref tangent_basis) = tangent;
|
||||
assert_eq!(tangent_basis.len(), 5);
|
||||
assert_eq!(tangent.basis_std.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
|
||||
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
|
||||
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 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(
|
||||
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,
|
||||
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(
|
||||
|(k, v)| tangent.proj(&v, k)
|
||||
).sum();
|
||||
assert!((motion - motion_proj).norm_squared() < tol_sq);
|
||||
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| {
|
||||
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() {
|
||||
// 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;
|
||||
const SCALED_TOL: f64 = 1.0e-12;
|
||||
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]
|
||||
|
||||
Loading…
Reference in New Issue
Block a user