diff --git a/app-proto/Cargo.toml b/app-proto/Cargo.toml index 04ea271..920469a 100644 --- a/app-proto/Cargo.toml +++ b/app-proto/Cargo.toml @@ -10,7 +10,6 @@ default = ["console_error_panic_hook"] [dependencies] itertools = "0.13.0" js-sys = "0.3.70" -lazy_static = "1.5.0" nalgebra = "0.33.0" rustc-hash = "2.0.0" slab = "0.4.9" diff --git a/app-proto/run-examples b/app-proto/run-examples deleted file mode 100755 index 6a5e3ae..0000000 --- a/app-proto/run-examples +++ /dev/null @@ -1,8 +0,0 @@ -# based on "Enabling print statements in Cargo tests", by Jon Almeida -# -# 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/add_remove.rs b/app-proto/src/add_remove.rs index 7066089..ab5db70 100644 --- a/app-proto/src/add_remove.rs +++ b/app-proto/src/add_remove.rs @@ -12,8 +12,7 @@ fn load_gen_assemb(assembly: &Assembly) { label: String::from("Castor"), color: [1.00_f32, 0.25_f32, 0.00_f32], rep: engine::sphere(0.5, 0.5, 0.0, 1.0), - constraints: BTreeSet::default(), - index: 0 + constraints: BTreeSet::default() } ); let _ = assembly.try_insert_element( @@ -22,8 +21,7 @@ fn load_gen_assemb(assembly: &Assembly) { label: String::from("Pollux"), color: [0.00_f32, 0.25_f32, 1.00_f32], rep: engine::sphere(-0.5, -0.5, 0.0, 1.0), - constraints: BTreeSet::default(), - index: 0 + constraints: BTreeSet::default() } ); let _ = assembly.try_insert_element( @@ -32,8 +30,7 @@ fn load_gen_assemb(assembly: &Assembly) { label: String::from("Ursa major"), color: [0.25_f32, 0.00_f32, 1.00_f32], rep: engine::sphere(-0.5, 0.5, 0.0, 0.75), - constraints: BTreeSet::default(), - index: 0 + constraints: BTreeSet::default() } ); let _ = assembly.try_insert_element( @@ -42,8 +39,7 @@ fn load_gen_assemb(assembly: &Assembly) { label: String::from("Ursa minor"), color: [0.25_f32, 1.00_f32, 0.00_f32], rep: engine::sphere(0.5, -0.5, 0.0, 0.5), - constraints: BTreeSet::default(), - index: 0 + constraints: BTreeSet::default() } ); let _ = assembly.try_insert_element( @@ -52,8 +48,7 @@ fn load_gen_assemb(assembly: &Assembly) { label: String::from("Deimos"), color: [0.75_f32, 0.75_f32, 0.00_f32], rep: engine::sphere(0.0, 0.15, 1.0, 0.25), - constraints: BTreeSet::default(), - index: 0 + constraints: BTreeSet::default() } ); let _ = assembly.try_insert_element( @@ -62,8 +57,17 @@ fn load_gen_assemb(assembly: &Assembly) { label: String::from("Phobos"), color: [0.00_f32, 0.75_f32, 0.50_f32], rep: engine::sphere(0.0, -0.15, -1.0, 0.25), - constraints: BTreeSet::default(), - index: 0 + constraints: BTreeSet::default() + } + ); + assembly.insert_constraint( + Constraint { + args: ( + assembly.elements_by_id.with_untracked(|elts_by_id| elts_by_id["gemini_a"]), + assembly.elements_by_id.with_untracked(|elts_by_id| elts_by_id["gemini_b"]) + ), + rep: 0.5, + active: create_signal(true) } ); } @@ -77,8 +81,7 @@ fn load_low_curv_assemb(assembly: &Assembly) { label: "Central".to_string(), color: [0.75_f32, 0.75_f32, 0.75_f32], rep: engine::sphere(0.0, 0.0, 0.0, 1.0), - constraints: BTreeSet::default(), - index: 0 + constraints: BTreeSet::default() } ); let _ = assembly.try_insert_element( @@ -87,8 +90,7 @@ fn load_low_curv_assemb(assembly: &Assembly) { label: "Assembly plane".to_string(), color: [0.75_f32, 0.75_f32, 0.75_f32], rep: engine::sphere_with_offset(0.0, 0.0, 1.0, 0.0, 0.0), - constraints: BTreeSet::default(), - index: 0 + constraints: BTreeSet::default() } ); let _ = assembly.try_insert_element( @@ -97,8 +99,7 @@ fn load_low_curv_assemb(assembly: &Assembly) { label: "Side 1".to_string(), color: [1.00_f32, 0.00_f32, 0.25_f32], rep: engine::sphere_with_offset(1.0, 0.0, 0.0, 1.0, 0.0), - constraints: BTreeSet::default(), - index: 0 + constraints: BTreeSet::default() } ); let _ = assembly.try_insert_element( @@ -107,8 +108,7 @@ fn load_low_curv_assemb(assembly: &Assembly) { label: "Side 2".to_string(), color: [0.25_f32, 1.00_f32, 0.00_f32], rep: engine::sphere_with_offset(-0.5, a, 0.0, 1.0, 0.0), - constraints: BTreeSet::default(), - index: 0 + constraints: BTreeSet::default() } ); let _ = assembly.try_insert_element( @@ -117,8 +117,7 @@ fn load_low_curv_assemb(assembly: &Assembly) { label: "Side 3".to_string(), color: [0.00_f32, 0.25_f32, 1.00_f32], rep: engine::sphere_with_offset(-0.5, -a, 0.0, 1.0, 0.0), - constraints: BTreeSet::default(), - index: 0 + constraints: BTreeSet::default() } ); let _ = assembly.try_insert_element( @@ -127,8 +126,7 @@ fn load_low_curv_assemb(assembly: &Assembly) { label: "Corner 1".to_string(), color: [0.75_f32, 0.75_f32, 0.75_f32], rep: engine::sphere(-4.0/3.0, 0.0, 0.0, 1.0/3.0), - constraints: BTreeSet::default(), - index: 0 + constraints: BTreeSet::default() } ); let _ = assembly.try_insert_element( @@ -137,8 +135,7 @@ fn load_low_curv_assemb(assembly: &Assembly) { label: "Corner 2".to_string(), color: [0.75_f32, 0.75_f32, 0.75_f32], rep: engine::sphere(2.0/3.0, -4.0/3.0 * a, 0.0, 1.0/3.0), - constraints: BTreeSet::default(), - index: 0 + constraints: BTreeSet::default() } ); let _ = assembly.try_insert_element( @@ -147,8 +144,7 @@ fn load_low_curv_assemb(assembly: &Assembly) { label: String::from("Corner 3"), color: [0.75_f32, 0.75_f32, 0.75_f32], rep: engine::sphere(2.0/3.0, 4.0/3.0 * a, 0.0, 1.0/3.0), - constraints: BTreeSet::default(), - index: 0 + constraints: BTreeSet::default() } ); } @@ -219,7 +215,6 @@ pub fn AddRemove() -> View { rep: 0.0, active: create_signal(true) }); - state.assembly.realize(); state.selection.update(|sel| sel.clear()); /* DEBUG */ diff --git a/app-proto/src/assembly.rs b/app-proto/src/assembly.rs index 228357e..e8dab79 100644 --- a/app-proto/src/assembly.rs +++ b/app-proto/src/assembly.rs @@ -1,11 +1,8 @@ -use nalgebra::{DMatrix, DVector}; +use nalgebra::DVector; use rustc_hash::FxHashMap; use slab::Slab; use std::collections::BTreeSet; use sycamore::prelude::*; -use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */ - -use crate::engine::{realize_gram, PartialMatrix}; #[derive(Clone, PartialEq)] pub struct Element { @@ -13,10 +10,7 @@ pub struct Element { pub label: String, pub color: [f32; 3], pub rep: DVector, - pub constraints: BTreeSet, - - // internal properties, not reflected in any view - pub index: usize + pub constraints: BTreeSet } #[derive(Clone)] @@ -46,8 +40,6 @@ impl Assembly { } } - // --- inserting elements and constraints --- - // insert an element into the assembly without checking whether we already // have an element with the same identifier. any element that does have the // same identifier will get kicked out of the `elements_by_id` index @@ -85,8 +77,7 @@ impl Assembly { label: format!("Sphere {}", id_num), color: [0.75_f32, 0.75_f32, 0.75_f32], rep: DVector::::from_column_slice(&[0.0, 0.0, 0.0, 0.5, -0.5]), - constraints: BTreeSet::default(), - index: 0 + constraints: BTreeSet::default() } ); } @@ -99,83 +90,4 @@ impl Assembly { elts[args.1].constraints.insert(key); }) } - - // --- realization --- - - pub fn realize(&self) { - // index the elements - self.elements.update_silent(|elts| { - for (index, (_, elt)) in elts.into_iter().enumerate() { - elt.index = index; - } - }); - - // set up the Gram matrix and the initial configuration matrix - let (gram, guess) = self.elements.with_untracked(|elts| { - // set up the off-diagonal part of the Gram matrix - let mut gram_to_be = PartialMatrix::new(); - self.constraints.with_untracked(|csts| { - for (_, cst) in csts { - let args = cst.args; - let row = elts[args.0].index; - let col = elts[args.1].index; - gram_to_be.push_sym(row, col, cst.rep); - } - }); - - // set up the initial configuration matrix and the diagonal of the - // Gram matrix - let mut guess_to_be = DMatrix::::zeros(5, elts.len()); - for (_, elt) in elts { - let index = elt.index; - gram_to_be.push_sym(index, index, 1.0); - guess_to_be.set_column(index, &elt.rep); - } - - (gram_to_be, guess_to_be) - }); - - /* DEBUG */ - // log the Gram matrix - console::log_1(&JsValue::from("Gram matrix:")); - gram.log_to_console(); - - /* DEBUG */ - // log the initial configuration matrix - console::log_1(&JsValue::from("old configuration:")); - for j in 0..guess.nrows() { - let mut row_str = String::new(); - for k in 0..guess.ncols() { - row_str.push_str(format!(" {:>8.3}", guess[(j, k)]).as_str()); - } - console::log_1(&JsValue::from(row_str)); - } - - // look for a configuration with the given Gram matrix - let (config, success, history) = realize_gram( - &gram, guess, &[], - 1.0e-12, 0.5, 0.9, 1.1, 200, 110 - ); - - /* DEBUG */ - // report the outcome of the search - console::log_1(&JsValue::from( - if success { - "Target accuracy achieved!" - } else { - "Failed to reach target accuracy" - } - )); - console::log_2(&JsValue::from("Steps:"), &JsValue::from(history.scaled_loss.len() - 1)); - console::log_2(&JsValue::from("Loss:"), &JsValue::from(*history.scaled_loss.last().unwrap())); - - if success { - // read out the solution - self.elements.update(|elts| { - for (_, elt) in elts.iter_mut() { - elt.rep.set_column(0, &config.column(elt.index)); - } - }); - } - } } \ No newline at end of file diff --git a/app-proto/src/engine.rs b/app-proto/src/engine.rs index 2971750..79668bb 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -1,12 +1,4 @@ -use lazy_static::lazy_static; -use nalgebra::{Const, DMatrix, DVector, Dyn}; -use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */ - -// --- elements --- - -pub fn point(x: f64, y: f64, z: f64) -> DVector { - DVector::from_column_slice(&[x, y, z, 0.5, 0.5*(x*x + y*y + z*z)]) -} +use nalgebra::DVector; // the sphere with the given center and radius, with inward-pointing normals pub fn sphere(center_x: f64, center_y: f64, center_z: f64, radius: f64) -> DVector { @@ -32,471 +24,4 @@ pub fn sphere_with_offset(dir_x: f64, dir_y: f64, dir_z: f64, off: f64, curv: f6 0.5 * curv, off * (1.0 + 0.5 * off * curv) ]) -} - -// --- partial matrices --- - -struct MatrixEntry { - index: (usize, usize), - value: f64 -} - -pub struct PartialMatrix(Vec); - -impl PartialMatrix { - pub fn new() -> PartialMatrix { - PartialMatrix(Vec::::new()) - } - - pub fn push_sym(&mut self, row: usize, col: usize, value: f64) { - let PartialMatrix(entries) = self; - entries.push(MatrixEntry { index: (row, col), value: value }); - if row != col { - entries.push(MatrixEntry { index: (col, row), value: value }); - } - } - - /* DEBUG */ - pub fn log_to_console(&self) { - let PartialMatrix(entries) = self; - for ent in entries { - let ent_str = format!("{} {} {}", ent.index.0, ent.index.1, ent.value); - console::log_1(&JsValue::from(ent_str.as_str())); - } - } - - fn proj(&self, a: &DMatrix) -> DMatrix { - let mut result = DMatrix::::zeros(a.nrows(), a.ncols()); - let PartialMatrix(entries) = self; - for ent in entries { - result[ent.index] = a[ent.index]; - } - result - } - - fn sub_proj(&self, rhs: &DMatrix) -> DMatrix { - let mut result = DMatrix::::zeros(rhs.nrows(), rhs.ncols()); - let PartialMatrix(entries) = self; - for ent in entries { - result[ent.index] = ent.value - rhs[ent.index]; - } - result - } -} - -// --- descent history --- - -pub struct DescentHistory { - pub config: Vec>, - pub scaled_loss: Vec, - pub neg_grad: Vec>, - pub min_eigval: Vec, - pub base_step: Vec>, - pub backoff_steps: Vec -} - -impl DescentHistory { - fn new() -> DescentHistory { - DescentHistory { - config: Vec::>::new(), - scaled_loss: Vec::::new(), - neg_grad: Vec::>::new(), - min_eigval: Vec::::new(), - base_step: Vec::>::new(), - backoff_steps: Vec::::new(), - } - } -} - -// --- gram matrix realization --- - -// the Lorentz form -lazy_static! { - static ref Q: DMatrix = DMatrix::from_row_slice(5, 5, &[ - 1.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 1.0, 0.0, 0.0, 0.0, - 0.0, 0.0, 1.0, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, -2.0, - 0.0, 0.0, 0.0, -2.0, 0.0 - ]); -} - -struct SearchState { - config: DMatrix, - err_proj: DMatrix, - loss: f64 -} - -impl SearchState { - fn from_config(gram: &PartialMatrix, config: DMatrix) -> SearchState { - let err_proj = gram.sub_proj(&(config.tr_mul(&*Q) * &config)); - let loss = err_proj.norm_squared(); - SearchState { - config: config, - err_proj: err_proj, - loss: loss - } - } -} - -fn basis_matrix(index: (usize, usize), nrows: usize, ncols: usize) -> DMatrix { - let mut result = DMatrix::::zeros(nrows, ncols); - result[index] = 1.0; - result -} - -// use backtracking line search to find a better configuration -fn seek_better_config( - gram: &PartialMatrix, - state: &SearchState, - base_step: &DMatrix, - base_target_improvement: f64, - min_efficiency: f64, - backoff: f64, - max_backoff_steps: i32 -) -> Option<(SearchState, i32)> { - let mut rate = 1.0; - for backoff_steps in 0..max_backoff_steps { - let trial_config = &state.config + rate * base_step; - let trial_state = SearchState::from_config(gram, trial_config); - let improvement = state.loss - trial_state.loss; - if improvement >= min_efficiency * rate * base_target_improvement { - return Some((trial_state, backoff_steps)); - } - rate *= backoff; - } - None -} - -// seek a matrix `config` for which `config' * Q * config` matches the partial -// matrix `gram`. use gradient descent starting from `guess` -pub fn realize_gram( - gram: &PartialMatrix, - guess: DMatrix, - frozen: &[(usize, usize)], - scaled_tol: f64, - min_efficiency: f64, - backoff: f64, - reg_scale: f64, - max_descent_steps: i32, - max_backoff_steps: i32 -) -> (DMatrix, bool, DescentHistory) { - // start the descent history - let mut history = DescentHistory::new(); - - // find the dimension of the search space - let element_dim = guess.nrows(); - let assembly_dim = guess.ncols(); - let total_dim = element_dim * assembly_dim; - - // scale the tolerance - let scale_adjustment = (gram.0.len() as f64).sqrt(); - let tol = scale_adjustment * scaled_tol; - - // convert the frozen indices to stacked format - let frozen_stacked: Vec = frozen.into_iter().map( - |index| index.1*element_dim + index.0 - ).collect(); - - // use Newton's method with backtracking and gradient descent backup - let mut state = SearchState::from_config(gram, guess); - for _ in 0..max_descent_steps { - // stop if the loss is tolerably low - history.config.push(state.config.clone()); - history.scaled_loss.push(state.loss / scale_adjustment); - if state.loss < tol { break; } - - // find the negative gradient of the loss function - let neg_grad = 4.0 * &*Q * &state.config * &state.err_proj; - let mut neg_grad_stacked = neg_grad.clone().reshape_generic(Dyn(total_dim), Const::<1>); - history.neg_grad.push(neg_grad.clone()); - - // find the negative Hessian of the loss function - let mut hess_cols = Vec::>::with_capacity(total_dim); - for col in 0..assembly_dim { - for row in 0..element_dim { - let index = (row, col); - let basis_mat = basis_matrix(index, element_dim, assembly_dim); - let neg_d_err = - basis_mat.tr_mul(&*Q) * &state.config - + state.config.tr_mul(&*Q) * &basis_mat; - let neg_d_err_proj = gram.proj(&neg_d_err); - let deriv_grad = 4.0 * &*Q * ( - -&basis_mat * &state.err_proj - + &state.config * &neg_d_err_proj - ); - hess_cols.push(deriv_grad.reshape_generic(Dyn(total_dim), Const::<1>)); - } - } - let mut hess = DMatrix::from_columns(hess_cols.as_slice()); - - // regularize the Hessian - let min_eigval = hess.symmetric_eigenvalues().min(); - if min_eigval <= 0.0 { - hess -= reg_scale * min_eigval * DMatrix::identity(total_dim, total_dim); - } - history.min_eigval.push(min_eigval); - - // project the negative gradient and negative Hessian onto the - // orthogonal complement of the frozen subspace - let zero_col = DVector::zeros(total_dim); - let zero_row = zero_col.transpose(); - for &k in &frozen_stacked { - neg_grad_stacked[k] = 0.0; - hess.set_row(k, &zero_row); - hess.set_column(k, &zero_col); - hess[(k, k)] = 1.0; - } - - // 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)); - history.base_step.push(base_step.clone()); - - // use backtracking line search to find a better configuration - match seek_better_config( - gram, &state, &base_step, neg_grad.dot(&base_step), - min_efficiency, backoff, max_backoff_steps - ) { - Some((better_state, backoff_steps)) => { - state = better_state; - history.backoff_steps.push(backoff_steps); - }, - None => return (state.config, false, history) - }; - } - (state.config, state.loss < tol, history) -} - -// --- tests --- - -#[cfg(test)] -mod tests { - use std::{array, f64::consts::PI}; - - use super::*; - - #[test] - fn sub_proj_test() { - let target = PartialMatrix(vec![ - 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, - 4.0, 5.0, 6.0 - ]); - let expected_result = DMatrix::::from_row_slice(2, 3, &[ - 18.0, 0.0, 36.0, - 0.0, 54.0, 63.0 - ]); - assert_eq!(target.sub_proj(&attempt), expected_result); - } - - #[test] - fn zero_loss_test() { - let gram = PartialMatrix({ - let mut entries = Vec::::new(); - for j in 0..3 { - for k in 0..3 { - entries.push(MatrixEntry { - index: (j, k), - value: if j == k { 1.0 } else { -1.0 } - }); - } - } - entries - }); - let config = { - let a: f64 = 0.75_f64.sqrt(); - DMatrix::from_columns(&[ - sphere(1.0, 0.0, 0.0, a), - sphere(-0.5, a, 0.0, a), - sphere(-0.5, -a, 0.0, a) - ]) - }; - let state = SearchState::from_config(&gram, config); - 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, history) = realize_gram( - &gram, guess, &frozen, - SCALED_TOL, 0.5, 0.9, 1.1, 200, 110 - ); - 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!((config[(3, k)] - 1.0 / diam).abs() < entry_tol); - } - print!("\nCompleted Gram matrix:{}", config.tr_mul(&*Q) * &config); - if success { - println!("Target accuracy achieved!"); - } else { - println!("Failed to reach target accuracy"); - } - println!("Steps: {}", history.scaled_loss.len() - 1); - println!("Loss: {}", history.scaled_loss.last().unwrap()); - if success { - println!("\nChain diameters:"); - println!(" {} sun (given)", 1.0 / config[(3, 3)]); - for k in 4..9 { - println!(" {} sun", 1.0 / config[(3, k)]); - } - } - println!("\nStep │ Loss\n─────┼────────────────────────────────"); - for (step, scaled_loss) in history.scaled_loss.into_iter().enumerate() { - println!("{:<4} │ {}", step, scaled_loss); - } - } - - // --- process inspection examples --- - - // these tests are meant for human inspection, not automated use. run them - // one at a time in `--nocapture` mode and read through the results and - // optimization histories that they print out. the `run-examples` script - // will run all of them - - #[test] - fn three_spheres_example() { - let gram = PartialMatrix({ - let mut entries = Vec::::new(); - for j in 0..3 { - for k in 0..3 { - entries.push(MatrixEntry { - index: (j, k), - value: if j == k { 1.0 } else { -1.0 } - }); - } - } - entries - }); - let guess = { - let a: f64 = 0.75_f64.sqrt(); - DMatrix::from_columns(&[ - sphere(1.0, 0.0, 0.0, 1.0), - sphere(-0.5, a, 0.0, 1.0), - sphere(-0.5, -a, 0.0, 1.0) - ]) - }; - println!(); - let (config, success, history) = realize_gram( - &gram, guess, &[], - 1.0e-12, 0.5, 0.9, 1.1, 200, 110 - ); - print!("\nCompleted Gram matrix:{}", config.tr_mul(&*Q) * &config); - if success { - println!("Target accuracy achieved!"); - } else { - println!("Failed to reach target accuracy"); - } - println!("Steps: {}", history.scaled_loss.len() - 1); - println!("Loss: {}", history.scaled_loss.last().unwrap()); - println!("\nStep │ Loss\n─────┼────────────────────────────────"); - for (step, scaled_loss) in history.scaled_loss.into_iter().enumerate() { - println!("{:<4} │ {}", step, scaled_loss); - } - } - - #[test] - fn point_on_sphere_example() { - let gram = PartialMatrix({ - let mut entries = Vec::::new(); - for j in 0..2 { - for k in 0..2 { - entries.push(MatrixEntry { - index: (j, k), - value: if (j, k) == (1, 1) { 1.0 } else { 0.0 } - }); - } - } - entries - }); - let guess = DMatrix::from_columns(&[ - point(0.0, 0.0, 2.0), - sphere(0.0, 0.0, 0.0, 1.0) - ]); - let frozen = [(3, 0)]; - println!(); - let (config, success, history) = realize_gram( - &gram, guess, &frozen, - 1.0e-12, 0.5, 0.9, 1.1, 200, 110 - ); - print!("\nCompleted 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: {}", history.scaled_loss.last().unwrap()); - println!("\nStep │ Loss\n─────┼────────────────────────────────"); - for (step, scaled_loss) in history.scaled_loss.into_iter().enumerate() { - println!("{:<4} │ {}", step, scaled_loss); - } - } } \ No newline at end of file diff --git a/engine-proto/gram-test/Engine.jl b/engine-proto/gram-test/Engine.jl index 6dfb6e9..22f5914 100644 --- a/engine-proto/gram-test/Engine.jl +++ b/engine-proto/gram-test/Engine.jl @@ -8,8 +8,7 @@ using Optim export rand_on_shell, Q, DescentHistory, - realize_gram_gradient, realize_gram_newton, realize_gram_optim, - realize_gram_alt_proj, realize_gram + realize_gram_gradient, realize_gram_newton, realize_gram_optim, realize_gram # === guessing === @@ -144,7 +143,7 @@ function realize_gram_gradient( break end - # find the negative gradient of the loss function + # find negative gradient of loss function neg_grad = 4*Q*L*Δ_proj slope = norm(neg_grad) dir = neg_grad / slope @@ -233,7 +232,7 @@ function realize_gram_newton( break end - # find the negative gradient of the loss function + # find the negative gradient of loss function neg_grad = 4*Q*L*Δ_proj # find the negative Hessian of the loss function @@ -314,129 +313,6 @@ function realize_gram_optim( ) end -# seek a matrix `L` for which `L'QL` matches the sparse matrix `gram` at every -# explicit entry of `gram`. use gradient descent starting from `guess`, with an -# alternate technique for finding the projected base step from the unprojected -# Hessian -function realize_gram_alt_proj( - gram::SparseMatrixCSC{T, <:Any}, - guess::Matrix{T}, - frozen = CartesianIndex[]; - scaled_tol = 1e-30, - min_efficiency = 0.5, - backoff = 0.9, - reg_scale = 1.1, - max_descent_steps = 200, - max_backoff_steps = 110 -) where T <: Number - # start history - history = DescentHistory{T}() - - # find the dimension of the search space - dims = size(guess) - element_dim, construction_dim = dims - total_dim = element_dim * construction_dim - - # list the constrained entries of the gram matrix - J, K, _ = findnz(gram) - constrained = zip(J, K) - - # scale the tolerance - scale_adjustment = sqrt(T(length(constrained))) - tol = scale_adjustment * scaled_tol - - # convert the frozen indices to stacked format - frozen_stacked = [(index[2]-1)*element_dim + index[1] for index in frozen] - - # initialize search state - L = copy(guess) - Δ_proj = proj_diff(gram, L'*Q*L) - loss = dot(Δ_proj, Δ_proj) - - # use Newton's method with backtracking and gradient descent backup - for step in 1:max_descent_steps - # stop if the loss is tolerably low - if loss < tol - break - end - - # find the negative gradient of the loss function - neg_grad = 4*Q*L*Δ_proj - - # find the negative Hessian of the loss function - hess = Matrix{T}(undef, total_dim, total_dim) - indices = [(j, k) for k in 1:construction_dim for j in 1:element_dim] - for (j, k) in indices - basis_mat = basis_matrix(T, j, k, dims) - neg_dΔ = basis_mat'*Q*L + L'*Q*basis_mat - neg_dΔ_proj = proj_to_entries(neg_dΔ, constrained) - deriv_grad = 4*Q*(-basis_mat*Δ_proj + L*neg_dΔ_proj) - hess[:, (k-1)*element_dim + j] = reshape(deriv_grad, total_dim) - end - hess_sym = Hermitian(hess) - push!(history.hess, hess_sym) - - # regularize the Hessian - min_eigval = minimum(eigvals(hess_sym)) - push!(history.positive, min_eigval > 0) - if min_eigval <= 0 - hess -= reg_scale * min_eigval * I - end - - # compute the Newton step - neg_grad_stacked = reshape(neg_grad, total_dim) - for k in frozen_stacked - neg_grad_stacked[k] = 0 - hess[k, :] .= 0 - hess[:, k] .= 0 - hess[k, k] = 1 - end - base_step_stacked = Hermitian(hess) \ neg_grad_stacked - base_step = reshape(base_step_stacked, dims) - push!(history.base_step, base_step) - - # store the current position, loss, and slope - L_last = L - loss_last = loss - push!(history.scaled_loss, loss / scale_adjustment) - push!(history.neg_grad, neg_grad) - push!(history.slope, norm(neg_grad)) - - # find a good step size using backtracking line search - push!(history.stepsize, 0) - push!(history.backoff_steps, max_backoff_steps) - empty!(history.last_line_L) - empty!(history.last_line_loss) - rate = one(T) - step_success = false - base_target_improvement = dot(neg_grad, base_step) - for backoff_steps in 0:max_backoff_steps - history.stepsize[end] = rate - L = L_last + rate * base_step - Δ_proj = proj_diff(gram, L'*Q*L) - loss = dot(Δ_proj, Δ_proj) - improvement = loss_last - loss - push!(history.last_line_L, L) - push!(history.last_line_loss, loss / scale_adjustment) - if improvement >= min_efficiency * rate * base_target_improvement - history.backoff_steps[end] = backoff_steps - step_success = true - break - end - rate *= backoff - end - - # if we've hit a wall, quit - if !step_success - return L_last, false, history - end - end - - # return the factorization and its history - push!(history.scaled_loss, loss / scale_adjustment) - L, loss < tol, history -end - # seek a matrix `L` for which `L'QL` matches the sparse matrix `gram` at every # explicit entry of `gram`. use gradient descent starting from `guess` function realize_gram( @@ -445,6 +321,7 @@ function realize_gram( frozen = nothing; scaled_tol = 1e-30, min_efficiency = 0.5, + init_rate = 1.0, backoff = 0.9, reg_scale = 1.1, max_descent_steps = 200, @@ -475,19 +352,20 @@ function realize_gram( unfrozen_stacked = reshape(is_unfrozen, total_dim) end - # initialize search state + # initialize variables + grad_rate = init_rate L = copy(guess) - Δ_proj = proj_diff(gram, L'*Q*L) - loss = dot(Δ_proj, Δ_proj) # use Newton's method with backtracking and gradient descent backup + Δ_proj = proj_diff(gram, L'*Q*L) + loss = dot(Δ_proj, Δ_proj) for step in 1:max_descent_steps # stop if the loss is tolerably low if loss < tol break end - # find the negative gradient of the loss function + # find the negative gradient of loss function neg_grad = 4*Q*L*Δ_proj # find the negative Hessian of the loss function @@ -542,7 +420,6 @@ function realize_gram( empty!(history.last_line_loss) rate = one(T) step_success = false - base_target_improvement = dot(neg_grad, base_step) for backoff_steps in 0:max_backoff_steps history.stepsize[end] = rate L = L_last + rate * base_step @@ -551,7 +428,7 @@ function realize_gram( improvement = loss_last - loss push!(history.last_line_L, L) push!(history.last_line_loss, loss / scale_adjustment) - if improvement >= min_efficiency * rate * base_target_improvement + if improvement >= min_efficiency * rate * dot(neg_grad, base_step) history.backoff_steps[end] = backoff_steps step_success = true break diff --git a/engine-proto/gram-test/irisawa-hexlet.jl b/engine-proto/gram-test/irisawa-hexlet.jl index 607db61..67def8c 100644 --- a/engine-proto/gram-test/irisawa-hexlet.jl +++ b/engine-proto/gram-test/irisawa-hexlet.jl @@ -74,13 +74,4 @@ if success for k in 5:9 println(" ", 1 / L[4,k], " sun") end -end - -# test an alternate technique for finding the projected base step from the -# unprojected Hessian -L_alt, success_alt, history_alt = Engine.realize_gram_alt_proj(gram, guess, frozen) -completed_gram_alt = L_alt'*Engine.Q*L_alt -println("\nDifference in result using alternate projection:\n") -display(completed_gram_alt - completed_gram) -println("\nDifference in steps: ", size(history_alt.scaled_loss, 1) - size(history.scaled_loss, 1)) -println("Difference in loss: ", history_alt.scaled_loss[end] - history.scaled_loss[end], "\n") \ No newline at end of file +end \ No newline at end of file diff --git a/engine-proto/gram-test/sphere-in-tetrahedron.jl b/engine-proto/gram-test/sphere-in-tetrahedron.jl index 5d479cf..97f0720 100644 --- a/engine-proto/gram-test/sphere-in-tetrahedron.jl +++ b/engine-proto/gram-test/sphere-in-tetrahedron.jl @@ -64,13 +64,4 @@ else println("\nFailed to reach target accuracy") end println("Steps: ", size(history.scaled_loss, 1)) -println("Loss: ", history.scaled_loss[end], "\n") - -# test an alternate technique for finding the projected base step from the -# unprojected Hessian -L_alt, success_alt, history_alt = Engine.realize_gram_alt_proj(gram, guess, frozen) -completed_gram_alt = L_alt'*Engine.Q*L_alt -println("\nDifference in result using alternate projection:\n") -display(completed_gram_alt - completed_gram) -println("\nDifference in steps: ", size(history_alt.scaled_loss, 1) - size(history.scaled_loss, 1)) -println("Difference in loss: ", history_alt.scaled_loss[end] - history.scaled_loss[end], "\n") \ No newline at end of file +println("Loss: ", history.scaled_loss[end], "\n") \ No newline at end of file diff --git a/engine-proto/gram-test/tetrahedron-radius-ratio.jl b/engine-proto/gram-test/tetrahedron-radius-ratio.jl index 9fec28e..7ceb794 100644 --- a/engine-proto/gram-test/tetrahedron-radius-ratio.jl +++ b/engine-proto/gram-test/tetrahedron-radius-ratio.jl @@ -93,13 +93,4 @@ if success infty = BigFloat[0, 0, 0, 0, 1] radius_ratio = dot(infty, Engine.Q * L[:,5]) / dot(infty, Engine.Q * L[:,6]) println("\nCircumradius / inradius: ", radius_ratio) -end - -# test an alternate technique for finding the projected base step from the -# unprojected Hessian -L_alt, success_alt, history_alt = Engine.realize_gram_alt_proj(gram, guess, frozen) -completed_gram_alt = L_alt'*Engine.Q*L_alt -println("\nDifference in result using alternate projection:\n") -display(completed_gram_alt - completed_gram) -println("\nDifference in steps: ", size(history_alt.scaled_loss, 1) - size(history.scaled_loss, 1)) -println("Difference in loss: ", history_alt.scaled_loss[end] - history.scaled_loss[end], "\n") \ No newline at end of file +end \ No newline at end of file