diff --git a/app-proto/run-examples b/app-proto/run-examples index b5c3de5..6a5e3ae 100755 --- a/app-proto/run-examples +++ b/app-proto/run-examples @@ -3,5 +3,6 @@ # https://jonalmeida.com/posts/2015/01/23/print-cargo/ # +cargo test -- --nocapture engine::tests::irisawa_hexlet_test cargo test -- --nocapture engine::tests::three_spheres_example cargo test -- --nocapture engine::tests::point_on_sphere_example diff --git a/app-proto/src/engine.rs b/app-proto/src/engine.rs index f5ad19a..6f2952c 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -37,7 +37,7 @@ pub fn sphere_with_offset(dir_x: f64, dir_y: f64, dir_z: f64, off: f64, curv: f6 struct MatrixEntry { index: (usize, usize), - val: f64 + value: f64 } struct PartialMatrix(Vec); @@ -56,7 +56,7 @@ impl PartialMatrix { let mut result = DMatrix::::zeros(rhs.nrows(), rhs.ncols()); let PartialMatrix(entries) = self; for ent in entries { - result[ent.index] = ent.val - rhs[ent.index]; + result[ent.index] = ent.value - rhs[ent.index]; } result } @@ -141,7 +141,7 @@ fn realize_gram( let total_dim = element_dim * assembly_dim; // scale the tolerance - let scale_adjustment = ((guess.ncols() - frozen.len()) as f64).sqrt(); + let scale_adjustment = (gram.0.len() as f64).sqrt(); let tol = scale_adjustment * scaled_tol; // convert the frozen indices to stacked format @@ -153,8 +153,8 @@ fn realize_gram( let mut state = SearchState::from_config(gram, guess); for _ in 0..max_descent_steps { // stop if the loss is tolerably low - println!("loss: {}", state.loss); - /*println!("projected error: {}", state.err_proj);*/ + println!("scaled loss: {}", state.loss / scale_adjustment); + /* println!("projected error: {}", state.err_proj); */ if state.loss < tol { break; } // find the negative gradient of the loss function @@ -182,6 +182,7 @@ fn realize_gram( // regularize the Hessian let min_eigval = hess.symmetric_eigenvalues().min(); + /* println!("lowest eigenvalue: {}", min_eigval); */ if min_eigval <= 0.0 { hess -= reg_scale * min_eigval * DMatrix::identity(total_dim, total_dim); } @@ -198,6 +199,12 @@ fn realize_gram( } // compute the Newton step + /* + we need to either handle or eliminate the case where the minimum + eigenvalue of the Hessian is zero, so the regularized Hessian is + singular. right now, this causes the Cholesky decomposition to return + `None`, leading to a panic when we unrap + */ let base_step_stacked = hess.cholesky().unwrap().solve(&neg_grad_stacked); let base_step = base_step_stacked.reshape_generic(Dyn(element_dim), Dyn(assembly_dim)); @@ -217,17 +224,17 @@ fn realize_gram( #[cfg(test)] mod tests { - use std::f64; + use std::{array, f64::consts::PI}; use super::*; #[test] fn sub_proj_test() { let target = PartialMatrix(vec![ - MatrixEntry { index: (0, 0), val: 19.0 }, - MatrixEntry { index: (0, 2), val: 39.0 }, - MatrixEntry { index: (1, 1), val: 59.0 }, - MatrixEntry { index: (1, 2), val: 69.0 } + MatrixEntry { index: (0, 0), value: 19.0 }, + MatrixEntry { index: (0, 2), value: 39.0 }, + MatrixEntry { index: (1, 1), value: 59.0 }, + MatrixEntry { index: (1, 2), value: 69.0 } ]); let attempt = DMatrix::::from_row_slice(2, 3, &[ 1.0, 2.0, 3.0, @@ -248,7 +255,7 @@ mod tests { for k in 0..3 { entries.push(MatrixEntry { index: (j, k), - val: if j == k { 1.0 } else { -1.0 } + value: if j == k { 1.0 } else { -1.0 } }); } } @@ -266,6 +273,88 @@ mod tests { assert!(state.loss.abs() < f64::EPSILON); } + // this problem is from a sangaku by Irisawa Shintarō Hiroatsu. the article + // below includes a nice translation of the problem statement, which was + // recorded in Uchida Itsumi's book _Kokon sankan_ (_Mathematics, Past and + // Present_) + // + // "Japan's 'Wasan' Mathematical Tradition", by Abe Haruki + // https://www.nippon.com/en/japan-topics/c12801/ + // + #[test] + fn irisawa_hexlet_test() { + let gram = PartialMatrix({ + let mut entries = Vec::::new(); + for s in 0..9 { + // each sphere is represented by a spacelike vector + entries.push(MatrixEntry { index: (s, s), value: 1.0 }); + + // the circumscribing sphere is tangent to all of the other + // spheres, with matching orientation + if s > 0 { + entries.push(MatrixEntry { index: (0, s), value: 1.0 }); + entries.push(MatrixEntry { index: (s, 0), value: 1.0 }); + } + + if s > 2 { + // each chain sphere is tangent to the "sun" and "moon" + // spheres, with opposing orientation + for n in 1..3 { + entries.push(MatrixEntry { index: (s, n), value: -1.0 }); + entries.push(MatrixEntry { index: (n, s), value: -1.0 }); + } + + // each chain sphere is tangent to the next chain sphere, + // with opposing orientation + let s_next = 3 + (s-2) % 6; + entries.push(MatrixEntry { index: (s, s_next), value: -1.0 }); + entries.push(MatrixEntry { index: (s_next, s), value: -1.0 }); + } + } + entries + }); + let guess = DMatrix::from_columns( + [ + sphere(0.0, 0.0, 0.0, 15.0), + sphere(0.0, 0.0, -9.0, 5.0), + sphere(0.0, 0.0, 11.0, 3.0) + ].into_iter().chain( + (1..=6).map( + |k| { + let ang = (k as f64) * PI/3.0; + sphere(9.0 * ang.cos(), 9.0 * ang.sin(), 0.0, 2.5) + } + ) + ).collect::>().as_slice() + ); + let frozen: [(usize, usize); 4] = array::from_fn(|k| (3, k)); + const SCALED_TOL: f64 = 1.0e-12; + let (config, success) = realize_gram( + &gram, guess, &frozen, + SCALED_TOL, 0.5, 0.9, 1.1, 200, 110 + ); + print!("\nCompleted Gram matrix:{}", config.tr_mul(&*Q) * &config); + let final_state = SearchState::from_config(&gram, config); + if success { + println!("Target accuracy achieved!"); + } else { + println!("Failed to reach target accuracy"); + } + println!("Loss: {}", final_state.loss); + if success { + println!("\nChain diameters:"); + println!(" {} sun (given)", 1.0 / final_state.config[(3, 3)]); + for k in 4..9 { + println!(" {} sun", 1.0 / final_state.config[(3, k)]); + } + } + let entry_tol = SCALED_TOL.sqrt(); + let solution_diams = [30.0, 10.0, 6.0, 5.0, 15.0, 10.0, 3.75, 2.5, 2.0 + 8.0/11.0]; + for (k, diam) in solution_diams.into_iter().enumerate() { + assert!((final_state.config[(3, k)] - 1.0 / diam).abs() < entry_tol); + } + } + // --- process inspection examples --- // these tests are meant for human inspection, not automated use. run them @@ -281,7 +370,7 @@ mod tests { for k in 0..3 { entries.push(MatrixEntry { index: (j, k), - val: if j == k { 1.0 } else { -1.0 } + value: if j == k { 1.0 } else { -1.0 } }); } } @@ -318,7 +407,7 @@ mod tests { for k in 0..2 { entries.push(MatrixEntry { index: (j, k), - val: if (j, k) == (1, 1) { 1.0 } else { 0.0 } + value: if (j, k) == (1, 1) { 1.0 } else { 0.0 } }); } } @@ -328,9 +417,10 @@ mod tests { point(0.0, 0.0, 2.0), sphere(0.0, 0.0, 0.0, 1.0) ]); + let frozen = [(3, 0)]; println!(); let (config, success) = realize_gram( - &gram, guess, &[(3, 0)], + &gram, guess, &frozen, 1.0e-12, 0.5, 0.9, 1.1, 200, 110 ); print!("\nCompleted Gram matrix:{}", config.tr_mul(&*Q) * &config);