From a05a2e1d5498f1d2263d3622133cece4d748773d Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Mon, 20 Jan 2025 15:44:20 -0800 Subject: [PATCH] Implement the uniform inner product for points --- app-proto/examples/kaleidocycle.rs | 81 ++++++++++++++++++++++++++++++ app-proto/run-examples | 1 + app-proto/src/engine.rs | 40 +++++++++++---- 3 files changed, 112 insertions(+), 10 deletions(-) create mode 100644 app-proto/examples/kaleidocycle.rs diff --git a/app-proto/examples/kaleidocycle.rs b/app-proto/examples/kaleidocycle.rs new file mode 100644 index 0000000..effb226 --- /dev/null +++ b/app-proto/examples/kaleidocycle.rs @@ -0,0 +1,81 @@ +use nalgebra::DMatrix; +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 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) + ] + } + ).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/engine.rs b/app-proto/src/engine.rs index c6b6150..9db6db0 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -132,6 +132,9 @@ impl ConfigSubspace { // 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; + + // 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) )); @@ -236,19 +239,36 @@ fn basis_matrix(index: (usize, usize), nrows: usize, ncols: usize) -> DMatrix) -> DMatrix { const ELEMENT_DIM: usize = 5; let curv = 2.0*v[3]; - DMatrix::from_column_slice(ELEMENT_DIM, ELEMENT_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] - ]) + 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, &[ + 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, &[ + 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] + ]) + } } // use backtracking line search to find a better configuration