From 16df161fe761e7c550282ff02916e6f089ff17b6 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 24 Oct 2024 19:51:10 -0700 Subject: [PATCH 01/24] Test alternate projection technique --- engine-proto/gram-test/Engine.jl | 133 +++++++++++++++++- engine-proto/gram-test/irisawa-hexlet.jl | 11 +- .../gram-test/sphere-in-tetrahedron.jl | 11 +- .../gram-test/tetrahedron-radius-ratio.jl | 11 +- 4 files changed, 159 insertions(+), 7 deletions(-) diff --git a/engine-proto/gram-test/Engine.jl b/engine-proto/gram-test/Engine.jl index 22f5914..1eb72f5 100644 --- a/engine-proto/gram-test/Engine.jl +++ b/engine-proto/gram-test/Engine.jl @@ -8,7 +8,8 @@ using Optim export rand_on_shell, Q, DescentHistory, - realize_gram_gradient, realize_gram_newton, realize_gram_optim, realize_gram + realize_gram_gradient, realize_gram_newton, realize_gram_optim, + realize_gram_alt_proj, realize_gram # === guessing === @@ -143,7 +144,7 @@ function realize_gram_gradient( break end - # find negative gradient of loss function + # find the negative gradient of the loss function neg_grad = 4*Q*L*Δ_proj slope = norm(neg_grad) dir = neg_grad / slope @@ -232,7 +233,7 @@ function realize_gram_newton( break end - # find the negative gradient of loss function + # find the negative gradient of the loss function neg_grad = 4*Q*L*Δ_proj # find the negative Hessian of the loss function @@ -313,6 +314,130 @@ 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, + init_rate = 1.0, + 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 variables + grad_rate = init_rate + L = copy(guess) + + # 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 + 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 + 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 * dot(neg_grad, base_step) + 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( @@ -365,7 +490,7 @@ function realize_gram( break end - # find the negative gradient of loss function + # find the negative gradient of the loss function neg_grad = 4*Q*L*Δ_proj # find the negative Hessian of the loss function diff --git a/engine-proto/gram-test/irisawa-hexlet.jl b/engine-proto/gram-test/irisawa-hexlet.jl index 67def8c..607db61 100644 --- a/engine-proto/gram-test/irisawa-hexlet.jl +++ b/engine-proto/gram-test/irisawa-hexlet.jl @@ -74,4 +74,13 @@ if success for k in 5:9 println(" ", 1 / L[4,k], " sun") end -end \ No newline at end of file +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 diff --git a/engine-proto/gram-test/sphere-in-tetrahedron.jl b/engine-proto/gram-test/sphere-in-tetrahedron.jl index 97f0720..5d479cf 100644 --- a/engine-proto/gram-test/sphere-in-tetrahedron.jl +++ b/engine-proto/gram-test/sphere-in-tetrahedron.jl @@ -64,4 +64,13 @@ else println("\nFailed to reach target accuracy") end println("Steps: ", size(history.scaled_loss, 1)) -println("Loss: ", history.scaled_loss[end], "\n") \ No newline at end of file +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 diff --git a/engine-proto/gram-test/tetrahedron-radius-ratio.jl b/engine-proto/gram-test/tetrahedron-radius-ratio.jl index 7ceb794..9fec28e 100644 --- a/engine-proto/gram-test/tetrahedron-radius-ratio.jl +++ b/engine-proto/gram-test/tetrahedron-radius-ratio.jl @@ -93,4 +93,13 @@ 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 \ No newline at end of file +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 -- 2.34.1 From e59d60bf7745dfb25ca054a5beecb08bc8f29b2f Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Fri, 25 Oct 2024 17:17:49 -0700 Subject: [PATCH 02/24] Reorganize search state; remove unused variables --- engine-proto/gram-test/Engine.jl | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/engine-proto/gram-test/Engine.jl b/engine-proto/gram-test/Engine.jl index 1eb72f5..6dfb6e9 100644 --- a/engine-proto/gram-test/Engine.jl +++ b/engine-proto/gram-test/Engine.jl @@ -324,7 +324,6 @@ function realize_gram_alt_proj( frozen = CartesianIndex[]; scaled_tol = 1e-30, min_efficiency = 0.5, - init_rate = 1.0, backoff = 0.9, reg_scale = 1.1, max_descent_steps = 200, @@ -349,13 +348,12 @@ function realize_gram_alt_proj( # convert the frozen indices to stacked format frozen_stacked = [(index[2]-1)*element_dim + index[1] for index in frozen] - # initialize variables - grad_rate = init_rate + # initialize search state L = copy(guess) - - # use Newton's method with backtracking and gradient descent backup Δ_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 @@ -411,6 +409,7 @@ function realize_gram_alt_proj( 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 @@ -419,7 +418,7 @@ function realize_gram_alt_proj( improvement = loss_last - loss push!(history.last_line_L, L) push!(history.last_line_loss, loss / scale_adjustment) - if improvement >= min_efficiency * rate * dot(neg_grad, base_step) + if improvement >= min_efficiency * rate * base_target_improvement history.backoff_steps[end] = backoff_steps step_success = true break @@ -446,7 +445,6 @@ 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, @@ -477,13 +475,12 @@ function realize_gram( unfrozen_stacked = reshape(is_unfrozen, total_dim) end - # initialize variables - grad_rate = init_rate + # initialize search state L = copy(guess) - - # use Newton's method with backtracking and gradient descent backup Δ_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 @@ -545,6 +542,7 @@ 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 @@ -553,7 +551,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 * dot(neg_grad, base_step) + if improvement >= min_efficiency * rate * base_target_improvement history.backoff_steps[end] = backoff_steps step_success = true break -- 2.34.1 From 9fe03264ab15201180c6759d03fe6094d0a25eac Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Fri, 25 Oct 2024 17:34:29 -0700 Subject: [PATCH 03/24] Port the Gram matrix realization routine to Rust Validate with the process inspection example tests, which print out their results and optimization histories when run one at a time in `--nocapture` mode. --- app-proto/Cargo.toml | 1 + app-proto/run-examples | 7 + app-proto/src/engine.rs | 321 +++++++++++++++++++++++++++++++++++++++- 3 files changed, 328 insertions(+), 1 deletion(-) create mode 100755 app-proto/run-examples diff --git a/app-proto/Cargo.toml b/app-proto/Cargo.toml index 920469a..04ea271 100644 --- a/app-proto/Cargo.toml +++ b/app-proto/Cargo.toml @@ -10,6 +10,7 @@ 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 new file mode 100755 index 0000000..b5c3de5 --- /dev/null +++ b/app-proto/run-examples @@ -0,0 +1,7 @@ +# 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::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 79668bb..f5ad19a 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -1,4 +1,11 @@ -use nalgebra::DVector; +use lazy_static::lazy_static; +use nalgebra::{Const, DMatrix, DVector, Dyn}; + +// --- 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)]) +} // 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 { @@ -24,4 +31,316 @@ 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), + val: f64 +} + +struct PartialMatrix(Vec); + +impl PartialMatrix { + 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.val - rhs[ent.index]; + } + result + } +} + +// --- 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 { + let mut rate = 1.0; + for _ 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); + } + rate *= backoff; + } + None +} + +// seek a matrix `config` for which `config' * Q * config` matches the partial +// matrix `gram`. use gradient descent starting from `guess` +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) { + // 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 = ((guess.ncols() - frozen.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 + println!("loss: {}", state.loss); + /*println!("projected error: {}", state.err_proj);*/ + 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>); + + // 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); + } + + // 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 + 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)); + + // 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) => state = better_state, + None => return (state.config, false) + }; + } + (state.config, state.loss < tol) +} + +// --- tests --- + +#[cfg(test)] +mod tests { + use std::f64; + + 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 } + ]); + 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), + val: 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); + } + + // --- 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), + val: 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) = realize_gram( + &gram, guess, &[], + 1.0e-12, 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); + } + + #[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), + val: 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) + ]); + println!(); + let (config, success) = realize_gram( + &gram, guess, &[(3, 0)], + 1.0e-12, 0.5, 0.9, 1.1, 200, 110 + ); + print!("\nCompleted Gram matrix:{}", config.tr_mul(&*Q) * &config); + print!("Configuration:{}", 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); + } } \ No newline at end of file -- 2.34.1 From 9f8632efb34977393ccb9fc315317f5d69b81bf4 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Fri, 25 Oct 2024 21:43:53 -0700 Subject: [PATCH 04/24] Port the Irisawa hexlet test to Rust In the process, notice that the tolerance scale adjustment was ported wrong, and correct it. --- app-proto/run-examples | 1 + app-proto/src/engine.rs | 118 +++++++++++++++++++++++++++++++++++----- 2 files changed, 105 insertions(+), 14 deletions(-) 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); -- 2.34.1 From ce33bbf41877094c6a2985e6b34f73ff50c55bd1 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sat, 26 Oct 2024 00:02:06 -0700 Subject: [PATCH 05/24] Record optimization history --- app-proto/src/engine.rs | 94 ++++++++++++++++++++++++++++++----------- 1 file changed, 69 insertions(+), 25 deletions(-) diff --git a/app-proto/src/engine.rs b/app-proto/src/engine.rs index 6f2952c..a4d89a1 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -62,6 +62,30 @@ impl PartialMatrix { } } +// --- descent history --- + +struct DescentHistory { + config: Vec>, + scaled_loss: Vec, + neg_grad: Vec>, + min_eigval: Vec, + base_step: Vec>, + 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 @@ -108,14 +132,14 @@ fn seek_better_config( min_efficiency: f64, backoff: f64, max_backoff_steps: i32 -) -> Option { +) -> Option<(SearchState, i32)> { let mut rate = 1.0; - for _ in 0..max_backoff_steps { + 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); + return Some((trial_state, backoff_steps)); } rate *= backoff; } @@ -134,7 +158,10 @@ fn realize_gram( reg_scale: f64, max_descent_steps: i32, max_backoff_steps: i32 -) -> (DMatrix, bool) { +) -> (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(); @@ -153,13 +180,14 @@ 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!("scaled loss: {}", state.loss / scale_adjustment); - /* println!("projected error: {}", state.err_proj); */ + 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); @@ -182,10 +210,10 @@ 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); } + history.min_eigval.push(min_eigval); // project the negative gradient and negative Hessian onto the // orthogonal complement of the frozen subspace @@ -207,17 +235,21 @@ fn realize_gram( */ 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) => state = better_state, - None => return (state.config, false) + 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) + (state.config, state.loss < tol, history) } // --- tests --- @@ -329,29 +361,33 @@ mod tests { ); let frozen: [(usize, usize); 4] = array::from_fn(|k| (3, k)); const SCALED_TOL: f64 = 1.0e-12; - let (config, success) = realize_gram( + 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); - 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); + println!("Steps: {}", history.scaled_loss.len() - 1); + println!("Loss: {}", history.scaled_loss.last().unwrap()); if success { println!("\nChain diameters:"); - println!(" {} sun (given)", 1.0 / final_state.config[(3, 3)]); + println!(" {} sun (given)", 1.0 / config[(3, 3)]); for k in 4..9 { - println!(" {} sun", 1.0 / final_state.config[(3, k)]); + println!(" {} sun", 1.0 / 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); + println!("\nStep │ Loss\n─────┼────────────────────────────────"); + for (step, scaled_loss) in history.scaled_loss.into_iter().enumerate() { + println!("{:<4} │ {}", step, scaled_loss); } } @@ -385,18 +421,22 @@ mod tests { ]) }; println!(); - let (config, success) = realize_gram( + 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); - 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); + 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] @@ -419,18 +459,22 @@ mod tests { ]); let frozen = [(3, 0)]; println!(); - let (config, success) = realize_gram( + 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); - 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); + 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 -- 2.34.1 From a37c71153d198e64895ced18e71cd39ebee3bc64 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sat, 26 Oct 2024 23:51:27 -0700 Subject: [PATCH 06/24] Enforce constraints in the editor --- app-proto/src/add_remove.rs | 53 +++++++++++---------- app-proto/src/assembly.rs | 94 +++++++++++++++++++++++++++++++++++-- app-proto/src/engine.rs | 40 ++++++++++++---- 3 files changed, 151 insertions(+), 36 deletions(-) diff --git a/app-proto/src/add_remove.rs b/app-proto/src/add_remove.rs index ab5db70..7066089 100644 --- a/app-proto/src/add_remove.rs +++ b/app-proto/src/add_remove.rs @@ -12,7 +12,8 @@ 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() + constraints: BTreeSet::default(), + index: 0 } ); let _ = assembly.try_insert_element( @@ -21,7 +22,8 @@ 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() + constraints: BTreeSet::default(), + index: 0 } ); let _ = assembly.try_insert_element( @@ -30,7 +32,8 @@ 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() + constraints: BTreeSet::default(), + index: 0 } ); let _ = assembly.try_insert_element( @@ -39,7 +42,8 @@ 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() + constraints: BTreeSet::default(), + index: 0 } ); let _ = assembly.try_insert_element( @@ -48,7 +52,8 @@ 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() + constraints: BTreeSet::default(), + index: 0 } ); let _ = assembly.try_insert_element( @@ -57,17 +62,8 @@ 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() - } - ); - 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) + constraints: BTreeSet::default(), + index: 0 } ); } @@ -81,7 +77,8 @@ 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() + constraints: BTreeSet::default(), + index: 0 } ); let _ = assembly.try_insert_element( @@ -90,7 +87,8 @@ 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() + constraints: BTreeSet::default(), + index: 0 } ); let _ = assembly.try_insert_element( @@ -99,7 +97,8 @@ 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() + constraints: BTreeSet::default(), + index: 0 } ); let _ = assembly.try_insert_element( @@ -108,7 +107,8 @@ 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() + constraints: BTreeSet::default(), + index: 0 } ); let _ = assembly.try_insert_element( @@ -117,7 +117,8 @@ 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() + constraints: BTreeSet::default(), + index: 0 } ); let _ = assembly.try_insert_element( @@ -126,7 +127,8 @@ 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() + constraints: BTreeSet::default(), + index: 0 } ); let _ = assembly.try_insert_element( @@ -135,7 +137,8 @@ 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() + constraints: BTreeSet::default(), + index: 0 } ); let _ = assembly.try_insert_element( @@ -144,7 +147,8 @@ 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() + constraints: BTreeSet::default(), + index: 0 } ); } @@ -215,6 +219,7 @@ 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 e8dab79..228357e 100644 --- a/app-proto/src/assembly.rs +++ b/app-proto/src/assembly.rs @@ -1,8 +1,11 @@ -use nalgebra::DVector; +use nalgebra::{DMatrix, 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 { @@ -10,7 +13,10 @@ pub struct Element { pub label: String, pub color: [f32; 3], pub rep: DVector, - pub constraints: BTreeSet + pub constraints: BTreeSet, + + // internal properties, not reflected in any view + pub index: usize } #[derive(Clone)] @@ -40,6 +46,8 @@ 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 @@ -77,7 +85,8 @@ 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() + constraints: BTreeSet::default(), + index: 0 } ); } @@ -90,4 +99,83 @@ 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 a4d89a1..2971750 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -1,5 +1,6 @@ use lazy_static::lazy_static; use nalgebra::{Const, DMatrix, DVector, Dyn}; +use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */ // --- elements --- @@ -40,9 +41,30 @@ struct MatrixEntry { value: f64 } -struct PartialMatrix(Vec); +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; @@ -64,13 +86,13 @@ impl PartialMatrix { // --- descent history --- -struct DescentHistory { - config: Vec>, - scaled_loss: Vec, - neg_grad: Vec>, - min_eigval: Vec, - base_step: Vec>, - backoff_steps: Vec +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 { @@ -148,7 +170,7 @@ fn seek_better_config( // seek a matrix `config` for which `config' * Q * config` matches the partial // matrix `gram`. use gradient descent starting from `guess` -fn realize_gram( +pub fn realize_gram( gram: &PartialMatrix, guess: DMatrix, frozen: &[(usize, usize)], -- 2.34.1 From e5f4d523f90161a4615c5d4c51d83cb6e2a6cb1c Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Tue, 29 Oct 2024 13:46:15 -0700 Subject: [PATCH 07/24] Update the realization when a constraint is activated Sycamore probably has a better way to do this, but this way works for now. --- app-proto/src/add_remove.rs | 10 +++++++++- app-proto/src/assembly.rs | 12 +++++++----- 2 files changed, 16 insertions(+), 6 deletions(-) diff --git a/app-proto/src/add_remove.rs b/app-proto/src/add_remove.rs index 7066089..00b63f8 100644 --- a/app-proto/src/add_remove.rs +++ b/app-proto/src/add_remove.rs @@ -214,10 +214,11 @@ pub fn AddRemove() -> View { (arg_vec[0].clone(), arg_vec[1].clone()) } ); + let active = create_signal(true); state.assembly.insert_constraint(Constraint { args: args, rep: 0.0, - active: create_signal(true) + active: active }); state.assembly.realize(); state.selection.update(|sel| sel.clear()); @@ -236,6 +237,13 @@ pub fn AddRemove() -> View { ); } }); + + // make constraint activation trigger a realization update + create_effect(move || { + if active.get() { + state.assembly.realize(); + } + }); } ) { "🔗" } select(bind:value=assembly_name) { /* DEBUG */ diff --git a/app-proto/src/assembly.rs b/app-proto/src/assembly.rs index 228357e..648d0ef 100644 --- a/app-proto/src/assembly.rs +++ b/app-proto/src/assembly.rs @@ -97,7 +97,7 @@ impl Assembly { self.elements.update(|elts| { elts[args.0].constraints.insert(key); elts[args.1].constraints.insert(key); - }) + }); } // --- realization --- @@ -116,10 +116,12 @@ impl Assembly { 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); + if cst.active.get_untracked() { + 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); + } } }); -- 2.34.1 From e0880d2ad2d314c06fbc324bcf2bbe8aa7ba0b8d Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Tue, 29 Oct 2024 22:32:00 -0700 Subject: [PATCH 08/24] Make constraints editable --- app-proto/Cargo.toml | 1 + app-proto/main.css | 11 +++++++++-- app-proto/src/add_remove.rs | 17 +++++++++++++---- app-proto/src/assembly.rs | 8 +++++--- app-proto/src/outline.rs | 25 +++++++++++++++++++++---- 5 files changed, 49 insertions(+), 13 deletions(-) diff --git a/app-proto/Cargo.toml b/app-proto/Cargo.toml index 04ea271..e5bc05e 100644 --- a/app-proto/Cargo.toml +++ b/app-proto/Cargo.toml @@ -26,6 +26,7 @@ console_error_panic_hook = { version = "0.1.7", optional = true } version = "0.3.69" features = [ 'HtmlCanvasElement', + 'HtmlInputElement', 'Performance', 'WebGl2RenderingContext', 'WebGlBuffer', diff --git a/app-proto/main.css b/app-proto/main.css index bdbacfb..44dc7a1 100644 --- a/app-proto/main.css +++ b/app-proto/main.css @@ -93,7 +93,7 @@ details[open]:has(li) .elt-switch::after { display: flex; } -.elt-rep > div, .cst-rep { +.elt-rep > div { padding: 2px 0px 0px 0px; font-size: 10pt; text-align: center; @@ -104,10 +104,17 @@ details[open]:has(li) .elt-switch::after { font-style: italic; } -.cst > input { +.cst > input[type=checkbox] { margin: 0px 8px 0px 0px; } +.cst > input[type=number] { + color: #fcfcfc; + background-color: inherit; + border: 1px solid #555; + border-radius: 2px; +} + /* display */ canvas { diff --git a/app-proto/src/add_remove.rs b/app-proto/src/add_remove.rs index 00b63f8..5435418 100644 --- a/app-proto/src/add_remove.rs +++ b/app-proto/src/add_remove.rs @@ -214,11 +214,14 @@ pub fn AddRemove() -> View { (arg_vec[0].clone(), arg_vec[1].clone()) } ); + let rep = create_signal(0.0); let active = create_signal(true); state.assembly.insert_constraint(Constraint { args: args, - rep: 0.0, - active: active + rep: rep, + rep_text: create_signal(String::new()), + rep_valid: create_signal(false), + active: active, }); state.assembly.realize(); state.selection.update(|sel| sel.clear()); @@ -233,13 +236,19 @@ pub fn AddRemove() -> View { &JsValue::from(cst.args.0), &JsValue::from(cst.args.1), &JsValue::from(":"), - &JsValue::from(cst.rep) + &JsValue::from(cst.rep.get_untracked()) ); } }); - // make constraint activation trigger a realization update + // update the realization when the constraint activated, or + // edited while active create_effect(move || { + rep.track(); + console::log_2( + &JsValue::from("Constraint rep updated to"), + &JsValue::from(rep.get_untracked()) + ); if active.get() { state.assembly.realize(); } diff --git a/app-proto/src/assembly.rs b/app-proto/src/assembly.rs index 648d0ef..62f5405 100644 --- a/app-proto/src/assembly.rs +++ b/app-proto/src/assembly.rs @@ -22,7 +22,9 @@ pub struct Element { #[derive(Clone)] pub struct Constraint { pub args: (usize, usize), - pub rep: f64, + pub rep: Signal, + pub rep_text: Signal, + pub rep_valid: Signal, pub active: Signal } @@ -116,11 +118,11 @@ impl Assembly { let mut gram_to_be = PartialMatrix::new(); self.constraints.with_untracked(|csts| { for (_, cst) in csts { - if cst.active.get_untracked() { + if cst.active.get_untracked() && cst.rep_valid.get_untracked() { 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); + gram_to_be.push_sym(row, col, cst.rep.get_untracked()); } } }); diff --git a/app-proto/src/outline.rs b/app-proto/src/outline.rs index 4e4de9c..d6b0390 100644 --- a/app-proto/src/outline.rs +++ b/app-proto/src/outline.rs @@ -1,6 +1,7 @@ use itertools::Itertools; use sycamore::{prelude::*, web::tags::div}; -use web_sys::{Element, KeyboardEvent, MouseEvent, wasm_bindgen::JsCast}; +use web_sys::{Element, Event, HtmlInputElement, KeyboardEvent, MouseEvent, wasm_bindgen::JsCast}; +use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */ use crate::AppState; @@ -51,8 +52,6 @@ pub fn Outline() -> View { let constrained = elt.constraints.len() > 0; let details_node = create_node_ref(); view! { - /* [TO DO] switch to integer-valued parameters whenever - that becomes possible again */ li { details(ref=details_node) { summary( @@ -138,7 +137,25 @@ pub fn Outline() -> View { li(class="cst") { input(r#type="checkbox", bind:checked=cst.active) div(class="cst-label") { (other_arg_label) } - div(class="cst-rep") { (cst.rep) } + input( + r#type="number", + step="0.01", + bind:value=cst.rep_text, + on:change=move |event: Event| { + let target: HtmlInputElement = event.target().unwrap().unchecked_into(); + let rep_valid = target.check_validity() && !target.value().is_empty(); + batch(|| { + cst.rep_valid.set(rep_valid); + if rep_valid { + console::log_2( + &JsValue::from("Constraint rep parsed to"), + &JsValue::from(target.value_as_number()) + ); + cst.rep.set(target.value_as_number()); + } + }); + } + ) } } }, -- 2.34.1 From a46ef2c8d63eaca964fceea6c1bf547a89944261 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Tue, 29 Oct 2024 22:53:48 -0700 Subject: [PATCH 09/24] Work around data binding bug in number input Setting `bind:value` or `bind:valueAsNumber` for a number input seems to restrict what you can type in it. We work around this by switching to text inputs for now. We should probably switch back to number inputs if we can, though, because they let us take advantage of the browser's parsing and validation. --- app-proto/main.css | 2 +- app-proto/src/outline.rs | 19 +++++++++---------- 2 files changed, 10 insertions(+), 11 deletions(-) diff --git a/app-proto/main.css b/app-proto/main.css index 44dc7a1..32ae5bf 100644 --- a/app-proto/main.css +++ b/app-proto/main.css @@ -108,7 +108,7 @@ details[open]:has(li) .elt-switch::after { margin: 0px 8px 0px 0px; } -.cst > input[type=number] { +.cst > input[type=text] { color: #fcfcfc; background-color: inherit; border: 1px solid #555; diff --git a/app-proto/src/outline.rs b/app-proto/src/outline.rs index d6b0390..62bc529 100644 --- a/app-proto/src/outline.rs +++ b/app-proto/src/outline.rs @@ -138,22 +138,21 @@ pub fn Outline() -> View { input(r#type="checkbox", bind:checked=cst.active) div(class="cst-label") { (other_arg_label) } input( - r#type="number", - step="0.01", + r#type="text", bind:value=cst.rep_text, on:change=move |event: Event| { let target: HtmlInputElement = event.target().unwrap().unchecked_into(); - let rep_valid = target.check_validity() && !target.value().is_empty(); - batch(|| { - cst.rep_valid.set(rep_valid); - if rep_valid { + match target.value().parse::() { + Ok(rep) => batch(|| { + cst.rep.set(rep); + cst.rep_valid.set(true); console::log_2( &JsValue::from("Constraint rep parsed to"), - &JsValue::from(target.value_as_number()) + &JsValue::from(rep) ); - cst.rep.set(target.value_as_number()); - } - }); + }), + Err(_) => cst.rep_valid.set(false) + }; } ) } -- 2.34.1 From 76ad4245d5346c4e6bafac5eeee2544608452723 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Tue, 29 Oct 2024 23:43:41 -0700 Subject: [PATCH 10/24] Factor out Lorentz product input --- app-proto/src/outline.rs | 52 +++++++++++++++++++++++----------------- 1 file changed, 30 insertions(+), 22 deletions(-) diff --git a/app-proto/src/outline.rs b/app-proto/src/outline.rs index 62bc529..fcf983f 100644 --- a/app-proto/src/outline.rs +++ b/app-proto/src/outline.rs @@ -3,11 +3,36 @@ use sycamore::{prelude::*, web::tags::div}; use web_sys::{Element, Event, HtmlInputElement, KeyboardEvent, MouseEvent, wasm_bindgen::JsCast}; use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */ -use crate::AppState; +use crate::{AppState, assembly::Constraint}; -// this component lists the elements of the assembly, showing the constraints -// on each element as a collapsible sub-list. its implementation is based on -// Kate Morley's HTML + CSS tree views: +// an editable view of the Lorentz product representing a constraint +#[component(inline_props)] +fn LorentzProductInput(constraint: Constraint) -> View { + view! { + input( + r#type="text", + bind:value=constraint.rep_text, + on:change=move |event: Event| { + let target: HtmlInputElement = event.target().unwrap().unchecked_into(); + match target.value().parse::() { + Ok(rep) => batch(|| { + constraint.rep.set(rep); + constraint.rep_valid.set(true); + console::log_2( + &JsValue::from("Constraint rep parsed to"), + &JsValue::from(rep) + ); + }), + Err(_) => constraint.rep_valid.set(false) + }; + } + ) + } +} + +// a component that lists the elements of the current assembly, showing the +// constraints on each element as a collapsible sub-list. its implementation +// is based on Kate Morley's HTML + CSS tree views: // // https://iamkate.com/code/tree-views/ // @@ -137,24 +162,7 @@ pub fn Outline() -> View { li(class="cst") { input(r#type="checkbox", bind:checked=cst.active) div(class="cst-label") { (other_arg_label) } - input( - r#type="text", - bind:value=cst.rep_text, - on:change=move |event: Event| { - let target: HtmlInputElement = event.target().unwrap().unchecked_into(); - match target.value().parse::() { - Ok(rep) => batch(|| { - cst.rep.set(rep); - cst.rep_valid.set(true); - console::log_2( - &JsValue::from("Constraint rep parsed to"), - &JsValue::from(rep) - ); - }), - Err(_) => cst.rep_valid.set(false) - }; - } - ) + LorentzProductInput(constraint=cst) } } }, -- 2.34.1 From c2e3c64d4a3750c03fceb789bc8111cd71a67d5c Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Wed, 30 Oct 2024 00:16:34 -0700 Subject: [PATCH 11/24] Remove debug log from Lorentz product input --- app-proto/src/outline.rs | 5 ----- 1 file changed, 5 deletions(-) diff --git a/app-proto/src/outline.rs b/app-proto/src/outline.rs index fcf983f..6497fdd 100644 --- a/app-proto/src/outline.rs +++ b/app-proto/src/outline.rs @@ -1,7 +1,6 @@ use itertools::Itertools; use sycamore::{prelude::*, web::tags::div}; use web_sys::{Element, Event, HtmlInputElement, KeyboardEvent, MouseEvent, wasm_bindgen::JsCast}; -use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */ use crate::{AppState, assembly::Constraint}; @@ -18,10 +17,6 @@ fn LorentzProductInput(constraint: Constraint) -> View { Ok(rep) => batch(|| { constraint.rep.set(rep); constraint.rep_valid.set(true); - console::log_2( - &JsValue::from("Constraint rep parsed to"), - &JsValue::from(rep) - ); }), Err(_) => constraint.rep_valid.set(false) }; -- 2.34.1 From 9e31037e17d1b3f21ac95b8fa7407487a257e18f Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Wed, 30 Oct 2024 00:19:44 -0700 Subject: [PATCH 12/24] Spread web-sys imports over multiple lines --- app-proto/src/outline.rs | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/app-proto/src/outline.rs b/app-proto/src/outline.rs index 6497fdd..8edbe07 100644 --- a/app-proto/src/outline.rs +++ b/app-proto/src/outline.rs @@ -1,6 +1,13 @@ use itertools::Itertools; use sycamore::{prelude::*, web::tags::div}; -use web_sys::{Element, Event, HtmlInputElement, KeyboardEvent, MouseEvent, wasm_bindgen::JsCast}; +use web_sys::{ + Element, + Event, + HtmlInputElement, + KeyboardEvent, + MouseEvent, + wasm_bindgen::JsCast +}; use crate::{AppState, assembly::Constraint}; -- 2.34.1 From 9c191ae586369e1217b9c809b52ad7f7253c8cf4 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Wed, 30 Oct 2024 00:27:16 -0700 Subject: [PATCH 13/24] Polish log messages --- app-proto/src/add_remove.rs | 4 ++-- app-proto/src/assembly.rs | 2 +- app-proto/src/engine.rs | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/app-proto/src/add_remove.rs b/app-proto/src/add_remove.rs index 5435418..92ae4be 100644 --- a/app-proto/src/add_remove.rs +++ b/app-proto/src/add_remove.rs @@ -228,7 +228,7 @@ pub fn AddRemove() -> View { /* DEBUG */ // print updated constraint list - console::log_1(&JsValue::from("constraints:")); + console::log_1(&JsValue::from("Constraints:")); state.assembly.constraints.with(|csts| { for (_, cst) in csts.into_iter() { console::log_5( @@ -246,7 +246,7 @@ pub fn AddRemove() -> View { create_effect(move || { rep.track(); console::log_2( - &JsValue::from("Constraint rep updated to"), + &JsValue::from("Lorentz product updated to"), &JsValue::from(rep.get_untracked()) ); if active.get() { diff --git a/app-proto/src/assembly.rs b/app-proto/src/assembly.rs index 62f5405..0970932 100644 --- a/app-proto/src/assembly.rs +++ b/app-proto/src/assembly.rs @@ -146,7 +146,7 @@ impl Assembly { /* DEBUG */ // log the initial configuration matrix - console::log_1(&JsValue::from("old configuration:")); + 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() { diff --git a/app-proto/src/engine.rs b/app-proto/src/engine.rs index 2971750..2978a9a 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -60,7 +60,7 @@ impl PartialMatrix { 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); + let ent_str = format!(" {} {} {}", ent.index.0, ent.index.1, ent.value); console::log_1(&JsValue::from(ent_str.as_str())); } } -- 2.34.1 From 933f05661d7c69595791d61c79e0651ef4c2d2e7 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sun, 10 Nov 2024 16:31:29 -0800 Subject: [PATCH 14/24] Only compile `engine::point` when it's used This function will eventually be used in the application, but right now it's only used in tests. --- app-proto/src/engine.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/app-proto/src/engine.rs b/app-proto/src/engine.rs index 2978a9a..ca05646 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -4,6 +4,7 @@ use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */ // --- elements --- +#[cfg(test)] 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)]) } -- 2.34.1 From da008fd09075e27e99e2a2d30a8826d26f701434 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sun, 10 Nov 2024 19:24:26 -0800 Subject: [PATCH 15/24] Write out `representation` in `Element` structure --- app-proto/src/add_remove.rs | 28 ++++++++++++++-------------- app-proto/src/assembly.rs | 8 ++++---- app-proto/src/display.rs | 4 +++- app-proto/src/outline.rs | 2 +- 4 files changed, 22 insertions(+), 20 deletions(-) diff --git a/app-proto/src/add_remove.rs b/app-proto/src/add_remove.rs index 92ae4be..6fee142 100644 --- a/app-proto/src/add_remove.rs +++ b/app-proto/src/add_remove.rs @@ -11,7 +11,7 @@ fn load_gen_assemb(assembly: &Assembly) { id: String::from("gemini_a"), 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), + representation: engine::sphere(0.5, 0.5, 0.0, 1.0), constraints: BTreeSet::default(), index: 0 } @@ -21,7 +21,7 @@ fn load_gen_assemb(assembly: &Assembly) { id: String::from("gemini_b"), 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), + representation: engine::sphere(-0.5, -0.5, 0.0, 1.0), constraints: BTreeSet::default(), index: 0 } @@ -31,7 +31,7 @@ fn load_gen_assemb(assembly: &Assembly) { id: String::from("ursa_major"), 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), + representation: engine::sphere(-0.5, 0.5, 0.0, 0.75), constraints: BTreeSet::default(), index: 0 } @@ -41,7 +41,7 @@ fn load_gen_assemb(assembly: &Assembly) { id: String::from("ursa_minor"), 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), + representation: engine::sphere(0.5, -0.5, 0.0, 0.5), constraints: BTreeSet::default(), index: 0 } @@ -51,7 +51,7 @@ fn load_gen_assemb(assembly: &Assembly) { id: String::from("moon_deimos"), 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), + representation: engine::sphere(0.0, 0.15, 1.0, 0.25), constraints: BTreeSet::default(), index: 0 } @@ -61,7 +61,7 @@ fn load_gen_assemb(assembly: &Assembly) { id: String::from("moon_phobos"), 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), + representation: engine::sphere(0.0, -0.15, -1.0, 0.25), constraints: BTreeSet::default(), index: 0 } @@ -76,7 +76,7 @@ fn load_low_curv_assemb(assembly: &Assembly) { id: "central".to_string(), 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), + representation: engine::sphere(0.0, 0.0, 0.0, 1.0), constraints: BTreeSet::default(), index: 0 } @@ -86,7 +86,7 @@ fn load_low_curv_assemb(assembly: &Assembly) { id: "assemb_plane".to_string(), 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), + representation: engine::sphere_with_offset(0.0, 0.0, 1.0, 0.0, 0.0), constraints: BTreeSet::default(), index: 0 } @@ -96,7 +96,7 @@ fn load_low_curv_assemb(assembly: &Assembly) { id: "side1".to_string(), 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), + representation: engine::sphere_with_offset(1.0, 0.0, 0.0, 1.0, 0.0), constraints: BTreeSet::default(), index: 0 } @@ -106,7 +106,7 @@ fn load_low_curv_assemb(assembly: &Assembly) { id: "side2".to_string(), 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), + representation: engine::sphere_with_offset(-0.5, a, 0.0, 1.0, 0.0), constraints: BTreeSet::default(), index: 0 } @@ -116,7 +116,7 @@ fn load_low_curv_assemb(assembly: &Assembly) { id: "side3".to_string(), 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), + representation: engine::sphere_with_offset(-0.5, -a, 0.0, 1.0, 0.0), constraints: BTreeSet::default(), index: 0 } @@ -126,7 +126,7 @@ fn load_low_curv_assemb(assembly: &Assembly) { id: "corner1".to_string(), 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), + representation: engine::sphere(-4.0/3.0, 0.0, 0.0, 1.0/3.0), constraints: BTreeSet::default(), index: 0 } @@ -136,7 +136,7 @@ fn load_low_curv_assemb(assembly: &Assembly) { id: "corner2".to_string(), 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), + representation: engine::sphere(2.0/3.0, -4.0/3.0 * a, 0.0, 1.0/3.0), constraints: BTreeSet::default(), index: 0 } @@ -146,7 +146,7 @@ fn load_low_curv_assemb(assembly: &Assembly) { id: String::from("corner3"), 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), + representation: engine::sphere(2.0/3.0, 4.0/3.0 * a, 0.0, 1.0/3.0), constraints: BTreeSet::default(), index: 0 } diff --git a/app-proto/src/assembly.rs b/app-proto/src/assembly.rs index 0970932..27774a1 100644 --- a/app-proto/src/assembly.rs +++ b/app-proto/src/assembly.rs @@ -12,7 +12,7 @@ pub struct Element { pub id: String, pub label: String, pub color: [f32; 3], - pub rep: DVector, + pub representation: DVector, pub constraints: BTreeSet, // internal properties, not reflected in any view @@ -86,7 +86,7 @@ impl Assembly { id: id, 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]), + representation: DVector::::from_column_slice(&[0.0, 0.0, 0.0, 0.5, -0.5]), constraints: BTreeSet::default(), index: 0 } @@ -133,7 +133,7 @@ impl Assembly { 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); + guess_to_be.set_column(index, &elt.representation); } (gram_to_be, guess_to_be) @@ -177,7 +177,7 @@ impl Assembly { // read out the solution self.elements.update(|elts| { for (_, elt) in elts.iter_mut() { - elt.rep.set_column(0, &config.column(elt.index)); + elt.representation.set_column(0, &config.column(elt.index)); } }); } diff --git a/app-proto/src/display.rs b/app-proto/src/display.rs index c32b470..79199ec 100644 --- a/app-proto/src/display.rs +++ b/app-proto/src/display.rs @@ -297,7 +297,9 @@ pub fn Display() -> View { // get the assembly let elements = state.assembly.elements.get_clone(); let element_iter = (&elements).into_iter(); - let reps_world: Vec<_> = element_iter.clone().map(|(_, elt)| &assembly_to_world * &elt.rep).collect(); + let reps_world: Vec<_> = element_iter.clone().map( + |(_, elt)| &assembly_to_world * &elt.representation + ).collect(); let colors: Vec<_> = element_iter.clone().map(|(key, elt)| if state.selection.with(|sel| sel.contains(&key)) { elt.color.map(|ch| 0.2 + 0.8*ch) diff --git a/app-proto/src/outline.rs b/app-proto/src/outline.rs index 8edbe07..9ebf34d 100644 --- a/app-proto/src/outline.rs +++ b/app-proto/src/outline.rs @@ -72,7 +72,7 @@ pub fn Outline() -> View { } }); let label = elt.label.clone(); - let rep_components = elt.rep.iter().map(|u| { + let rep_components = elt.representation.iter().map(|u| { let u_coord = u.to_string().replace("-", "\u{2212}"); View::from(div().children(u_coord)) }).collect::>(); -- 2.34.1 From ed1890bffc9886a79f7c9eea7acb268f660bf5df Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sun, 10 Nov 2024 19:36:40 -0800 Subject: [PATCH 16/24] Improve naming of constraint subjects --- app-proto/src/add_remove.rs | 12 ++++++------ app-proto/src/assembly.rs | 14 +++++++------- app-proto/src/outline.rs | 6 +++--- 3 files changed, 16 insertions(+), 16 deletions(-) diff --git a/app-proto/src/add_remove.rs b/app-proto/src/add_remove.rs index 6fee142..df91d3c 100644 --- a/app-proto/src/add_remove.rs +++ b/app-proto/src/add_remove.rs @@ -208,16 +208,16 @@ pub fn AddRemove() -> View { }, on:click=|_| { let state = use_context::(); - let args = state.selection.with( + let subjects = state.selection.with( |sel| { - let arg_vec: Vec<_> = sel.into_iter().collect(); - (arg_vec[0].clone(), arg_vec[1].clone()) + let subject_vec: Vec<_> = sel.into_iter().collect(); + (subject_vec[0].clone(), subject_vec[1].clone()) } ); let rep = create_signal(0.0); let active = create_signal(true); state.assembly.insert_constraint(Constraint { - args: args, + subjects: subjects, rep: rep, rep_text: create_signal(String::new()), rep_valid: create_signal(false), @@ -233,8 +233,8 @@ pub fn AddRemove() -> View { for (_, cst) in csts.into_iter() { console::log_5( &JsValue::from(" "), - &JsValue::from(cst.args.0), - &JsValue::from(cst.args.1), + &JsValue::from(cst.subjects.0), + &JsValue::from(cst.subjects.1), &JsValue::from(":"), &JsValue::from(cst.rep.get_untracked()) ); diff --git a/app-proto/src/assembly.rs b/app-proto/src/assembly.rs index 27774a1..8217922 100644 --- a/app-proto/src/assembly.rs +++ b/app-proto/src/assembly.rs @@ -21,7 +21,7 @@ pub struct Element { #[derive(Clone)] pub struct Constraint { - pub args: (usize, usize), + pub subjects: (usize, usize), pub rep: Signal, pub rep_text: Signal, pub rep_valid: Signal, @@ -94,11 +94,11 @@ impl Assembly { } pub fn insert_constraint(&self, constraint: Constraint) { - let args = constraint.args; + let subjects = constraint.subjects; let key = self.constraints.update(|csts| csts.insert(constraint)); self.elements.update(|elts| { - elts[args.0].constraints.insert(key); - elts[args.1].constraints.insert(key); + elts[subjects.0].constraints.insert(key); + elts[subjects.1].constraints.insert(key); }); } @@ -119,9 +119,9 @@ impl Assembly { self.constraints.with_untracked(|csts| { for (_, cst) in csts { if cst.active.get_untracked() && cst.rep_valid.get_untracked() { - let args = cst.args; - let row = elts[args.0].index; - let col = elts[args.1].index; + let subjects = cst.subjects; + let row = elts[subjects.0].index; + let col = elts[subjects.1].index; gram_to_be.push_sym(row, col, cst.rep.get_untracked()); } } diff --git a/app-proto/src/outline.rs b/app-proto/src/outline.rs index 9ebf34d..98b422f 100644 --- a/app-proto/src/outline.rs +++ b/app-proto/src/outline.rs @@ -154,10 +154,10 @@ pub fn Outline() -> View { let c_state = use_context::(); let assembly = &c_state.assembly; let cst = assembly.constraints.with(|csts| csts[c_key].clone()); - let other_arg = if cst.args.0 == key { - cst.args.1 + let other_arg = if cst.subjects.0 == key { + cst.subjects.1 } else { - cst.args.0 + cst.subjects.0 }; let other_arg_label = assembly.elements.with(|elts| elts[other_arg].label.clone()); view! { -- 2.34.1 From ced001bbfe2597bebbda2bcc8659d38d7ab1d3ff Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sun, 10 Nov 2024 20:13:40 -0800 Subject: [PATCH 17/24] Alias the types of element and constraint keys This will make it easier to change the key types if we change how we store and access elements and constraints. --- app-proto/src/assembly.rs | 10 +++++++--- app-proto/src/main.rs | 4 ++-- app-proto/src/outline.rs | 2 +- 3 files changed, 10 insertions(+), 6 deletions(-) diff --git a/app-proto/src/assembly.rs b/app-proto/src/assembly.rs index 8217922..6c4aeb7 100644 --- a/app-proto/src/assembly.rs +++ b/app-proto/src/assembly.rs @@ -7,13 +7,17 @@ use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */ use crate::engine::{realize_gram, PartialMatrix}; +// the types of the keys we use to access an assembly's elements and constraints +pub type ElementKey = usize; +pub type ConstraintKey = usize; + #[derive(Clone, PartialEq)] pub struct Element { pub id: String, pub label: String, pub color: [f32; 3], pub representation: DVector, - pub constraints: BTreeSet, + pub constraints: BTreeSet, // internal properties, not reflected in any view pub index: usize @@ -21,7 +25,7 @@ pub struct Element { #[derive(Clone)] pub struct Constraint { - pub subjects: (usize, usize), + pub subjects: (ElementKey, ElementKey), pub rep: Signal, pub rep_text: Signal, pub rep_valid: Signal, @@ -36,7 +40,7 @@ pub struct Assembly { pub constraints: Signal>, // indexing - pub elements_by_id: Signal> + pub elements_by_id: Signal> } impl Assembly { diff --git a/app-proto/src/main.rs b/app-proto/src/main.rs index 2c71a83..897f9d4 100644 --- a/app-proto/src/main.rs +++ b/app-proto/src/main.rs @@ -8,14 +8,14 @@ use rustc_hash::FxHashSet; use sycamore::prelude::*; use add_remove::AddRemove; -use assembly::Assembly; +use assembly::{Assembly, ElementKey}; use display::Display; use outline::Outline; #[derive(Clone)] struct AppState { assembly: Assembly, - selection: Signal> + selection: Signal> } impl AppState { diff --git a/app-proto/src/outline.rs b/app-proto/src/outline.rs index 98b422f..0a8816d 100644 --- a/app-proto/src/outline.rs +++ b/app-proto/src/outline.rs @@ -150,7 +150,7 @@ pub fn Outline() -> View { ul(class="constraints") { Keyed( list=elt.constraints.into_iter().collect::>(), - view=move |c_key: usize| { + view=move |c_key| { let c_state = use_context::(); let assembly = &c_state.assembly; let cst = assembly.constraints.with(|csts| csts[c_key].clone()); -- 2.34.1 From b8ca1139d5108a4d7873d4d44060f2daf403d472 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sun, 10 Nov 2024 23:22:30 -0800 Subject: [PATCH 18/24] Explain what the `Element::index` field holds Also, remind us to make the field private when that becomes possible. --- app-proto/src/assembly.rs | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/app-proto/src/assembly.rs b/app-proto/src/assembly.rs index 6c4aeb7..da54210 100644 --- a/app-proto/src/assembly.rs +++ b/app-proto/src/assembly.rs @@ -19,7 +19,11 @@ pub struct Element { pub representation: DVector, pub constraints: BTreeSet, - // internal properties, not reflected in any view + // the configuration matrix column index that was assigned to this element + // last time the assembly was realized + /* TO DO */ + // this is public, as a kludge, because `Element` doesn't have a constructor + // yet. it should be made private as soon as the constructor is written pub index: usize } -- 2.34.1 From a170492e3dc1a3556062f6174551195787ce014c Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sun, 10 Nov 2024 23:47:34 -0800 Subject: [PATCH 19/24] Add engine conventions to inversive coordinates notes --- notes/inversive.md | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/notes/inversive.md b/notes/inversive.md index 5ee6329..4ca8091 100644 --- a/notes/inversive.md +++ b/notes/inversive.md @@ -41,3 +41,23 @@ I will have to work out formulas for the Euclidean distance between two entities In this vein, it seems as though if J1 and J2 are the reps of two points, then Q(J1,J2) = d^2/2. So then the sphere centered at J1 through J2 is (J1-(2Q(J1,J2),0,0,0,0))/sqrt(2Q(J1,J2)). Ugh has a sqrt in it. Similarly for sphere centered at J3 through J2, (J3-(2Q(J3,J2),0000))/sqrt(2Q(J3,J2)). J1,J2,J3 are collinear if these spheres are tangent, i.e. if those vectors have Q-inner-product 1, which is to say Q(J1,J3) - Q(J1,J2) - Q(J3,J2) = 2sqrt(Q(J1,J2)Q(J2,J3)). But maybe that's not the simplest way of putting it. After all, we can just say that the cross-product of the two differences is 0; that has no square roots in it. One conceivable way to canonicalize lines is to use the *perpendicular* plane that goes through the origin, that's uniquely defined, and anyway just amounts to I = (0,0,d) where d is the ordinary direction vector of the line; and a point J in that plane that the line goes through, which just amounts to J=(r^2,1,E) with Q(I,J) = 0, i.e. E\dot d = 0. It's also the point on the line closest to the origin. The reason that we don't usually use that point as the companion to the direction vector is that the resulting set of six coordinates is not homogeneous. But here that's not an issue, since we have our standard point coordinates and plane coordinates; and for a plane through the origin, only two of the direction coordinates are really free, and then we have the one dot-product relation, so only two of the point coordinates are really free, giving us the correct dimensionality of 4 for the set of lines. So in some sense this says that we could take naively as coordinates for a line the projection of the unit direction vector to the xy plane and the projection of the line's closest point to the origin to the xy plane. That doesn't seem to have any weird gimbal locks or discontinuities or anything. And with these coordinates, you can test if the point E=x,y,z is on the line (dx,dy,cx,cy) by extending (dx,dy) to d via dz = sqrt(1-dx^2 - dy^2), extending (cx,cy) to c by determining cz via d\dot c = 0, and then checking if d\cross(E-c) = 0. And you can see if two lines are parallel just by checking if they have the same direction vector, and if not, you can see if they are coplanar by projecting both of their closest points perpendicularly onto the line in the direction of the cross product of their directions, and if the projections match they are coplanar. + +#### Engine Conventions + +The coordinate conventions used in the engine are different from the ones used in these notes. Marking the engine vectors and coordinates with $'$, we have +$$I' = (x', y', z', b', c'),$$ +where +$$\begin{align*} +x' & = x & b' & = b/2 \\ +y' & = y & c' & = c/2. \\ +z' & = z +\end{align*}$$ +The engine uses the quadratic form $Q' = -Q$, which is expressed in engine coordinates as +$$Q'(I'_1, I'_2) = x'_1 x'_2 + y'_1 y'_2 + z'_1 z'_2 - 2(b'_1c'_2 + c'_1 b'_2).$$ +In the `engine` module, the matrix of $Q'$ is encoded in the lazy static variable `Q`. + +In the engine's coordinate conventions, a sphere with radius $r > 0$ centered on $P = (P_x, P_y, P_z)$ is represented by the vector +$$I'_s = \left(\frac{P_x}{r}, \frac{P_y}{r}, \frac{P_z}{r}, \frac1{2r}, \frac{\|P\|^2 - r^2}{2r}\right),$$ +which has the normalization $Q'(I'_s, I'_s) = 1$. The point $P$ is represented by the vector +$$I'_P = \left(P_x, P_y, P_z, \frac{1}{2}, \frac{\|P\|^2}{2}\right).$$ +In the `engine` module, these formulas are encoded in the `sphere` and `point` functions. \ No newline at end of file -- 2.34.1 From abb9d353358ed21a9112233ef8d5324af56f92bf Mon Sep 17 00:00:00 2001 From: Vectornaut Date: Mon, 11 Nov 2024 07:56:27 +0000 Subject: [PATCH 20/24] Correct `align` environment in notes --- notes/inversive.md | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/notes/inversive.md b/notes/inversive.md index 4ca8091..933eb35 100644 --- a/notes/inversive.md +++ b/notes/inversive.md @@ -47,11 +47,13 @@ One conceivable way to canonicalize lines is to use the *perpendicular* plane th The coordinate conventions used in the engine are different from the ones used in these notes. Marking the engine vectors and coordinates with $'$, we have $$I' = (x', y', z', b', c'),$$ where -$$\begin{align*} +$$ +\begin{align*} x' & = x & b' & = b/2 \\ y' & = y & c' & = c/2. \\ z' & = z -\end{align*}$$ +\end{align*} +$$ The engine uses the quadratic form $Q' = -Q$, which is expressed in engine coordinates as $$Q'(I'_1, I'_2) = x'_1 x'_2 + y'_1 y'_2 + z'_1 z'_2 - 2(b'_1c'_2 + c'_1 b'_2).$$ In the `engine` module, the matrix of $Q'$ is encoded in the lazy static variable `Q`. -- 2.34.1 From a4ec52a4e7b6ac1f797b14e7a38bb659a588d6c5 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Mon, 11 Nov 2024 00:04:48 -0800 Subject: [PATCH 21/24] Alias the type of an element's color --- app-proto/src/assembly.rs | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/app-proto/src/assembly.rs b/app-proto/src/assembly.rs index da54210..485448e 100644 --- a/app-proto/src/assembly.rs +++ b/app-proto/src/assembly.rs @@ -11,11 +11,13 @@ use crate::engine::{realize_gram, PartialMatrix}; pub type ElementKey = usize; pub type ConstraintKey = usize; +pub type ElementColor = [f32; 3]; + #[derive(Clone, PartialEq)] pub struct Element { pub id: String, pub label: String, - pub color: [f32; 3], + pub color: ElementColor, pub representation: DVector, pub constraints: BTreeSet, -- 2.34.1 From 5839882ed7a1cd391031978ac14d9d58bcad8173 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Mon, 11 Nov 2024 00:26:53 -0800 Subject: [PATCH 22/24] Name constraint Lorentz product more descriptively --- app-proto/src/add_remove.rs | 14 +++++++------- app-proto/src/assembly.rs | 10 +++++----- app-proto/src/outline.rs | 10 +++++----- 3 files changed, 17 insertions(+), 17 deletions(-) diff --git a/app-proto/src/add_remove.rs b/app-proto/src/add_remove.rs index df91d3c..19b4b8d 100644 --- a/app-proto/src/add_remove.rs +++ b/app-proto/src/add_remove.rs @@ -214,13 +214,13 @@ pub fn AddRemove() -> View { (subject_vec[0].clone(), subject_vec[1].clone()) } ); - let rep = create_signal(0.0); + let lorentz_prod = create_signal(0.0); let active = create_signal(true); state.assembly.insert_constraint(Constraint { subjects: subjects, - rep: rep, - rep_text: create_signal(String::new()), - rep_valid: create_signal(false), + lorentz_prod: lorentz_prod, + lorentz_prod_text: create_signal(String::new()), + lorentz_prod_valid: create_signal(false), active: active, }); state.assembly.realize(); @@ -236,7 +236,7 @@ pub fn AddRemove() -> View { &JsValue::from(cst.subjects.0), &JsValue::from(cst.subjects.1), &JsValue::from(":"), - &JsValue::from(cst.rep.get_untracked()) + &JsValue::from(cst.lorentz_prod.get_untracked()) ); } }); @@ -244,10 +244,10 @@ pub fn AddRemove() -> View { // update the realization when the constraint activated, or // edited while active create_effect(move || { - rep.track(); + lorentz_prod.track(); console::log_2( &JsValue::from("Lorentz product updated to"), - &JsValue::from(rep.get_untracked()) + &JsValue::from(lorentz_prod.get_untracked()) ); if active.get() { state.assembly.realize(); diff --git a/app-proto/src/assembly.rs b/app-proto/src/assembly.rs index 485448e..0cdf61b 100644 --- a/app-proto/src/assembly.rs +++ b/app-proto/src/assembly.rs @@ -32,9 +32,9 @@ pub struct Element { #[derive(Clone)] pub struct Constraint { pub subjects: (ElementKey, ElementKey), - pub rep: Signal, - pub rep_text: Signal, - pub rep_valid: Signal, + pub lorentz_prod: Signal, + pub lorentz_prod_text: Signal, + pub lorentz_prod_valid: Signal, pub active: Signal } @@ -128,11 +128,11 @@ impl Assembly { let mut gram_to_be = PartialMatrix::new(); self.constraints.with_untracked(|csts| { for (_, cst) in csts { - if cst.active.get_untracked() && cst.rep_valid.get_untracked() { + if cst.active.get_untracked() && cst.lorentz_prod_valid.get_untracked() { let subjects = cst.subjects; let row = elts[subjects.0].index; let col = elts[subjects.1].index; - gram_to_be.push_sym(row, col, cst.rep.get_untracked()); + gram_to_be.push_sym(row, col, cst.lorentz_prod.get_untracked()); } } }); diff --git a/app-proto/src/outline.rs b/app-proto/src/outline.rs index 0a8816d..f7c975c 100644 --- a/app-proto/src/outline.rs +++ b/app-proto/src/outline.rs @@ -17,15 +17,15 @@ fn LorentzProductInput(constraint: Constraint) -> View { view! { input( r#type="text", - bind:value=constraint.rep_text, + bind:value=constraint.lorentz_prod_text, on:change=move |event: Event| { let target: HtmlInputElement = event.target().unwrap().unchecked_into(); match target.value().parse::() { - Ok(rep) => batch(|| { - constraint.rep.set(rep); - constraint.rep_valid.set(true); + Ok(lorentz_prod) => batch(|| { + constraint.lorentz_prod.set(lorentz_prod); + constraint.lorentz_prod_valid.set(true); }), - Err(_) => constraint.rep_valid.set(false) + Err(_) => constraint.lorentz_prod_valid.set(false) }; } ) -- 2.34.1 From 22a93bee28f2ace712cc292cb7447851a4bf6fe6 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Mon, 11 Nov 2024 15:34:51 -0800 Subject: [PATCH 23/24] Confirm that frozen entries are frozen exactly --- app-proto/src/engine.rs | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/app-proto/src/engine.rs b/app-proto/src/engine.rs index ca05646..1c53ac8 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -414,6 +414,40 @@ mod tests { } } + // at the frozen indices, the optimization steps should have exact zeros, + // and the realized configuration should match the initial guess + #[test] + fn frozen_entry_test() { + let gram = { + let mut gram_to_be = PartialMatrix::new(); + for j in 0..2 { + for k in j..2 { + gram_to_be.push_sym(j, k, if (j, k) == (1, 1) { 1.0 } else { 0.0 }); + } + } + gram_to_be + }; + 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), (3, 1)]; + println!(); + let (config, success, history) = realize_gram( + &gram, guess.clone(), &frozen, + 1.0e-12, 0.5, 0.9, 1.1, 200, 110 + ); + assert_eq!(success, true); + for base_step in history.base_step.into_iter() { + for index in frozen { + assert_eq!(base_step[index], 0.0); + } + } + for index in frozen { + assert_eq!(config[index], guess[index]); + } + } + // --- process inspection examples --- // these tests are meant for human inspection, not automated use. run them -- 2.34.1 From 5332fda6e42b7e9c52d8ec48663012d4e12054dc Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Mon, 11 Nov 2024 15:41:42 -0800 Subject: [PATCH 24/24] Move new test to avoid merge conflict --- app-proto/src/engine.rs | 71 +++++++++++++++++++++-------------------- 1 file changed, 37 insertions(+), 34 deletions(-) diff --git a/app-proto/src/engine.rs b/app-proto/src/engine.rs index 1c53ac8..343b96e 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -414,40 +414,6 @@ mod tests { } } - // at the frozen indices, the optimization steps should have exact zeros, - // and the realized configuration should match the initial guess - #[test] - fn frozen_entry_test() { - let gram = { - let mut gram_to_be = PartialMatrix::new(); - for j in 0..2 { - for k in j..2 { - gram_to_be.push_sym(j, k, if (j, k) == (1, 1) { 1.0 } else { 0.0 }); - } - } - gram_to_be - }; - 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), (3, 1)]; - println!(); - let (config, success, history) = realize_gram( - &gram, guess.clone(), &frozen, - 1.0e-12, 0.5, 0.9, 1.1, 200, 110 - ); - assert_eq!(success, true); - for base_step in history.base_step.into_iter() { - for index in frozen { - assert_eq!(base_step[index], 0.0); - } - } - for index in frozen { - assert_eq!(config[index], guess[index]); - } - } - // --- process inspection examples --- // these tests are meant for human inspection, not automated use. run them @@ -534,4 +500,41 @@ mod tests { println!("{:<4} │ {}", step, scaled_loss); } } + + /* TO DO */ + // --- new test placed here to avoid merge conflict --- + + // at the frozen indices, the optimization steps should have exact zeros, + // and the realized configuration should match the initial guess + #[test] + fn frozen_entry_test() { + let gram = { + let mut gram_to_be = PartialMatrix::new(); + for j in 0..2 { + for k in j..2 { + gram_to_be.push_sym(j, k, if (j, k) == (1, 1) { 1.0 } else { 0.0 }); + } + } + gram_to_be + }; + 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), (3, 1)]; + println!(); + let (config, success, history) = realize_gram( + &gram, guess.clone(), &frozen, + 1.0e-12, 0.5, 0.9, 1.1, 200, 110 + ); + assert_eq!(success, true); + for base_step in history.base_step.into_iter() { + for index in frozen { + assert_eq!(base_step[index], 0.0); + } + } + for index in frozen { + assert_eq!(config[index], guess[index]); + } + } } \ No newline at end of file -- 2.34.1