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

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<MatrixEntry>);
@@ -56,7 +56,7 @@ impl PartialMatrix {
         let mut result = DMatrix::<f64>::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!("scaled loss: {}", state.loss / scale_adjustment);
-        /*println!("projected error: {}", state.err_proj);*/
+        /* 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, 0), value: 19.0 },
-            MatrixEntry { index: (0, 2), val: 39.0 },
+            MatrixEntry { index: (0, 2), value: 39.0 },
-            MatrixEntry { index: (1, 1), val: 59.0 },
+            MatrixEntry { index: (1, 1), value: 59.0 },
-            MatrixEntry { index: (1, 2), val: 69.0 }
+            MatrixEntry { index: (1, 2), value: 69.0 }
         ]);
         let attempt = DMatrix::<f64>::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::<MatrixEntry>::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::<Vec<_>>().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);