diff --git a/app-proto/examples/kaleidocycle.rs b/app-proto/examples/kaleidocycle.rs new file mode 100644 index 0000000..88116d3 --- /dev/null +++ b/app-proto/examples/kaleidocycle.rs @@ -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::>(); + 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); +} \ No newline at end of file diff --git a/app-proto/run-examples b/app-proto/run-examples index bc7e933..52173b0 100755 --- a/app-proto/run-examples +++ b/app-proto/run-examples @@ -9,3 +9,4 @@ cargo run --example irisawa-hexlet cargo run --example three-spheres cargo run --example point-on-sphere +cargo run --example kaleidocycle \ No newline at end of file diff --git a/app-proto/src/assembly.rs b/app-proto/src/assembly.rs index 278a8f9..9b0e065 100644 --- a/app-proto/src/assembly.rs +++ b/app-proto/src/assembly.rs @@ -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 } +// 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; } } diff --git a/app-proto/src/display.rs b/app-proto/src/display.rs index d8ee23c..4e0c7e4 100644 --- a/app-proto/src/display.rs +++ b/app-proto/src/display.rs @@ -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 = 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 { diff --git a/app-proto/src/engine.rs b/app-proto/src/engine.rs index 6a0202e..0f1fd36 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -90,32 +90,34 @@ impl PartialMatrix { #[derive(Clone)] pub struct ConfigSubspace { assembly_dim: usize, - basis: Vec> + basis_std: Vec>, + basis_proj: Vec> } 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, 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, proj_to_std: DMatrix, 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::>::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::>().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::>::into( + v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim)) + ) + ).collect(), + basis_proj: basis_proj.column_iter().map( + |v| Into::>::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, column_index: usize) -> DMatrix { 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) -> DMatrix { + 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::::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) };