glen/dyna3
1
0
Fork 0
Code Issues 18 Pull Requests 1 Actions Projects Releases 2 Wiki Activity

Do symmetric_kernel in projection coordinates

Browse Source 
Instead of finding the kernel in the standard coordinates and then
expressing it in the projection coordinates, work in the projection
coordinates from the beginning by applying a change of basis to the
input matrix.
This commit is contained in:
Aaron Fenyes 2025-01-21 18:52:01 -08:00
parent 01261aed91
commit 21cefa9f8a
3 changed files with 25 additions and 25 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
app-proto/examples/kaleidocycle.rs
View File

@ -59,7 +59,7 @@ fn main() {
println!("Loss: {}\n", history.scaled_loss.last().unwrap()); println!("Loss: {}\n", history.scaled_loss.last().unwrap());
// find the kaleidocycle's twist motion // find the kaleidocycle's twist motion
let up = DVector::from_column_slice(&[0.0, 0.0, 1.0, 0.0, 0.0]); let up = DVector::from_column_slice(&[0.0, 0.0, 1.0, 0.0]);
let down = -&up; let down = -&up;
let twist_motion: DMatrix<_> = (0..N_POINTS).step_by(4).flat_map( let twist_motion: DMatrix<_> = (0..N_POINTS).step_by(4).flat_map(
|n| [ |n| [

2
app-proto/src/display.rs
View File

@ -346,7 +346,7 @@ pub fn Display() -> View {
Vector3::zeros() Vector3::zeros()
}; };
time_step * DVector::from_column_slice( time_step * DVector::from_column_slice(
&[u[0], u[1], u[2], SHRINKING_SPEED * shrink, 0.0] &[u[0], u[1], u[2], SHRINKING_SPEED * shrink]
) )
}; };
assembly_for_raf.deform( assembly_for_raf.deform(

46
app-proto/src/engine.rs
View File

@ -107,13 +107,14 @@ impl ConfigSubspace {
// space for `assembly_dim` elements. we consider an eigenvector to be part // space for `assembly_dim` elements. we consider an eigenvector to be part
// of the kernel if its eigenvalue is smaller than the constant `THRESHOLD` // 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 { 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 // find a basis for the kernel. the basis is expressed in the projection
const ELEMENT_DIM: usize = 5; // coordinates, and it's orthonormal with respect to the projection
const THRESHOLD: f64 = 1.0e-4; // inner product
let eig = SymmetricEigen::new(a); const THRESHOLD: f64 = 0.1;
let eig = SymmetricEigen::new(proj_to_std.tr_mul(&a) * &proj_to_std);
let eig_vecs = eig.eigenvectors.column_iter(); let eig_vecs = eig.eigenvectors.column_iter();
let eig_pairs = eig.eigenvalues.iter().zip(eig_vecs); let eig_pairs = eig.eigenvalues.iter().zip(eig_vecs);
let basis_std = DMatrix::from_columns( let basis_proj = DMatrix::from_columns(
eig_pairs.filter_map( eig_pairs.filter_map(
|(λ, v)| (λ.abs() < THRESHOLD).then_some(v) |(λ, v)| (λ.abs() < THRESHOLD).then_some(v)
).collect::<Vec<_>>().as_slice() ).collect::<Vec<_>>().as_slice()
@ -126,29 +127,27 @@ impl ConfigSubspace {
format!("Eigenvalues used to find kernel:{}", eig.eigenvalues) format!("Eigenvalues used to find kernel:{}", eig.eigenvalues)
)); ));
// express the basis in the projection coordinates // express the basis in the standard coordinates
let basis_proj = proj_to_std.clone().qr().solve(&basis_std).unwrap(); let basis_std = proj_to_std * &basis_proj;
// 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 // print the projection basis in projection coordinates
#[cfg(all(target_family = "wasm", target_os = "unknown"))] #[cfg(all(target_family = "wasm", target_os = "unknown"))]
console::log_1(&JsValue::from( console::log_1(&JsValue::from(
format!("Basis in projection coordinates:{}", basis_proj_orth) format!("Basis in projection coordinates:{}", basis_proj)
)); ));
const ELEMENT_DIM: usize = 5;
const UNIFORM_DIM: usize = 4;
ConfigSubspace { ConfigSubspace {
assembly_dim: assembly_dim, assembly_dim: assembly_dim,
basis_std: basis_std_orth.column_iter().map( basis_std: basis_std.column_iter().map(
|v| Into::<DMatrix<f64>>::into( |v| Into::<DMatrix<f64>>::into(
v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim)) v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim))
) )
).collect(), ).collect(),
basis_proj: basis_proj_orth.column_iter().map( basis_proj: basis_proj.column_iter().map(
|v| Into::<DMatrix<f64>>::into( |v| Into::<DMatrix<f64>>::into(
v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim)) v.reshape_generic(Dyn(UNIFORM_DIM), Dyn(assembly_dim))
) )
).collect() ).collect()
} }
@ -247,26 +246,25 @@ fn basis_matrix(index: (usize, usize), nrows: usize, ncols: usize) -> DMatrix<f6
// normalization variety // normalization variety
pub fn local_unif_to_std(v: DVectorView<f64>) -> DMatrix<f64> { pub fn local_unif_to_std(v: DVectorView<f64>) -> DMatrix<f64> {
const ELEMENT_DIM: usize = 5; const ELEMENT_DIM: usize = 5;
const UNIFORM_DIM: usize = 4;
let curv = 2.0*v[3]; let curv = 2.0*v[3];
if v.dot(&(&*Q * v)) < 0.5 { if v.dot(&(&*Q * v)) < 0.5 {
// `v` represents a point. the normalization condition says that the // `v` represents a point. the normalization condition says that the
// curvature component of `v` is 1/2 // curvature component of `v` is 1/2
DMatrix::from_column_slice(ELEMENT_DIM, ELEMENT_DIM, &[ DMatrix::from_column_slice(ELEMENT_DIM, UNIFORM_DIM, &[
curv, 0.0, 0.0, 0.0, v[0], curv, 0.0, 0.0, 0.0, v[0],
0.0, curv, 0.0, 0.0, v[1], 0.0, curv, 0.0, 0.0, v[1],
0.0, 0.0, curv, 0.0, v[2], 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 0.0, 0.0, 0.0, 0.0, 1.0
]) ])
} else { } else {
// `v` represents a sphere. the normalization condition says that the // `v` represents a sphere. the normalization condition says that the
// Lorentz product of `v` with itself is 1 // Lorentz product of `v` with itself is 1
DMatrix::from_column_slice(ELEMENT_DIM, ELEMENT_DIM, &[ DMatrix::from_column_slice(ELEMENT_DIM, UNIFORM_DIM, &[
curv, 0.0, 0.0, 0.0, v[0], curv, 0.0, 0.0, 0.0, v[0],
0.0, curv, 0.0, 0.0, v[1], 0.0, curv, 0.0, 0.0, v[1],
0.0, 0.0, curv, 0.0, v[2], 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, 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]
]) ])
} }
} }
@ -401,11 +399,13 @@ pub fn realize_gram(
let success = state.loss < tol; let success = state.loss < tol;
let tangent = if success { let tangent = if success {
// express the uniform basis in the standard basis // express the uniform basis in the standard basis
let mut unif_to_std = DMatrix::<f64>::zeros(total_dim, total_dim); const UNIFORM_DIM: usize = 4;
let total_dim_unif = UNIFORM_DIM * assembly_dim;
let mut unif_to_std = DMatrix::<f64>::zeros(total_dim, total_dim_unif);
for n in 0..assembly_dim { for n in 0..assembly_dim {
let block_start = element_dim * n; let block_start = (element_dim * n, UNIFORM_DIM * n);
unif_to_std unif_to_std
.view_mut((block_start, block_start), (element_dim, element_dim)) .view_mut(block_start, (element_dim, UNIFORM_DIM))
.copy_from(&local_unif_to_std(state.config.column(n))); .copy_from(&local_unif_to_std(state.config.column(n)));
} }