glen/dyna3
1
0
Fork 0
Code Issues 18 Pull Requests 1 Actions Projects Releases 2 Wiki Activity

Switch to Euclidean-invariant projection onto tangent space of solution variety #34

Open
Vectornaut wants to merge 4 commits from uniform-projection into main
pull from: uniform-projection
merge into: glen:main
Conversation 1 Commits 4 Files Changed 5 +193 -36
5 changed files with 193 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

72
app-proto/examples/kaleidocycle.rs Normal file
View 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);
}

1
app-proto/run-examples
View File

@ -9,3 +9,4 @@
cargo run --example irisawa-hexlet
cargo run --example three-spheres
cargo run --example point-on-sphere
cargo run --example kaleidocycle

12
app-proto/src/assembly.rs
View File

@ -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;
}
}

38
app-proto/src/display.rs
View File

@ -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 {

106
app-proto/src/engine.rs
View File

@ -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,51 @@ 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;
// print the projection basis in projection coordinates
#[cfg(all(target_family = "wasm", target_os = "unknown"))]
console::log_1(&JsValue::from(
format!("Basis in projection coordinates:{}", 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 +238,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 +398,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)
};