diff --git a/app-proto/src/engine.rs b/app-proto/src/engine.rs index 57243b8..0f1fd36 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -490,9 +490,6 @@ pub mod irisawa { #[cfg(test)] mod tests { - use nalgebra::Vector3; - use std::{array, f64::consts::{FRAC_1_SQRT_2, PI}, iter}; - use super::{*, irisawa::realize_irisawa_hexlet}; #[test] @@ -555,8 +552,10 @@ mod tests { } #[test] - fn tangent_test_three_spheres() { + fn tangent_test() { const SCALED_TOL: f64 = 1.0e-12; + const ELEMENT_DIM: usize = 5; + const ASSEMBLY_DIM: usize = 3; let gram = { let mut gram_to_be = PartialMatrix::new(); for j in 0..3 { @@ -581,258 +580,32 @@ mod tests { assert_eq!(history.scaled_loss.len(), 1); // confirm that the tangent space has dimension five or less - assert_eq!(tangent.basis_std.len(), 5); + let ConfigSubspace(ref tangent_basis) = tangent; + assert_eq!(tangent_basis.len(), 5); // confirm that the tangent space contains all the motions we expect it // to. since we've already bounded the dimension of the tangent space, // this confirms that the tangent space is what we expect it to be - const UNIFORM_DIM: usize = 4; - let element_dim = guess.nrows(); - let assembly_dim = guess.ncols(); - let tangent_motions_unif = vec![ - basis_matrix((0, 1), UNIFORM_DIM, assembly_dim), - basis_matrix((1, 1), UNIFORM_DIM, assembly_dim), - basis_matrix((0, 2), UNIFORM_DIM, assembly_dim), - basis_matrix((1, 2), UNIFORM_DIM, assembly_dim), - DMatrix::::from_column_slice(UNIFORM_DIM, assembly_dim, &[ - 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, -0.5, -0.5, - 0.0, 0.0, -0.5, 0.5 - ]) - ]; - let tangent_motions_std = vec![ - basis_matrix((0, 1), element_dim, assembly_dim), - basis_matrix((1, 1), element_dim, assembly_dim), - basis_matrix((0, 2), element_dim, assembly_dim), - basis_matrix((1, 2), element_dim, assembly_dim), - DMatrix::::from_column_slice(element_dim, assembly_dim, &[ - 0.0, 0.0, 0.0, 0.00, 0.0, + let tangent_motions = vec![ + basis_matrix((0, 1), ELEMENT_DIM, ASSEMBLY_DIM), + basis_matrix((1, 1), ELEMENT_DIM, ASSEMBLY_DIM), + basis_matrix((0, 2), ELEMENT_DIM, ASSEMBLY_DIM), + basis_matrix((1, 2), ELEMENT_DIM, ASSEMBLY_DIM), + DMatrix::::from_column_slice(ELEMENT_DIM, 3, &[ + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, -0.25, -1.0, - 0.0, 0.0, -1.0, 0.25, 1.0 + 0.0, 0.0, -1.0, 0.25, 1.0 ]) ]; - let tol_sq = ((element_dim * assembly_dim) as f64) * SCALED_TOL * SCALED_TOL; - for (motion_unif, motion_std) in tangent_motions_unif.into_iter().zip(tangent_motions_std) { - let motion_proj: DMatrix<_> = motion_unif.column_iter().enumerate().map( + let tol_sq = ((ELEMENT_DIM * ASSEMBLY_DIM) as f64) * SCALED_TOL * SCALED_TOL; + for motion in tangent_motions { + let motion_proj: DMatrix<_> = motion.column_iter().enumerate().map( |(k, v)| tangent.proj(&v, k) ).sum(); - assert!((motion_std - motion_proj).norm_squared() < tol_sq); + assert!((motion - motion_proj).norm_squared() < tol_sq); } } - fn translation_motion_unif(vel: &Vector3, assembly_dim: usize) -> Vec> { - let mut elt_motion = DVector::zeros(4); - elt_motion.fixed_rows_mut::<3>(0).copy_from(vel); - iter::repeat(elt_motion).take(assembly_dim).collect() - } - - fn rotation_motion_unif(ang_vel: &Vector3, points: Vec>) -> Vec> { - points.into_iter().map( - |pt| { - let vel = ang_vel.cross(&pt.fixed_rows::<3>(0)); - let mut elt_motion = DVector::zeros(4); - elt_motion.fixed_rows_mut::<3>(0).copy_from(&vel); - elt_motion - } - ).collect() - } - - #[test] - fn tangent_test_kaleidocycle() { - const SCALED_TOL: f64 = 1.0e-12; - // set up a kaleidocycle, made of points with fixed distances between - // them, and find its tangent space - const N_POINTS: usize = 12; - const N_HINGES: usize = 6; - 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 = { - 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.clone(), &frozen, - SCALED_TOL, 0.5, 0.9, 1.1, 200, 110 - ); - assert_eq!(config, guess); - assert_eq!(success, true); - assert_eq!(history.scaled_loss.len(), 1); - - // list some motions that should form a basis for the tangent space of - // the solution variety - let element_dim = guess.nrows(); - let assembly_dim = guess.ncols(); - let tangent_motions_unif = vec![ - // the translations along the coordinate axes - translation_motion_unif(&Vector3::new(1.0, 0.0, 0.0), assembly_dim), - translation_motion_unif(&Vector3::new(0.0, 1.0, 0.0), assembly_dim), - translation_motion_unif(&Vector3::new(0.0, 0.0, 1.0), assembly_dim), - - // the rotations about the coordinate axes - rotation_motion_unif(&Vector3::new(1.0, 0.0, 0.0), guess.column_iter().collect()), - rotation_motion_unif(&Vector3::new(0.0, 1.0, 0.0), guess.column_iter().collect()), - rotation_motion_unif(&Vector3::new(0.0, 0.0, 1.0), guess.column_iter().collect()), - - // the twist motion. more precisely: a motion that keeps the center - // of mass stationary and preserves the distances between the - // vertices to first order. this has to be the twist as long as: - // - twisting is the kaleidocycle's only internal degree of - // freedom - // - every first-order motion of the kaleidocycle comes from an - // actual motion - (0..N_HINGES).step_by(2).flat_map( - |n| { - let ang_vert = ((n + 1) as f64) * PI/3.0; - let vel_vert_x = 4.0 * ang_vert.cos(); - let vel_vert_y = 4.0 * ang_vert.sin(); - [ - DVector::from_column_slice(&[0.0, 0.0, 5.0, 0.0]), - DVector::from_column_slice(&[0.0, 0.0, 1.0, 0.0]), - DVector::from_column_slice(&[-vel_vert_x, -vel_vert_y, -3.0, 0.0]), - DVector::from_column_slice(&[vel_vert_x, vel_vert_y, -3.0, 0.0]) - ] - } - ).collect::>() - ]; - let tangent_motions_std = tangent_motions_unif.iter().map( - |motion| DMatrix::from_columns( - &guess.column_iter().zip(motion).map( - |(v, elt_motion)| local_unif_to_std(v) * elt_motion - ).collect::>() - ) - ).collect::>(); - - // confirm that the dimension of the tangent space is no greater than - // expected - assert_eq!(tangent.basis_std.len(), tangent_motions_unif.len()); - - // confirm that the tangent space contains all the motions we expect it - // to. since we've already bounded the dimension of the tangent space, - // this confirms that the tangent space is what we expect it to be - let tol_sq = ((element_dim * assembly_dim) as f64) * SCALED_TOL * SCALED_TOL; - for (motion_unif, motion_std) in tangent_motions_unif.into_iter().zip(tangent_motions_std) { - let motion_proj: DMatrix<_> = motion_unif.into_iter().enumerate().map( - |(k, v)| tangent.proj(&v.as_view(), k) - ).sum(); - assert!((motion_std - motion_proj).norm_squared() < tol_sq); - } - } - - fn translation(dis: Vector3) -> DMatrix { - const ELEMENT_DIM: usize = 5; - DMatrix::from_column_slice(ELEMENT_DIM, ELEMENT_DIM, &[ - 1.0, 0.0, 0.0, 0.0, dis[0], - 0.0, 1.0, 0.0, 0.0, dis[1], - 0.0, 0.0, 1.0, 0.0, dis[2], - 2.0*dis[0], 2.0*dis[1], 2.0*dis[2], 1.0, dis.norm_squared(), - 0.0, 0.0, 0.0, 0.0, 1.0 - ]) - } - - // confirm that projection onto a configuration subspace is equivariant with - // respect to Euclidean motions - #[test] - fn proj_equivar_test() { - // find a pair of spheres that meet at 120° - const SCALED_TOL: f64 = 1.0e-12; - let gram = { - let mut gram_to_be = PartialMatrix::new(); - gram_to_be.push_sym(0, 0, 1.0); - gram_to_be.push_sym(1, 1, 1.0); - gram_to_be.push_sym(0, 1, 0.5); - gram_to_be - }; - let guess_orig = DMatrix::from_columns(&[ - sphere(0.0, 0.0, 0.5, 1.0), - sphere(0.0, 0.0, -0.5, 1.0) - ]); - let (config_orig, tangent_orig, success_orig, history_orig) = realize_gram( - &gram, guess_orig.clone(), &[], - SCALED_TOL, 0.5, 0.9, 1.1, 200, 110 - ); - assert_eq!(config_orig, guess_orig); - assert_eq!(success_orig, true); - assert_eq!(history_orig.scaled_loss.len(), 1); - - // find another pair of spheres that meet at 120°. we'll think of this - // solution as a transformed version of the original one - let guess_tfm = { - let a = 0.5 * FRAC_1_SQRT_2; - DMatrix::from_columns(&[ - sphere(a, 0.0, 7.0 + a, 1.0), - sphere(-a, 0.0, 7.0 - a, 1.0) - ]) - }; - let (config_tfm, tangent_tfm, success_tfm, history_tfm) = realize_gram( - &gram, guess_tfm.clone(), &[], - SCALED_TOL, 0.5, 0.9, 1.1, 200, 110 - ); - assert_eq!(config_tfm, guess_tfm); - assert_eq!(success_tfm, true); - assert_eq!(history_tfm.scaled_loss.len(), 1); - - // project a nudge to the tangent space of the solution variety at the - // original solution - let motion_orig = DVector::from_column_slice(&[0.0, 0.0, 1.0, 0.0]); - let motion_orig_proj = tangent_orig.proj(&motion_orig.as_view(), 0); - - // project the equivalent nudge to the tangent space of the solution - // variety at the transformed solution - let motion_tfm = DVector::from_column_slice(&[FRAC_1_SQRT_2, 0.0, FRAC_1_SQRT_2, 0.0]); - let motion_tfm_proj = tangent_tfm.proj(&motion_tfm.as_view(), 0); - - // take the transformation that sends the original solution to the - // transformed solution and apply it to the motion that the original - // solution makes in response to the nudge - const ELEMENT_DIM: usize = 5; - let rot = DMatrix::from_column_slice(ELEMENT_DIM, ELEMENT_DIM, &[ - FRAC_1_SQRT_2, 0.0, -FRAC_1_SQRT_2, 0.0, 0.0, - 0.0, 1.0, 0.0, 0.0, 0.0, - FRAC_1_SQRT_2, 0.0, FRAC_1_SQRT_2, 0.0, 0.0, - 0.0, 0.0, 0.0, 1.0, 0.0, - 0.0, 0.0, 0.0, 0.0, 1.0 - ]); - let transl = translation(Vector3::new(0.0, 0.0, 7.0)); - let motion_proj_tfm = transl * rot * motion_orig_proj; - - // confirm that the projection of the nudge is equivariant. we loosen - // the comparison tolerance because the transformation seems to - // introduce some numerical error - const SCALED_TOL_TFM: f64 = 1.0e-9; - let tol_sq = ((guess_orig.nrows() * guess_orig.ncols()) as f64) * SCALED_TOL_TFM * SCALED_TOL_TFM; - assert!((motion_proj_tfm - motion_tfm_proj).norm_squared() < tol_sq); - } - // at the frozen indices, the optimization steps should have exact zeros, // and the realized configuration should match the initial guess #[test]