From 01261aed91d4b4dae8f898b904a29720ddb5cf9a Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Tue, 21 Jan 2025 18:04:21 -0800 Subject: [PATCH 1/2] Write kaleidocycle nudge in uniform coordinates In the previous commit, the nudge was written in standard coordinates, so the example shouldn't have worked. However, by sheer dumb luck, the standard and uniform coordinates match for this particular nudge. In fact, the expression used before to get the standard coordinates may even produce equivalent floating point values as the expression used here to get the uniform coordinates. That would explain why the example prints exactly the same output in this commit. --- app-proto/examples/kaleidocycle.rs | 23 +++++++---------------- 1 file changed, 7 insertions(+), 16 deletions(-) diff --git a/app-proto/examples/kaleidocycle.rs b/app-proto/examples/kaleidocycle.rs index effb226..3448b87 100644 --- a/app-proto/examples/kaleidocycle.rs +++ b/app-proto/examples/kaleidocycle.rs @@ -1,4 +1,4 @@ -use nalgebra::DMatrix; +use nalgebra::{DMatrix, DVector}; use std::{array, f64::consts::PI}; use dyna3::engine::{Q, point, realize_gram, PartialMatrix}; @@ -59,22 +59,13 @@ fn main() { 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, 0.0]); + let down = -&up; let twist_motion: DMatrix<_> = (0..N_POINTS).step_by(4).flat_map( - |n| { - let up_field = { - DMatrix::from_column_slice(5, 5, &[ - 0.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 1.0, - 0.0, 0.0, 2.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 0.0 - ]) - }; - [ - tangent.proj(&(&up_field * config.column(n)).as_view(), n), - tangent.proj(&(-&up_field * config.column(n+1)).as_view(), n+1) - ] - } + |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); From 21cefa9f8a2f297b6cf36fbc31ae5249cca3a88f Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Tue, 21 Jan 2025 18:52:01 -0800 Subject: [PATCH 2/2] Do `symmetric_kernel` in projection coordinates Instead of finding the kernel in the standard coordinates and then expressing it in the projection coordinates, work in the projection coordinates from the beginning by applying a change of basis to the input matrix. --- app-proto/examples/kaleidocycle.rs | 2 +- app-proto/src/display.rs | 2 +- app-proto/src/engine.rs | 46 +++++++++++++++--------------- 3 files changed, 25 insertions(+), 25 deletions(-) diff --git a/app-proto/examples/kaleidocycle.rs b/app-proto/examples/kaleidocycle.rs index 3448b87..88116d3 100644 --- a/app-proto/examples/kaleidocycle.rs +++ b/app-proto/examples/kaleidocycle.rs @@ -59,7 +59,7 @@ fn main() { 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, 0.0]); + 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| [ diff --git a/app-proto/src/display.rs b/app-proto/src/display.rs index bbfe059..4e0c7e4 100644 --- a/app-proto/src/display.rs +++ b/app-proto/src/display.rs @@ -346,7 +346,7 @@ pub fn Display() -> View { Vector3::zeros() }; time_step * DVector::from_column_slice( - &[u[0], u[1], u[2], SHRINKING_SPEED * shrink, 0.0] + &[u[0], u[1], u[2], SHRINKING_SPEED * shrink] ) }; assembly_for_raf.deform( diff --git a/app-proto/src/engine.rs b/app-proto/src/engine.rs index 9db6db0..0f1fd36 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -107,13 +107,14 @@ impl ConfigSubspace { // 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, proj_to_std: DMatrix, assembly_dim: usize) -> ConfigSubspace { - // find a basis for the kernel, expressed in the standard coordinates - const ELEMENT_DIM: usize = 5; - const THRESHOLD: f64 = 1.0e-4; - let eig = SymmetricEigen::new(a); + // 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_std = DMatrix::from_columns( + let basis_proj = DMatrix::from_columns( eig_pairs.filter_map( |(λ, v)| (λ.abs() < THRESHOLD).then_some(v) ).collect::>().as_slice() @@ -126,29 +127,27 @@ impl ConfigSubspace { format!("Eigenvalues used to find kernel:{}", eig.eigenvalues) )); - // express the basis in the projection coordinates - let basis_proj = proj_to_std.clone().qr().solve(&basis_std).unwrap(); - - // orthonormalize the basis with respect to the projection inner product - let basis_proj_orth = basis_proj.qr().q(); - let basis_std_orth = proj_to_std * &basis_proj_orth; + // 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_orth) + format!("Basis in projection coordinates:{}", basis_proj) )); + const ELEMENT_DIM: usize = 5; + const UNIFORM_DIM: usize = 4; ConfigSubspace { assembly_dim: assembly_dim, - basis_std: basis_std_orth.column_iter().map( + basis_std: basis_std.column_iter().map( |v| Into::>::into( v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim)) ) ).collect(), - basis_proj: basis_proj_orth.column_iter().map( + basis_proj: basis_proj.column_iter().map( |v| Into::>::into( - v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim)) + v.reshape_generic(Dyn(UNIFORM_DIM), Dyn(assembly_dim)) ) ).collect() } @@ -247,26 +246,25 @@ 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, ELEMENT_DIM, &[ + 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], - v[0], v[1], v[2], v[3], v[4], 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, ELEMENT_DIM, &[ + 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, - v[0], v[1], v[2], v[3], v[4] + curv*v[0], curv*v[1], curv*v[2], curv*v[3], curv*v[4] + 1.0 ]) } } @@ -401,11 +399,13 @@ pub fn realize_gram( let success = state.loss < tol; let tangent = if success { // express the uniform basis in the standard basis - let mut unif_to_std = DMatrix::::zeros(total_dim, total_dim); + 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; + let block_start = (element_dim * n, UNIFORM_DIM * n); unif_to_std - .view_mut((block_start, block_start), (element_dim, element_dim)) + .view_mut(block_start, (element_dim, UNIFORM_DIM)) .copy_from(&local_unif_to_std(state.config.column(n))); }