642 lines
No EOL
23 KiB
Rust
642 lines
No EOL
23 KiB
Rust
use lazy_static::lazy_static;
|
|
use nalgebra::{Const, DMatrix, DVector, DVectorView, Dyn, SymmetricEigen};
|
|
use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */
|
|
|
|
// --- elements ---
|
|
|
|
#[cfg(feature = "dev")]
|
|
pub fn point(x: f64, y: f64, z: f64) -> DVector<f64> {
|
|
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<f64> {
|
|
let center_norm_sq = center_x * center_x + center_y * center_y + center_z * center_z;
|
|
DVector::from_column_slice(&[
|
|
center_x / radius,
|
|
center_y / radius,
|
|
center_z / radius,
|
|
0.5 / radius,
|
|
0.5 * (center_norm_sq / radius - radius)
|
|
])
|
|
}
|
|
|
|
// the sphere of curvature `curv` whose closest point to the origin has position
|
|
// `off * dir` and normal `dir`, where `dir` is a unit vector. setting the
|
|
// curvature to zero gives a plane
|
|
pub fn sphere_with_offset(dir_x: f64, dir_y: f64, dir_z: f64, off: f64, curv: f64) -> DVector<f64> {
|
|
let norm_sp = 1.0 + off * curv;
|
|
DVector::from_column_slice(&[
|
|
norm_sp * dir_x,
|
|
norm_sp * dir_y,
|
|
norm_sp * dir_z,
|
|
0.5 * curv,
|
|
off * (1.0 + 0.5 * off * curv)
|
|
])
|
|
}
|
|
|
|
// --- partial matrices ---
|
|
|
|
struct MatrixEntry {
|
|
index: (usize, usize),
|
|
value: f64
|
|
}
|
|
|
|
pub struct PartialMatrix(Vec<MatrixEntry>);
|
|
|
|
impl PartialMatrix {
|
|
pub fn new() -> PartialMatrix {
|
|
PartialMatrix(Vec::<MatrixEntry>::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<f64>) -> DMatrix<f64> {
|
|
let mut result = DMatrix::<f64>::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<f64>) -> DMatrix<f64> {
|
|
let mut result = DMatrix::<f64>::zeros(rhs.nrows(), rhs.ncols());
|
|
let PartialMatrix(entries) = self;
|
|
for ent in entries {
|
|
result[ent.index] = ent.value - rhs[ent.index];
|
|
}
|
|
result
|
|
}
|
|
}
|
|
|
|
// --- configuration subspaces ---
|
|
|
|
#[derive(Clone)]
|
|
pub struct ConfigSubspace {
|
|
assembly_dim: usize,
|
|
basis_std: Vec<DMatrix<f64>>,
|
|
basis_proj: Vec<DMatrix<f64>>
|
|
}
|
|
|
|
impl ConfigSubspace {
|
|
pub fn zero(assembly_dim: usize) -> ConfigSubspace {
|
|
ConfigSubspace {
|
|
assembly_dim: assembly_dim,
|
|
basis_proj: Vec::new(),
|
|
basis_std: Vec::new()
|
|
}
|
|
}
|
|
|
|
// approximate the kernel of a symmetric endomorphism of the configuration
|
|
// space for `assembly_dim` elements. we consider an eigenvector to be part
|
|
// of the kernel if its eigenvalue is smaller than the constant `THRESHOLD`
|
|
fn symmetric_kernel(a: DMatrix<f64>, proj_to_std: DMatrix<f64>, assembly_dim: usize) -> ConfigSubspace {
|
|
// find a basis for the kernel, expressed in the standard coordinates
|
|
const ELEMENT_DIM: usize = 5;
|
|
const THRESHOLD: f64 = 1.0e-4;
|
|
let eig = SymmetricEigen::new(a);
|
|
let eig_vecs = eig.eigenvectors.column_iter();
|
|
let eig_pairs = eig.eigenvalues.iter().zip(eig_vecs);
|
|
let basis_std = DMatrix::from_columns(
|
|
eig_pairs.filter_map(
|
|
|(λ, v)| (λ.abs() < THRESHOLD).then_some(v)
|
|
).collect::<Vec<_>>().as_slice()
|
|
);
|
|
|
|
/* DEBUG */
|
|
// print the eigenvalues
|
|
#[cfg(all(target_family = "wasm", target_os = "unknown"))]
|
|
console::log_1(&JsValue::from(
|
|
format!("Eigenvalues used to find kernel:{}", eig.eigenvalues)
|
|
));
|
|
|
|
// express the basis in the projection coordinates
|
|
let basis_proj = proj_to_std.clone().qr().solve(&basis_std).unwrap();
|
|
|
|
// orthonormalize the basis with respect to the projection inner product
|
|
let basis_proj_orth = basis_proj.qr().q();
|
|
let basis_std_orth = proj_to_std * &basis_proj_orth;
|
|
|
|
// print the projection basis in projection coordinates
|
|
#[cfg(all(target_family = "wasm", target_os = "unknown"))]
|
|
console::log_1(&JsValue::from(
|
|
format!("Basis in projection coordinates:{}", basis_proj_orth)
|
|
));
|
|
|
|
ConfigSubspace {
|
|
assembly_dim: assembly_dim,
|
|
basis_std: basis_std_orth.column_iter().map(
|
|
|v| Into::<DMatrix<f64>>::into(
|
|
v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim))
|
|
)
|
|
).collect(),
|
|
basis_proj: basis_proj_orth.column_iter().map(
|
|
|v| Into::<DMatrix<f64>>::into(
|
|
v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim))
|
|
)
|
|
).collect()
|
|
}
|
|
}
|
|
|
|
pub fn dim(&self) -> usize {
|
|
self.basis_std.len()
|
|
}
|
|
|
|
pub fn assembly_dim(&self) -> usize {
|
|
self.assembly_dim
|
|
}
|
|
|
|
// find the projection onto this subspace of the motion where the element
|
|
// with the given column index has velocity `v`. the velocity is given in
|
|
// projection coordinates, and the projection is done with respect to the
|
|
// projection inner product
|
|
pub fn proj(&self, v: &DVectorView<f64>, column_index: usize) -> DMatrix<f64> {
|
|
if self.dim() == 0 {
|
|
const ELEMENT_DIM: usize = 5;
|
|
DMatrix::zeros(ELEMENT_DIM, self.assembly_dim)
|
|
} else {
|
|
self.basis_proj.iter().zip(self.basis_std.iter()).map(
|
|
|(b_proj, b_std)| b_proj.column(column_index).dot(&v) * b_std
|
|
).sum()
|
|
}
|
|
}
|
|
}
|
|
|
|
// --- descent history ---
|
|
|
|
pub struct DescentHistory {
|
|
pub config: Vec<DMatrix<f64>>,
|
|
pub scaled_loss: Vec<f64>,
|
|
pub neg_grad: Vec<DMatrix<f64>>,
|
|
pub min_eigval: Vec<f64>,
|
|
pub base_step: Vec<DMatrix<f64>>,
|
|
pub backoff_steps: Vec<i32>
|
|
}
|
|
|
|
impl DescentHistory {
|
|
fn new() -> DescentHistory {
|
|
DescentHistory {
|
|
config: Vec::<DMatrix<f64>>::new(),
|
|
scaled_loss: Vec::<f64>::new(),
|
|
neg_grad: Vec::<DMatrix<f64>>::new(),
|
|
min_eigval: Vec::<f64>::new(),
|
|
base_step: Vec::<DMatrix<f64>>::new(),
|
|
backoff_steps: Vec::<i32>::new(),
|
|
}
|
|
}
|
|
}
|
|
|
|
// --- gram matrix realization ---
|
|
|
|
// the Lorentz form
|
|
lazy_static! {
|
|
pub static ref Q: DMatrix<f64> = 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<f64>,
|
|
err_proj: DMatrix<f64>,
|
|
loss: f64
|
|
}
|
|
|
|
impl SearchState {
|
|
fn from_config(gram: &PartialMatrix, config: DMatrix<f64>) -> 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<f64> {
|
|
let mut result = DMatrix::<f64>::zeros(nrows, ncols);
|
|
result[index] = 1.0;
|
|
result
|
|
}
|
|
|
|
// given a normalized vector `v` representing an element, build a basis for the
|
|
// element's linear configuration space consisting of:
|
|
// - the unit translation motions of the element
|
|
// - the unit shrinking motion of the element, if it's a sphere
|
|
// - one or two vectors whose coefficients vanish on the tangent space of the
|
|
// normalization variety
|
|
pub fn local_unif_to_std(v: DVectorView<f64>) -> DMatrix<f64> {
|
|
const ELEMENT_DIM: usize = 5;
|
|
let curv = 2.0*v[3];
|
|
if v.dot(&(&*Q * v)) < 0.5 {
|
|
// `v` represents a point. the normalization condition says that the
|
|
// curvature component of `v` is 1/2
|
|
DMatrix::from_column_slice(ELEMENT_DIM, ELEMENT_DIM, &[
|
|
curv, 0.0, 0.0, 0.0, v[0],
|
|
0.0, curv, 0.0, 0.0, v[1],
|
|
0.0, 0.0, curv, 0.0, v[2],
|
|
v[0], v[1], v[2], v[3], v[4],
|
|
0.0, 0.0, 0.0, 0.0, 1.0
|
|
])
|
|
} else {
|
|
// `v` represents a sphere. the normalization condition says that the
|
|
// Lorentz product of `v` with itself is 1
|
|
DMatrix::from_column_slice(ELEMENT_DIM, ELEMENT_DIM, &[
|
|
curv, 0.0, 0.0, 0.0, v[0],
|
|
0.0, curv, 0.0, 0.0, v[1],
|
|
0.0, 0.0, curv, 0.0, v[2],
|
|
curv*v[0], curv*v[1], curv*v[2], curv*v[3], curv*v[4] + 1.0,
|
|
v[0], v[1], v[2], v[3], v[4]
|
|
])
|
|
}
|
|
}
|
|
|
|
// use backtracking line search to find a better configuration
|
|
fn seek_better_config(
|
|
gram: &PartialMatrix,
|
|
state: &SearchState,
|
|
base_step: &DMatrix<f64>,
|
|
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<f64>,
|
|
frozen: &[(usize, usize)],
|
|
scaled_tol: f64,
|
|
min_efficiency: f64,
|
|
backoff: f64,
|
|
reg_scale: f64,
|
|
max_descent_steps: i32,
|
|
max_backoff_steps: i32
|
|
) -> (DMatrix<f64>, ConfigSubspace, 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<usize> = 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);
|
|
let mut hess = DMatrix::zeros(element_dim, assembly_dim);
|
|
for _ in 0..max_descent_steps {
|
|
// 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::<DVector<f64>>::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>));
|
|
}
|
|
}
|
|
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;
|
|
}
|
|
|
|
// 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; }
|
|
|
|
// 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.clone().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, ConfigSubspace::zero(assembly_dim), false, history)
|
|
};
|
|
}
|
|
let success = state.loss < tol;
|
|
let tangent = if success {
|
|
// express the uniform basis in the standard basis
|
|
let mut unif_to_std = DMatrix::<f64>::zeros(total_dim, total_dim);
|
|
for n in 0..assembly_dim {
|
|
let block_start = element_dim * n;
|
|
unif_to_std
|
|
.view_mut((block_start, block_start), (element_dim, element_dim))
|
|
.copy_from(&local_unif_to_std(state.config.column(n)));
|
|
}
|
|
|
|
// find the kernel of the Hessian. give it the uniform inner product
|
|
ConfigSubspace::symmetric_kernel(hess, unif_to_std, assembly_dim)
|
|
} else {
|
|
ConfigSubspace::zero(assembly_dim)
|
|
};
|
|
(state.config, tangent, success, history)
|
|
}
|
|
|
|
// --- tests ---
|
|
|
|
// 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/
|
|
//
|
|
#[cfg(feature = "dev")]
|
|
pub mod irisawa {
|
|
use std::{array, f64::consts::PI};
|
|
|
|
use super::*;
|
|
|
|
pub fn realize_irisawa_hexlet(scaled_tol: f64) -> (DMatrix<f64>, ConfigSubspace, bool, DescentHistory) {
|
|
let gram = {
|
|
let mut gram_to_be = PartialMatrix::new();
|
|
for s in 0..9 {
|
|
// each sphere is represented by a spacelike vector
|
|
gram_to_be.push_sym(s, s, 1.0);
|
|
|
|
// the circumscribing sphere is tangent to all of the other
|
|
// spheres, with matching orientation
|
|
if s > 0 {
|
|
gram_to_be.push_sym(0, s, 1.0);
|
|
}
|
|
|
|
if s > 2 {
|
|
// each chain sphere is tangent to the "sun" and "moon"
|
|
// spheres, with opposing orientation
|
|
for n in 1..3 {
|
|
gram_to_be.push_sym(s, n, -1.0);
|
|
}
|
|
|
|
// each chain sphere is tangent to the next chain sphere,
|
|
// with opposing orientation
|
|
let s_next = 3 + (s-2) % 6;
|
|
gram_to_be.push_sym(s, s_next, -1.0);
|
|
}
|
|
}
|
|
gram_to_be
|
|
};
|
|
|
|
let guess = DMatrix::from_columns(
|
|
[
|
|
sphere(0.0, 0.0, 0.0, 15.0),
|
|
sphere(0.0, 0.0, -9.0, 5.0),
|
|
sphere(0.0, 0.0, 11.0, 3.0)
|
|
].into_iter().chain(
|
|
(1..=6).map(
|
|
|k| {
|
|
let ang = (k as f64) * PI/3.0;
|
|
sphere(9.0 * ang.cos(), 9.0 * ang.sin(), 0.0, 2.5)
|
|
}
|
|
)
|
|
).collect::<Vec<_>>().as_slice()
|
|
);
|
|
|
|
// the frozen entries fix the radii of the circumscribing sphere, the
|
|
// "sun" and "moon" spheres, and one of the chain spheres
|
|
let frozen: [(usize, usize); 4] = array::from_fn(|k| (3, k));
|
|
|
|
realize_gram(
|
|
&gram, guess, &frozen,
|
|
scaled_tol, 0.5, 0.9, 1.1, 200, 110
|
|
)
|
|
}
|
|
}
|
|
|
|
#[cfg(test)]
|
|
mod tests {
|
|
use super::{*, irisawa::realize_irisawa_hexlet};
|
|
|
|
#[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::<f64>::from_row_slice(2, 3, &[
|
|
1.0, 2.0, 3.0,
|
|
4.0, 5.0, 6.0
|
|
]);
|
|
let expected_result = DMatrix::<f64>::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::<MatrixEntry>::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);
|
|
}
|
|
|
|
#[test]
|
|
fn irisawa_hexlet_test() {
|
|
// solve Irisawa's problem
|
|
const SCALED_TOL: f64 = 1.0e-12;
|
|
let (config, _, _, _) = realize_irisawa_hexlet(SCALED_TOL);
|
|
|
|
// check against Irisawa's solution
|
|
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);
|
|
}
|
|
}
|
|
|
|
#[test]
|
|
fn tangent_test() {
|
|
const SCALED_TOL: f64 = 1.0e-12;
|
|
const ELEMENT_DIM: usize = 5;
|
|
const ASSEMBLY_DIM: usize = 3;
|
|
let gram = {
|
|
let mut gram_to_be = PartialMatrix::new();
|
|
for j in 0..3 {
|
|
for k in j..3 {
|
|
gram_to_be.push_sym(j, k, if j == k { 1.0 } else { -1.0 });
|
|
}
|
|
}
|
|
gram_to_be
|
|
};
|
|
let guess = DMatrix::from_columns(&[
|
|
sphere(0.0, 0.0, 0.0, -2.0),
|
|
sphere(0.0, 0.0, 1.0, 1.0),
|
|
sphere(0.0, 0.0, -1.0, 1.0)
|
|
]);
|
|
let frozen: [_; 5] = std::array::from_fn(|k| (k, 0));
|
|
let (config, tangent, success, history) = realize_gram(
|
|
&gram, guess.clone(), &frozen,
|
|
SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
|
|
);
|
|
assert_eq!(config, guess);
|
|
assert_eq!(success, true);
|
|
assert_eq!(history.scaled_loss.len(), 1);
|
|
|
|
// confirm that the tangent space has dimension five or less
|
|
let ConfigSubspace(ref tangent_basis) = tangent;
|
|
assert_eq!(tangent_basis.len(), 5);
|
|
|
|
// confirm that the tangent space contains all the motions we expect it
|
|
// to. since we've already bounded the dimension of the tangent space,
|
|
// this confirms that the tangent space is what we expect it to be
|
|
let tangent_motions = vec![
|
|
basis_matrix((0, 1), ELEMENT_DIM, ASSEMBLY_DIM),
|
|
basis_matrix((1, 1), ELEMENT_DIM, ASSEMBLY_DIM),
|
|
basis_matrix((0, 2), ELEMENT_DIM, ASSEMBLY_DIM),
|
|
basis_matrix((1, 2), ELEMENT_DIM, ASSEMBLY_DIM),
|
|
DMatrix::<f64>::from_column_slice(ELEMENT_DIM, 3, &[
|
|
0.0, 0.0, 0.0, 0.0, 0.0,
|
|
0.0, 0.0, -1.0, -0.25, -1.0,
|
|
0.0, 0.0, -1.0, 0.25, 1.0
|
|
])
|
|
];
|
|
let tol_sq = ((ELEMENT_DIM * ASSEMBLY_DIM) as f64) * SCALED_TOL * SCALED_TOL;
|
|
for motion in tangent_motions {
|
|
let motion_proj: DMatrix<_> = motion.column_iter().enumerate().map(
|
|
|(k, v)| tangent.proj(&v, k)
|
|
).sum();
|
|
assert!((motion - motion_proj).norm_squared() < tol_sq);
|
|
}
|
|
}
|
|
|
|
// at the frozen indices, the optimization steps should have exact zeros,
|
|
// and the realized configuration should match the initial guess
|
|
#[test]
|
|
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]);
|
|
}
|
|
}
|
|
} |