forked from glen/dyna3
Vectornaut
817a446fad
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: glen/dyna3#34 Co-authored-by: Vectornaut <vectornaut@nobody@nowhere.net> Co-committed-by: Vectornaut <vectornaut@nobody@nowhere.net>
72 lines
No EOL
2.5 KiB
Rust
72 lines
No EOL
2.5 KiB
Rust
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);
|
|
} |