Compare commits
2 Commits
main
...
outline-up
| Author | SHA1 | Date | |
|---|---|---|---|
|
b479e57446 | ||
|
|
8daef5f71e |
48
README.md
48
README.md
@ -17,51 +17,3 @@ Note that currently this is just the barest beginnings of the project, more of a
|
||||
* Able to run in browser (so implemented in WASM-compatible language)
|
||||
|
||||
* Produce scalable graphics of 3D diagrams, and maybe STL files (or other fabricatable file format) as well.
|
||||
|
||||
## Prototype
|
||||
|
||||
The latest prototype is in the folder `app-proto`. It includes both a user interface and a numerical constraint-solving engine.
|
||||
|
||||
### Install the prerequisites
|
||||
|
||||
1. Install [`rustup`](https://rust-lang.github.io/rustup/): the officially recommended Rust toolchain manager
|
||||
* It's available on Ubuntu as a [Snap](https://snapcraft.io/rustup)
|
||||
2. Call `rustup default stable` to "download the latest stable release of Rust and set it as your default toolchain"
|
||||
* If you forget, the `rustup` [help system](https://github.com/rust-lang/rustup/blob/d9b3601c3feb2e88cf3f8ca4f7ab4fdad71441fd/src/errors.rs#L109-L112) will remind you
|
||||
3. Call `rustup target add wasm32-unknown-unknown` to add the [most generic 32-bit WebAssembly target](https://doc.rust-lang.org/nightly/rustc/platform-support/wasm32-unknown-unknown.html)
|
||||
4. Call `cargo install wasm-pack` to install the [WebAssembly toolchain](https://rustwasm.github.io/docs/wasm-pack/)
|
||||
5. Call `cargo install trunk` to install the [Trunk](https://trunkrs.dev/) web-build tool
|
||||
6. Add the `.cargo/bin` folder in your home directory to your executable search path
|
||||
* This lets you call Trunk, and other tools installed by Cargo, without specifying their paths
|
||||
* On POSIX systems, the search path is stored in the `PATH` environment variable
|
||||
|
||||
### Play with the prototype
|
||||
|
||||
1. Go into the `app-proto` folder
|
||||
2. Call `trunk serve --release` to build and serve the prototype
|
||||
* *The crates the prototype depends on will be downloaded and served automatically*
|
||||
* *For a faster build, at the expense of a much slower prototype, you can call `trunk serve` without the `--release` flag*
|
||||
3. In a web browser, visit one of the URLs listed under the message `INFO 📡 server listening at:`
|
||||
* *Touching any file in the `app-proto` folder will make Trunk rebuild and live-reload the prototype*
|
||||
4. Press *ctrl+C* in the shell where Trunk is running to stop serving the prototype
|
||||
|
||||
### Run the engine on some example problems
|
||||
|
||||
1. Go into the `app-proto` folder
|
||||
2. Call `./run-examples`
|
||||
* *For each example problem, the engine will print the value of the loss function at each optimization step*
|
||||
* *The first example that prints is the same as the Irisawa hexlet example from the Julia version of the engine prototype. If you go into `engine-proto/gram-test`, launch Julia, and then*
|
||||
|
||||
```julia
|
||||
include("irisawa-hexlet.jl")
|
||||
for (step, scaled_loss) in enumerate(history_alt.scaled_loss)
|
||||
println(rpad(step-1, 4), " | ", scaled_loss)
|
||||
end
|
||||
```
|
||||
|
||||
*you should see that it prints basically the same loss history until the last few steps, when the lower default precision of the Rust engine really starts to show*
|
||||
|
||||
### Run the automated tests
|
||||
|
||||
1. Go into the `app-proto` folder
|
||||
2. Call `cargo test`
|
||||
|
||||
@ -6,7 +6,6 @@ edition = "2021"
|
||||
|
||||
[features]
|
||||
default = ["console_error_panic_hook"]
|
||||
dev = []
|
||||
|
||||
[dependencies]
|
||||
itertools = "0.13.0"
|
||||
@ -26,7 +25,6 @@ console_error_panic_hook = { version = "0.1.7", optional = true }
|
||||
[dependencies.web-sys]
|
||||
version = "0.3.69"
|
||||
features = [
|
||||
'DomRect',
|
||||
'HtmlCanvasElement',
|
||||
'HtmlInputElement',
|
||||
'Performance',
|
||||
@ -38,12 +36,7 @@ features = [
|
||||
'WebGlVertexArrayObject'
|
||||
]
|
||||
|
||||
# the self-dependency specifies features to use for tests and examples
|
||||
#
|
||||
# https://github.com/rust-lang/cargo/issues/2911#issuecomment-1483256987
|
||||
#
|
||||
[dev-dependencies]
|
||||
dyna3 = { path = ".", default-features = false, features = ["dev"] }
|
||||
wasm-bindgen-test = "0.3.34"
|
||||
|
||||
[profile.release]
|
||||
|
||||
@ -1,25 +0,0 @@
|
||||
use dyna3::engine::{Q, irisawa::realize_irisawa_hexlet};
|
||||
|
||||
fn main() {
|
||||
const SCALED_TOL: f64 = 1.0e-12;
|
||||
let (config, _, success, history) = realize_irisawa_hexlet(SCALED_TOL);
|
||||
print!("\nCompleted Gram matrix:{}", config.tr_mul(&*Q) * &config);
|
||||
if success {
|
||||
println!("Target accuracy achieved!");
|
||||
} else {
|
||||
println!("Failed to reach target accuracy");
|
||||
}
|
||||
println!("Steps: {}", history.scaled_loss.len() - 1);
|
||||
println!("Loss: {}", history.scaled_loss.last().unwrap());
|
||||
if success {
|
||||
println!("\nChain diameters:");
|
||||
println!(" {} sun (given)", 1.0 / config[(3, 3)]);
|
||||
for k in 4..9 {
|
||||
println!(" {} sun", 1.0 / config[(3, k)]);
|
||||
}
|
||||
}
|
||||
println!("\nStep │ Loss\n─────┼────────────────────────────────");
|
||||
for (step, scaled_loss) in history.scaled_loss.into_iter().enumerate() {
|
||||
println!("{:<4} │ {}", step, scaled_loss);
|
||||
}
|
||||
}
|
||||
@ -1,38 +0,0 @@
|
||||
use nalgebra::DMatrix;
|
||||
|
||||
use dyna3::engine::{Q, point, realize_gram, sphere, PartialMatrix};
|
||||
|
||||
fn main() {
|
||||
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)];
|
||||
println!();
|
||||
let (config, _, success, history) = realize_gram(
|
||||
&gram, guess, &frozen,
|
||||
1.0e-12, 0.5, 0.9, 1.1, 200, 110
|
||||
);
|
||||
print!("\nCompleted Gram matrix:{}", config.tr_mul(&*Q) * &config);
|
||||
print!("Configuration:{}", config);
|
||||
if success {
|
||||
println!("Target accuracy achieved!");
|
||||
} else {
|
||||
println!("Failed to reach target accuracy");
|
||||
}
|
||||
println!("Steps: {}", history.scaled_loss.len() - 1);
|
||||
println!("Loss: {}", history.scaled_loss.last().unwrap());
|
||||
println!("\nStep │ Loss\n─────┼────────────────────────────────");
|
||||
for (step, scaled_loss) in history.scaled_loss.into_iter().enumerate() {
|
||||
println!("{:<4} │ {}", step, scaled_loss);
|
||||
}
|
||||
}
|
||||
@ -1,40 +0,0 @@
|
||||
use nalgebra::DMatrix;
|
||||
|
||||
use dyna3::engine::{Q, realize_gram, sphere, PartialMatrix};
|
||||
|
||||
fn main() {
|
||||
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 = {
|
||||
let a: f64 = 0.75_f64.sqrt();
|
||||
DMatrix::from_columns(&[
|
||||
sphere(1.0, 0.0, 0.0, 1.0),
|
||||
sphere(-0.5, a, 0.0, 1.0),
|
||||
sphere(-0.5, -a, 0.0, 1.0)
|
||||
])
|
||||
};
|
||||
println!();
|
||||
let (config, _, success, history) = realize_gram(
|
||||
&gram, guess, &[],
|
||||
1.0e-12, 0.5, 0.9, 1.1, 200, 110
|
||||
);
|
||||
print!("\nCompleted Gram matrix:{}", config.tr_mul(&*Q) * &config);
|
||||
if success {
|
||||
println!("Target accuracy achieved!");
|
||||
} else {
|
||||
println!("Failed to reach target accuracy");
|
||||
}
|
||||
println!("Steps: {}", history.scaled_loss.len() - 1);
|
||||
println!("Loss: {}", history.scaled_loss.last().unwrap());
|
||||
println!("\nStep │ Loss\n─────┼────────────────────────────────");
|
||||
for (step, scaled_loss) in history.scaled_loss.into_iter().enumerate() {
|
||||
println!("{:<4} │ {}", step, scaled_loss);
|
||||
}
|
||||
}
|
||||
@ -1,11 +1,8 @@
|
||||
#!/bin/sh
|
||||
|
||||
# run all Cargo examples, as described here:
|
||||
# based on "Enabling print statements in Cargo tests", by Jon Almeida
|
||||
#
|
||||
# Karol Kuczmarski. "Add examples to your Rust libraries"
|
||||
# http://xion.io/post/code/rust-examples.html
|
||||
# https://jonalmeida.com/posts/2015/01/23/print-cargo/
|
||||
#
|
||||
|
||||
cargo run --example irisawa-hexlet
|
||||
cargo run --example three-spheres
|
||||
cargo run --example point-on-sphere
|
||||
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
|
||||
|
||||
@ -148,6 +148,9 @@ pub fn AddRemove() -> View {
|
||||
assembly.elements_by_id.update(|elts_by_id| elts_by_id.clear());
|
||||
state.selection.update(|sel| sel.clear());
|
||||
|
||||
// increment assembly serial number
|
||||
state.assembly_serial.set_fn_silent(|serial| serial.wrapping_add(1));
|
||||
|
||||
// load assembly
|
||||
match name.as_str() {
|
||||
"general" => load_gen_assemb(assembly),
|
||||
|
||||
@ -1,11 +1,11 @@
|
||||
use nalgebra::{DMatrix, DVector, DVectorView, Vector3};
|
||||
use nalgebra::{DMatrix, DVector};
|
||||
use rustc_hash::FxHashMap;
|
||||
use slab::Slab;
|
||||
use std::{collections::BTreeSet, sync::atomic::{AtomicU64, Ordering}};
|
||||
use std::collections::BTreeSet;
|
||||
use sycamore::prelude::*;
|
||||
use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */
|
||||
|
||||
use crate::engine::{realize_gram, ConfigSubspace, PartialMatrix, Q};
|
||||
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;
|
||||
@ -13,13 +13,6 @@ pub type ConstraintKey = usize;
|
||||
|
||||
pub type ElementColor = [f32; 3];
|
||||
|
||||
/* KLUDGE */
|
||||
// we should reconsider this design when we build a system for switching between
|
||||
// assemblies. at that point, we might want to switch to hierarchical keys,
|
||||
// where each each element has a key that identifies it within its assembly and
|
||||
// each assembly has a key that identifies it within the sesssion
|
||||
static NEXT_ELEMENT_SERIAL: AtomicU64 = AtomicU64::new(0);
|
||||
|
||||
#[derive(Clone, PartialEq)]
|
||||
pub struct Element {
|
||||
pub id: String,
|
||||
@ -27,15 +20,10 @@ pub struct Element {
|
||||
pub color: ElementColor,
|
||||
pub representation: Signal<DVector<f64>>,
|
||||
pub constraints: Signal<BTreeSet<ConstraintKey>>,
|
||||
|
||||
// a serial number, assigned by `Element::new`, that uniquely identifies
|
||||
// each element
|
||||
pub serial: u64,
|
||||
|
||||
// the configuration matrix column index that was assigned to this element
|
||||
// last time the assembly was realized, or `None` if the element has never
|
||||
// been through a realization
|
||||
column_index: Option<usize>
|
||||
// last time the assembly was realized
|
||||
column_index: usize
|
||||
}
|
||||
|
||||
impl Element {
|
||||
@ -45,71 +33,17 @@ impl Element {
|
||||
color: ElementColor,
|
||||
representation: DVector<f64>
|
||||
) -> Element {
|
||||
// take the next serial number, panicking if that was the last number we
|
||||
// had left. the technique we use to panic on overflow is taken from
|
||||
// _Rust Atomics and Locks_, by Mara Bos
|
||||
//
|
||||
// https://marabos.nl/atomics/atomics.html#example-handle-overflow
|
||||
//
|
||||
let serial = NEXT_ELEMENT_SERIAL.fetch_update(
|
||||
Ordering::SeqCst, Ordering::SeqCst,
|
||||
|serial| serial.checked_add(1)
|
||||
).expect("Out of serial numbers for elements");
|
||||
|
||||
Element {
|
||||
id: id,
|
||||
label: label,
|
||||
color: color,
|
||||
representation: create_signal(representation),
|
||||
constraints: create_signal(BTreeSet::default()),
|
||||
serial: serial,
|
||||
column_index: None
|
||||
}
|
||||
}
|
||||
|
||||
// the smallest positive depth, represented as a multiple of `dir`, where
|
||||
// the line generated by `dir` hits the element (which is assumed to be a
|
||||
// sphere). returns `None` if the line misses the sphere. this function
|
||||
// should be kept synchronized with `sphere_cast` in `inversive.frag`, which
|
||||
// does essentially the same thing on the GPU side
|
||||
pub fn cast(&self, dir: Vector3<f64>, assembly_to_world: &DMatrix<f64>) -> Option<f64> {
|
||||
// if `a/b` is less than this threshold, we approximate
|
||||
// `a*u^2 + b*u + c` by the linear function `b*u + c`
|
||||
const DEG_THRESHOLD: f64 = 1e-9;
|
||||
|
||||
let rep = self.representation.with_untracked(|rep| assembly_to_world * rep);
|
||||
let a = -rep[3] * dir.norm_squared();
|
||||
let b = rep.rows_range(..3).dot(&dir);
|
||||
let c = -rep[4];
|
||||
|
||||
let adjust = 4.0*a*c/(b*b);
|
||||
if adjust < 1.0 {
|
||||
// as long as `b` is non-zero, the linear approximation of
|
||||
//
|
||||
// a*u^2 + b*u + c
|
||||
//
|
||||
// at `u = 0` will reach zero at a finite depth `u_lin`. the root of
|
||||
// the quadratic adjacent to `u_lin` is stored in `lin_root`. if
|
||||
// both roots have the same sign, `lin_root` will be the one closer
|
||||
// to `u = 0`
|
||||
let square_rect_ratio = 1.0 + (1.0 - adjust).sqrt();
|
||||
let lin_root = -(2.0*c)/b / square_rect_ratio;
|
||||
if a.abs() > DEG_THRESHOLD * b.abs() {
|
||||
if lin_root > 0.0 {
|
||||
Some(lin_root)
|
||||
} else {
|
||||
let other_root = -b/(2.*a) * square_rect_ratio;
|
||||
(other_root > 0.0).then_some(other_root)
|
||||
}
|
||||
} else {
|
||||
(lin_root > 0.0).then_some(lin_root)
|
||||
}
|
||||
} else {
|
||||
// the line through `dir` misses the sphere completely
|
||||
None
|
||||
column_index: 0
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#[derive(Clone)]
|
||||
pub struct Constraint {
|
||||
@ -120,13 +54,6 @@ pub struct Constraint {
|
||||
pub active: Signal<bool>
|
||||
}
|
||||
|
||||
pub struct ElementMotion<'a> {
|
||||
pub key: ElementKey,
|
||||
pub velocity: DVectorView<'a, f64>
|
||||
}
|
||||
|
||||
type AssemblyMotion<'a> = Vec<ElementMotion<'a>>;
|
||||
|
||||
// a complete, view-independent description of an assembly
|
||||
#[derive(Clone)]
|
||||
pub struct Assembly {
|
||||
@ -134,18 +61,6 @@ pub struct Assembly {
|
||||
pub elements: Signal<Slab<Element>>,
|
||||
pub constraints: Signal<Slab<Constraint>>,
|
||||
|
||||
// solution variety tangent space. the basis vectors are stored in
|
||||
// configuration matrix format, ordered according to the elements' column
|
||||
// indices. when you realize the assembly, every element that's present
|
||||
// during realization gets a column index and is reflected in the tangent
|
||||
// space. since the methods in this module never assign column indices
|
||||
// without later realizing the assembly, we get the following invariant:
|
||||
//
|
||||
// (1) if an element has a column index, its tangent motions can be found
|
||||
// in that column of the tangent space basis matrices
|
||||
//
|
||||
pub tangent: Signal<ConfigSubspace>,
|
||||
|
||||
// indexing
|
||||
pub elements_by_id: Signal<FxHashMap<String, ElementKey>>
|
||||
}
|
||||
@ -155,7 +70,6 @@ impl Assembly {
|
||||
Assembly {
|
||||
elements: create_signal(Slab::new()),
|
||||
constraints: create_signal(Slab::new()),
|
||||
tangent: create_signal(ConfigSubspace::zero(0)),
|
||||
elements_by_id: create_signal(FxHashMap::default())
|
||||
}
|
||||
}
|
||||
@ -219,7 +133,7 @@ impl Assembly {
|
||||
// index the elements
|
||||
self.elements.update_silent(|elts| {
|
||||
for (index, (_, elt)) in elts.into_iter().enumerate() {
|
||||
elt.column_index = Some(index);
|
||||
elt.column_index = index;
|
||||
}
|
||||
});
|
||||
|
||||
@ -231,8 +145,8 @@ impl Assembly {
|
||||
for (_, cst) in csts {
|
||||
if cst.active.get_untracked() && cst.lorentz_prod_valid.get_untracked() {
|
||||
let subjects = cst.subjects;
|
||||
let row = elts[subjects.0].column_index.unwrap();
|
||||
let col = elts[subjects.1].column_index.unwrap();
|
||||
let row = elts[subjects.0].column_index;
|
||||
let col = elts[subjects.1].column_index;
|
||||
gram_to_be.push_sym(row, col, cst.lorentz_prod.get_untracked());
|
||||
}
|
||||
}
|
||||
@ -242,7 +156,7 @@ impl Assembly {
|
||||
// Gram matrix
|
||||
let mut guess_to_be = DMatrix::<f64>::zeros(5, elts.len());
|
||||
for (_, elt) in elts {
|
||||
let index = elt.column_index.unwrap();
|
||||
let index = elt.column_index;
|
||||
gram_to_be.push_sym(index, index, 1.0);
|
||||
guess_to_be.set_column(index, &elt.representation.get_clone_untracked());
|
||||
}
|
||||
@ -267,7 +181,7 @@ impl Assembly {
|
||||
}
|
||||
|
||||
// look for a configuration with the given Gram matrix
|
||||
let (config, tangent, success, history) = realize_gram(
|
||||
let (config, success, history) = realize_gram(
|
||||
&gram, guess, &[],
|
||||
1.0e-12, 0.5, 0.9, 1.1, 200, 110
|
||||
);
|
||||
@ -283,111 +197,14 @@ impl Assembly {
|
||||
));
|
||||
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()));
|
||||
console::log_2(&JsValue::from("Tangent dimension:"), &JsValue::from(tangent.dim()));
|
||||
|
||||
if success {
|
||||
// read out the solution
|
||||
for (_, elt) in self.elements.get_clone_untracked() {
|
||||
elt.representation.update(
|
||||
|rep| rep.set_column(0, &config.column(elt.column_index.unwrap()))
|
||||
|rep| rep.set_column(0, &config.column(elt.column_index))
|
||||
);
|
||||
}
|
||||
|
||||
// save the tangent space
|
||||
self.tangent.set_silent(tangent);
|
||||
}
|
||||
}
|
||||
|
||||
// --- deformation ---
|
||||
|
||||
// project the given motion to the tangent space of the solution variety and
|
||||
// move the assembly along it. the implementation is based on invariant (1)
|
||||
// from above and the following additional invariant:
|
||||
//
|
||||
// (2) if an element is affected by a constraint, it has a column index
|
||||
//
|
||||
// we have this invariant because the assembly gets realized each time you
|
||||
// add a constraint
|
||||
pub fn deform(&self, motion: AssemblyMotion) {
|
||||
/* KLUDGE */
|
||||
// when the tangent space is zero, deformation won't do anything, but
|
||||
// the attempt to deform should be registered in the UI. this console
|
||||
// message will do for now
|
||||
if self.tangent.with(|tan| tan.dim() <= 0 && tan.assembly_dim() > 0) {
|
||||
console::log_1(&JsValue::from("The assembly is rigid"));
|
||||
}
|
||||
|
||||
// give a column index to each moving element that doesn't have one yet.
|
||||
// this temporarily breaks invariant (1), but the invariant will be
|
||||
// restored when we realize the assembly at the end of the deformation.
|
||||
// in the process, we find out how many matrix columns we'll need to
|
||||
// hold the deformation
|
||||
let realized_dim = self.tangent.with(|tan| tan.assembly_dim());
|
||||
let motion_dim = self.elements.update_silent(|elts| {
|
||||
let mut next_column_index = realized_dim;
|
||||
for elt_motion in motion.iter() {
|
||||
let moving_elt = &mut elts[elt_motion.key];
|
||||
if moving_elt.column_index.is_none() {
|
||||
moving_elt.column_index = Some(next_column_index);
|
||||
next_column_index += 1;
|
||||
}
|
||||
}
|
||||
next_column_index
|
||||
});
|
||||
|
||||
// project the element motions onto the tangent space of the solution
|
||||
// variety and sum them to get a deformation of the whole assembly. the
|
||||
// matrix `motion_proj` that holds the deformation has extra columns for
|
||||
// any moving elements that aren't reflected in the saved tangent space
|
||||
const ELEMENT_DIM: usize = 5;
|
||||
let mut motion_proj = DMatrix::zeros(ELEMENT_DIM, motion_dim);
|
||||
for elt_motion in motion {
|
||||
// we can unwrap the column index because we know that every moving
|
||||
// element has one at this point
|
||||
let column_index = self.elements.with_untracked(
|
||||
|elts| elts[elt_motion.key].column_index.unwrap()
|
||||
);
|
||||
|
||||
if column_index < realized_dim {
|
||||
// this element had a column index when we started, so by
|
||||
// invariant (1), it's reflected in the tangent space
|
||||
let mut target_columns = motion_proj.columns_mut(0, realized_dim);
|
||||
target_columns += self.tangent.with(
|
||||
|tan| tan.proj(&elt_motion.velocity, column_index)
|
||||
);
|
||||
} else {
|
||||
// this element didn't have a column index when we started, so
|
||||
// by invariant (2), it's unconstrained
|
||||
let mut target_column = motion_proj.column_mut(column_index);
|
||||
target_column += elt_motion.velocity;
|
||||
}
|
||||
}
|
||||
|
||||
// step each element along the mass shell geodesic that matches its
|
||||
// velocity in the deformation found above
|
||||
/* KLUDGE */
|
||||
// since our test assemblies only include spheres, we assume that every
|
||||
// element is on the 1 mass shell
|
||||
for (_, elt) in self.elements.get_clone_untracked() {
|
||||
elt.representation.update_silent(|rep| {
|
||||
match elt.column_index {
|
||||
Some(column_index) => {
|
||||
let rep_next = &*rep + motion_proj.column(column_index);
|
||||
let normalizer = rep_next.dot(&(&*Q * &rep_next));
|
||||
rep.set_column(0, &(rep_next / normalizer));
|
||||
},
|
||||
None => {
|
||||
console::log_1(&JsValue::from(
|
||||
format!("No velocity to unpack for fresh element \"{}\"", elt.id)
|
||||
))
|
||||
}
|
||||
};
|
||||
});
|
||||
}
|
||||
|
||||
// bring the configuration back onto the solution variety. this also
|
||||
// gets the elements' column indices and the saved tangent space back in
|
||||
// sync
|
||||
self.realize();
|
||||
}
|
||||
}
|
||||
@ -1,12 +1,10 @@
|
||||
use core::array;
|
||||
use nalgebra::{DMatrix, DVector, Rotation3, Vector3};
|
||||
use nalgebra::{DMatrix, Rotation3, Vector3};
|
||||
use sycamore::{prelude::*, motion::create_raf};
|
||||
use web_sys::{
|
||||
console,
|
||||
window,
|
||||
Element,
|
||||
KeyboardEvent,
|
||||
MouseEvent,
|
||||
WebGl2RenderingContext,
|
||||
WebGlProgram,
|
||||
WebGlShader,
|
||||
@ -14,7 +12,7 @@ use web_sys::{
|
||||
wasm_bindgen::{JsCast, JsValue}
|
||||
};
|
||||
|
||||
use crate::{AppState, assembly::{ElementKey, ElementMotion}};
|
||||
use crate::AppState;
|
||||
|
||||
fn compile_shader(
|
||||
context: &WebGl2RenderingContext,
|
||||
@ -84,24 +82,6 @@ fn bind_vertex_attrib(
|
||||
);
|
||||
}
|
||||
|
||||
// the direction in camera space that a mouse event is pointing along
|
||||
fn event_dir(event: &MouseEvent) -> Vector3<f64> {
|
||||
let target: Element = event.target().unwrap().unchecked_into();
|
||||
let rect = target.get_bounding_client_rect();
|
||||
let width = rect.width();
|
||||
let height = rect.height();
|
||||
let shortdim = width.min(height);
|
||||
|
||||
// this constant should be kept synchronized with `inversive.frag`
|
||||
const FOCAL_SLOPE: f64 = 0.3;
|
||||
|
||||
Vector3::new(
|
||||
FOCAL_SLOPE * (2.0*(f64::from(event.client_x()) - rect.left()) - width) / shortdim,
|
||||
FOCAL_SLOPE * (2.0*(rect.bottom() - f64::from(event.client_y())) - height) / shortdim,
|
||||
-1.0
|
||||
)
|
||||
}
|
||||
|
||||
#[component]
|
||||
pub fn Display() -> View {
|
||||
let state = use_context::<AppState>();
|
||||
@ -109,9 +89,6 @@ pub fn Display() -> View {
|
||||
// canvas
|
||||
let display = create_node_ref();
|
||||
|
||||
// viewpoint
|
||||
let assembly_to_world = create_signal(DMatrix::<f64>::identity(5, 5));
|
||||
|
||||
// navigation
|
||||
let pitch_up = create_signal(0.0);
|
||||
let pitch_down = create_signal(0.0);
|
||||
@ -123,14 +100,6 @@ pub fn Display() -> View {
|
||||
let zoom_out = create_signal(0.0);
|
||||
let turntable = create_signal(false); /* BENCHMARKING */
|
||||
|
||||
// manipulation
|
||||
let translate_neg_x = create_signal(0.0);
|
||||
let translate_pos_x = create_signal(0.0);
|
||||
let translate_neg_y = create_signal(0.0);
|
||||
let translate_pos_y = create_signal(0.0);
|
||||
let translate_neg_z = create_signal(0.0);
|
||||
let translate_pos_z = create_signal(0.0);
|
||||
|
||||
// change listener
|
||||
let scene_changed = create_signal(true);
|
||||
create_effect(move || {
|
||||
@ -149,7 +118,6 @@ pub fn Display() -> View {
|
||||
let mut frames_since_last_sample = 0;
|
||||
let mean_frame_interval = create_signal(0.0);
|
||||
|
||||
let assembly_for_raf = state.assembly.clone();
|
||||
on_mount(move || {
|
||||
// timing
|
||||
let mut last_time = 0.0;
|
||||
@ -162,9 +130,6 @@ pub fn Display() -> View {
|
||||
let mut rotation = DMatrix::<f64>::identity(5, 5);
|
||||
let mut location_z: f64 = 5.0;
|
||||
|
||||
// manipulation
|
||||
const TRANSLATION_SPEED: f64 = 0.15; // in length units per second
|
||||
|
||||
// display parameters
|
||||
const OPACITY: f32 = 0.5; /* SCAFFOLDING */
|
||||
const HIGHLIGHT: f32 = 0.2; /* SCAFFOLDING */
|
||||
@ -285,14 +250,6 @@ pub fn Display() -> View {
|
||||
let zoom_out_val = zoom_out.get();
|
||||
let turntable_val = turntable.get(); /* BENCHMARKING */
|
||||
|
||||
// get the manipulation state
|
||||
let translate_neg_x_val = translate_neg_x.get();
|
||||
let translate_pos_x_val = translate_pos_x.get();
|
||||
let translate_neg_y_val = translate_neg_y.get();
|
||||
let translate_pos_y_val = translate_pos_y.get();
|
||||
let translate_neg_z_val = translate_neg_z.get();
|
||||
let translate_pos_z_val = translate_pos_z.get();
|
||||
|
||||
// update the assembly's orientation
|
||||
let ang_vel = {
|
||||
let pitch = pitch_up_val - pitch_down_val;
|
||||
@ -318,41 +275,6 @@ pub fn Display() -> View {
|
||||
let zoom = zoom_out_val - zoom_in_val;
|
||||
location_z *= (time_step * ZOOM_SPEED * zoom).exp();
|
||||
|
||||
// manipulate the assembly
|
||||
if state.selection.with(|sel| sel.len() == 1) {
|
||||
let sel = state.selection.with(
|
||||
|sel| *sel.into_iter().next().unwrap()
|
||||
);
|
||||
let rep = state.assembly.elements.with_untracked(
|
||||
|elts| elts[sel].representation.get_clone_untracked()
|
||||
);
|
||||
let translate_x = translate_pos_x_val - translate_neg_x_val;
|
||||
let translate_y = translate_pos_y_val - translate_neg_y_val;
|
||||
let translate_z = translate_pos_z_val - translate_neg_z_val;
|
||||
if translate_x != 0.0 || translate_y != 0.0 || translate_z != 0.0 {
|
||||
let vel_field = {
|
||||
let u = Vector3::new(translate_x, translate_y, translate_z).normalize();
|
||||
DMatrix::from_column_slice(5, 5, &[
|
||||
0.0, 0.0, 0.0, 0.0, u[0],
|
||||
0.0, 0.0, 0.0, 0.0, u[1],
|
||||
0.0, 0.0, 0.0, 0.0, u[2],
|
||||
2.0*u[0], 2.0*u[1], 2.0*u[2], 0.0, 0.0,
|
||||
0.0, 0.0, 0.0, 0.0, 0.0
|
||||
])
|
||||
};
|
||||
let elt_motion: DVector<f64> = time_step * TRANSLATION_SPEED * vel_field * rep;
|
||||
assembly_for_raf.deform(
|
||||
vec![
|
||||
ElementMotion {
|
||||
key: sel,
|
||||
velocity: elt_motion.as_view()
|
||||
}
|
||||
]
|
||||
);
|
||||
scene_changed.set(true);
|
||||
}
|
||||
}
|
||||
|
||||
if scene_changed.get() {
|
||||
/* INSTRUMENTS */
|
||||
// measure mean frame interval
|
||||
@ -374,7 +296,7 @@ pub fn Display() -> View {
|
||||
0.0, 0.0, 0.0, 0.0, 1.0
|
||||
])
|
||||
};
|
||||
let asm_to_world = &location * &orientation;
|
||||
let assembly_to_world = &location * &orientation;
|
||||
|
||||
// get the assembly
|
||||
let (
|
||||
@ -389,7 +311,7 @@ pub fn Display() -> View {
|
||||
|
||||
// representation vectors in world coordinates
|
||||
elts.iter().map(
|
||||
|(_, elt)| elt.representation.with(|rep| &asm_to_world * rep)
|
||||
|(_, elt)| elt.representation.with(|rep| &assembly_to_world * rep)
|
||||
).collect::<Vec<_>>(),
|
||||
|
||||
// colors
|
||||
@ -448,9 +370,6 @@ pub fn Display() -> View {
|
||||
// draw the scene
|
||||
ctx.draw_arrays(WebGl2RenderingContext::TRIANGLES, 0, VERTEX_CNT as i32);
|
||||
|
||||
// update the viewpoint
|
||||
assembly_to_world.set(asm_to_world);
|
||||
|
||||
// clear the scene change flag
|
||||
scene_changed.set(
|
||||
pitch_up_val != 0.0
|
||||
@ -471,7 +390,7 @@ pub fn Display() -> View {
|
||||
start_animation_loop();
|
||||
});
|
||||
|
||||
let set_nav_signal = move |event: &KeyboardEvent, value: f64| {
|
||||
let set_nav_signal = move |event: KeyboardEvent, value: f64| {
|
||||
let mut navigating = true;
|
||||
let shift = event.shift_key();
|
||||
match event.key().as_str() {
|
||||
@ -491,23 +410,6 @@ pub fn Display() -> View {
|
||||
}
|
||||
};
|
||||
|
||||
let set_manip_signal = move |event: &KeyboardEvent, value: f64| {
|
||||
let mut manipulating = true;
|
||||
let shift = event.shift_key();
|
||||
match event.key().as_str() {
|
||||
"d" | "D" => translate_pos_x.set(value),
|
||||
"a" | "A" => translate_neg_x.set(value),
|
||||
"w" | "W" if shift => translate_neg_z.set(value),
|
||||
"s" | "S" if shift => translate_pos_z.set(value),
|
||||
"w" | "W" => translate_pos_y.set(value),
|
||||
"s" | "S" => translate_neg_y.set(value),
|
||||
_ => manipulating = false
|
||||
};
|
||||
if manipulating {
|
||||
event.prevent_default();
|
||||
}
|
||||
};
|
||||
|
||||
view! {
|
||||
/* TO DO */
|
||||
// switch back to integer-valued parameters when that becomes possible
|
||||
@ -519,7 +421,6 @@ pub fn Display() -> View {
|
||||
tabindex="0",
|
||||
on:keydown=move |event: KeyboardEvent| {
|
||||
if event.key() == "Shift" {
|
||||
// swap navigation inputs
|
||||
roll_cw.set(yaw_right.get());
|
||||
roll_ccw.set(yaw_left.get());
|
||||
zoom_in.set(pitch_up.get());
|
||||
@ -528,24 +429,16 @@ pub fn Display() -> View {
|
||||
yaw_left.set(0.0);
|
||||
pitch_up.set(0.0);
|
||||
pitch_down.set(0.0);
|
||||
|
||||
// swap manipulation inputs
|
||||
translate_pos_z.set(translate_neg_y.get());
|
||||
translate_neg_z.set(translate_pos_y.get());
|
||||
translate_pos_y.set(0.0);
|
||||
translate_neg_y.set(0.0);
|
||||
} else {
|
||||
if event.key() == "Enter" { /* BENCHMARKING */
|
||||
turntable.set_fn(|turn| !turn);
|
||||
scene_changed.set(true);
|
||||
}
|
||||
set_nav_signal(&event, 1.0);
|
||||
set_manip_signal(&event, 1.0);
|
||||
set_nav_signal(event, 1.0);
|
||||
}
|
||||
},
|
||||
on:keyup=move |event: KeyboardEvent| {
|
||||
if event.key() == "Shift" {
|
||||
// swap navigation inputs
|
||||
yaw_right.set(roll_cw.get());
|
||||
yaw_left.set(roll_ccw.get());
|
||||
pitch_up.set(zoom_in.get());
|
||||
@ -554,15 +447,8 @@ pub fn Display() -> View {
|
||||
roll_ccw.set(0.0);
|
||||
zoom_in.set(0.0);
|
||||
zoom_out.set(0.0);
|
||||
|
||||
// swap manipulation inputs
|
||||
translate_pos_y.set(translate_neg_z.get());
|
||||
translate_neg_y.set(translate_pos_z.get());
|
||||
translate_pos_z.set(0.0);
|
||||
translate_neg_z.set(0.0);
|
||||
} else {
|
||||
set_nav_signal(&event, 0.0);
|
||||
set_manip_signal(&event, 0.0);
|
||||
set_nav_signal(event, 0.0);
|
||||
}
|
||||
},
|
||||
on:blur=move |_| {
|
||||
@ -572,31 +458,6 @@ pub fn Display() -> View {
|
||||
yaw_left.set(0.0);
|
||||
roll_ccw.set(0.0);
|
||||
roll_cw.set(0.0);
|
||||
},
|
||||
on:click=move |event: MouseEvent| {
|
||||
// find the nearest element along the pointer direction
|
||||
let dir = event_dir(&event);
|
||||
console::log_1(&JsValue::from(dir.to_string()));
|
||||
let mut clicked: Option<(ElementKey, f64)> = None;
|
||||
for (key, elt) in state.assembly.elements.get_clone_untracked() {
|
||||
match assembly_to_world.with(|asm_to_world| elt.cast(dir, asm_to_world)) {
|
||||
Some(depth) => match clicked {
|
||||
Some((_, best_depth)) => {
|
||||
if depth < best_depth {
|
||||
clicked = Some((key, depth))
|
||||
}
|
||||
},
|
||||
None => clicked = Some((key, depth))
|
||||
}
|
||||
None => ()
|
||||
};
|
||||
}
|
||||
|
||||
// if we clicked something, select it
|
||||
match clicked {
|
||||
Some((key, _)) => state.select(key, event.shift_key()),
|
||||
None => state.selection.update(|sel| sel.clear())
|
||||
};
|
||||
}
|
||||
)
|
||||
}
|
||||
|
||||
@ -1,10 +1,10 @@
|
||||
use lazy_static::lazy_static;
|
||||
use nalgebra::{Const, DMatrix, DVector, DVectorView, Dyn, SymmetricEigen};
|
||||
use nalgebra::{Const, DMatrix, DVector, Dyn};
|
||||
use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */
|
||||
|
||||
// --- elements ---
|
||||
|
||||
#[cfg(feature = "dev")]
|
||||
#[cfg(test)]
|
||||
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)])
|
||||
}
|
||||
@ -85,75 +85,6 @@ impl PartialMatrix {
|
||||
}
|
||||
}
|
||||
|
||||
// --- configuration subspaces ---
|
||||
|
||||
#[derive(Clone)]
|
||||
pub struct ConfigSubspace {
|
||||
assembly_dim: usize,
|
||||
basis: Vec<DMatrix<f64>>
|
||||
}
|
||||
|
||||
impl ConfigSubspace {
|
||||
pub fn zero(assembly_dim: usize) -> ConfigSubspace {
|
||||
ConfigSubspace {
|
||||
assembly_dim: assembly_dim,
|
||||
basis: 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>, assembly_dim: usize) -> ConfigSubspace {
|
||||
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 = eig_pairs.filter_map(
|
||||
|(λ, v)| (λ.abs() < THRESHOLD).then_some(
|
||||
Into::<DMatrix<f64>>::into(
|
||||
v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim))
|
||||
)
|
||||
)
|
||||
);
|
||||
|
||||
/* 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)
|
||||
));
|
||||
|
||||
ConfigSubspace {
|
||||
assembly_dim: assembly_dim,
|
||||
basis: basis.collect()
|
||||
}
|
||||
}
|
||||
|
||||
pub fn dim(&self) -> usize {
|
||||
self.basis.len()
|
||||
}
|
||||
|
||||
pub fn assembly_dim(&self) -> usize {
|
||||
self.assembly_dim
|
||||
}
|
||||
|
||||
// find the projection onto this subspace, with respect to the Euclidean
|
||||
// inner product, of the motion where the element with the given column
|
||||
// index has velocity `v`
|
||||
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.iter().map(
|
||||
|b| b.column(column_index).dot(&v) * b
|
||||
).sum()
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// --- descent history ---
|
||||
|
||||
pub struct DescentHistory {
|
||||
@ -182,7 +113,7 @@ impl DescentHistory {
|
||||
|
||||
// the Lorentz form
|
||||
lazy_static! {
|
||||
pub static ref Q: DMatrix<f64> = DMatrix::from_row_slice(5, 5, &[
|
||||
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,
|
||||
@ -250,7 +181,7 @@ pub fn realize_gram(
|
||||
reg_scale: f64,
|
||||
max_descent_steps: i32,
|
||||
max_backoff_steps: i32
|
||||
) -> (DMatrix<f64>, ConfigSubspace, bool, DescentHistory) {
|
||||
) -> (DMatrix<f64>, bool, DescentHistory) {
|
||||
// start the descent history
|
||||
let mut history = DescentHistory::new();
|
||||
|
||||
@ -270,8 +201,12 @@ pub fn realize_gram(
|
||||
|
||||
// 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 {
|
||||
// stop if the loss is tolerably low
|
||||
history.config.push(state.config.clone());
|
||||
history.scaled_loss.push(state.loss / scale_adjustment);
|
||||
if state.loss < tol { break; }
|
||||
|
||||
// find the negative gradient of the loss function
|
||||
let neg_grad = 4.0 * &*Q * &state.config * &state.err_proj;
|
||||
let mut neg_grad_stacked = neg_grad.clone().reshape_generic(Dyn(total_dim), Const::<1>);
|
||||
@ -294,7 +229,7 @@ pub fn realize_gram(
|
||||
hess_cols.push(deriv_grad.reshape_generic(Dyn(total_dim), Const::<1>));
|
||||
}
|
||||
}
|
||||
hess = DMatrix::from_columns(hess_cols.as_slice());
|
||||
let mut hess = DMatrix::from_columns(hess_cols.as_slice());
|
||||
|
||||
// regularize the Hessian
|
||||
let min_eigval = hess.symmetric_eigenvalues().min();
|
||||
@ -314,11 +249,6 @@ pub fn realize_gram(
|
||||
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
|
||||
@ -326,7 +256,7 @@ pub fn realize_gram(
|
||||
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_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());
|
||||
|
||||
@ -339,93 +269,20 @@ pub fn realize_gram(
|
||||
state = better_state;
|
||||
history.backoff_steps.push(backoff_steps);
|
||||
},
|
||||
None => return (state.config, ConfigSubspace::zero(assembly_dim), false, history)
|
||||
None => return (state.config, false, history)
|
||||
};
|
||||
}
|
||||
let success = state.loss < tol;
|
||||
let tangent = if success {
|
||||
ConfigSubspace::symmetric_kernel(hess, assembly_dim)
|
||||
} else {
|
||||
ConfigSubspace::zero(assembly_dim)
|
||||
};
|
||||
(state.config, tangent, success, history)
|
||||
(state.config, state.loss < tol, 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 {
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
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![
|
||||
@ -471,75 +328,182 @@ 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() {
|
||||
// solve Irisawa's problem
|
||||
let gram = PartialMatrix({
|
||||
let mut entries = Vec::<MatrixEntry>::new();
|
||||
for s in 0..9 {
|
||||
// each sphere is represented by a spacelike vector
|
||||
entries.push(MatrixEntry { index: (s, s), value: 1.0 });
|
||||
|
||||
// the circumscribing sphere is tangent to all of the other
|
||||
// spheres, with matching orientation
|
||||
if s > 0 {
|
||||
entries.push(MatrixEntry { index: (0, s), value: 1.0 });
|
||||
entries.push(MatrixEntry { index: (s, 0), value: 1.0 });
|
||||
}
|
||||
|
||||
if s > 2 {
|
||||
// each chain sphere is tangent to the "sun" and "moon"
|
||||
// spheres, with opposing orientation
|
||||
for n in 1..3 {
|
||||
entries.push(MatrixEntry { index: (s, n), value: -1.0 });
|
||||
entries.push(MatrixEntry { index: (n, s), value: -1.0 });
|
||||
}
|
||||
|
||||
// each chain sphere is tangent to the next chain sphere,
|
||||
// with opposing orientation
|
||||
let s_next = 3 + (s-2) % 6;
|
||||
entries.push(MatrixEntry { index: (s, s_next), value: -1.0 });
|
||||
entries.push(MatrixEntry { index: (s_next, s), value: -1.0 });
|
||||
}
|
||||
}
|
||||
entries
|
||||
});
|
||||
let guess = DMatrix::from_columns(
|
||||
[
|
||||
sphere(0.0, 0.0, 0.0, 15.0),
|
||||
sphere(0.0, 0.0, -9.0, 5.0),
|
||||
sphere(0.0, 0.0, 11.0, 3.0)
|
||||
].into_iter().chain(
|
||||
(1..=6).map(
|
||||
|k| {
|
||||
let ang = (k as f64) * PI/3.0;
|
||||
sphere(9.0 * ang.cos(), 9.0 * ang.sin(), 0.0, 2.5)
|
||||
}
|
||||
)
|
||||
).collect::<Vec<_>>().as_slice()
|
||||
);
|
||||
let frozen: [(usize, usize); 4] = array::from_fn(|k| (3, k));
|
||||
const SCALED_TOL: f64 = 1.0e-12;
|
||||
let (config, _, _, _) = realize_irisawa_hexlet(SCALED_TOL);
|
||||
|
||||
// check against Irisawa's solution
|
||||
let (config, success, history) = realize_gram(
|
||||
&gram, guess, &frozen,
|
||||
SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
|
||||
);
|
||||
let entry_tol = SCALED_TOL.sqrt();
|
||||
let solution_diams = [30.0, 10.0, 6.0, 5.0, 15.0, 10.0, 3.75, 2.5, 2.0 + 8.0/11.0];
|
||||
for (k, diam) in solution_diams.into_iter().enumerate() {
|
||||
assert!((config[(3, k)] - 1.0 / diam).abs() < entry_tol);
|
||||
}
|
||||
print!("\nCompleted Gram matrix:{}", config.tr_mul(&*Q) * &config);
|
||||
if success {
|
||||
println!("Target accuracy achieved!");
|
||||
} else {
|
||||
println!("Failed to reach target accuracy");
|
||||
}
|
||||
println!("Steps: {}", history.scaled_loss.len() - 1);
|
||||
println!("Loss: {}", history.scaled_loss.last().unwrap());
|
||||
if success {
|
||||
println!("\nChain diameters:");
|
||||
println!(" {} sun (given)", 1.0 / config[(3, 3)]);
|
||||
for k in 4..9 {
|
||||
println!(" {} sun", 1.0 / config[(3, k)]);
|
||||
}
|
||||
}
|
||||
println!("\nStep │ Loss\n─────┼────────────────────────────────");
|
||||
for (step, scaled_loss) in history.scaled_loss.into_iter().enumerate() {
|
||||
println!("{:<4} │ {}", step, scaled_loss);
|
||||
}
|
||||
}
|
||||
|
||||
// --- process inspection examples ---
|
||||
|
||||
// these tests are meant for human inspection, not automated use. run them
|
||||
// one at a time in `--nocapture` mode and read through the results and
|
||||
// optimization histories that they print out. the `run-examples` script
|
||||
// will run all of them
|
||||
|
||||
#[test]
|
||||
fn three_spheres_example() {
|
||||
let gram = PartialMatrix({
|
||||
let mut entries = Vec::<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 guess = {
|
||||
let a: f64 = 0.75_f64.sqrt();
|
||||
DMatrix::from_columns(&[
|
||||
sphere(1.0, 0.0, 0.0, 1.0),
|
||||
sphere(-0.5, a, 0.0, 1.0),
|
||||
sphere(-0.5, -a, 0.0, 1.0)
|
||||
])
|
||||
};
|
||||
println!();
|
||||
let (config, success, history) = realize_gram(
|
||||
&gram, guess, &[],
|
||||
1.0e-12, 0.5, 0.9, 1.1, 200, 110
|
||||
);
|
||||
print!("\nCompleted Gram matrix:{}", config.tr_mul(&*Q) * &config);
|
||||
if success {
|
||||
println!("Target accuracy achieved!");
|
||||
} else {
|
||||
println!("Failed to reach target accuracy");
|
||||
}
|
||||
println!("Steps: {}", history.scaled_loss.len() - 1);
|
||||
println!("Loss: {}", history.scaled_loss.last().unwrap());
|
||||
println!("\nStep │ Loss\n─────┼────────────────────────────────");
|
||||
for (step, scaled_loss) in history.scaled_loss.into_iter().enumerate() {
|
||||
println!("{:<4} │ {}", step, scaled_loss);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn 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 });
|
||||
fn point_on_sphere_example() {
|
||||
let gram = PartialMatrix({
|
||||
let mut entries = Vec::<MatrixEntry>::new();
|
||||
for j in 0..2 {
|
||||
for k in 0..2 {
|
||||
entries.push(MatrixEntry {
|
||||
index: (j, k),
|
||||
value: if (j, k) == (1, 1) { 1.0 } else { 0.0 }
|
||||
});
|
||||
}
|
||||
}
|
||||
gram_to_be
|
||||
};
|
||||
entries
|
||||
});
|
||||
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)
|
||||
point(0.0, 0.0, 2.0),
|
||||
sphere(0.0, 0.0, 0.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
|
||||
let frozen = [(3, 0)];
|
||||
println!();
|
||||
let (config, success, history) = realize_gram(
|
||||
&gram, guess, &frozen,
|
||||
1.0e-12, 0.5, 0.9, 1.1, 200, 110
|
||||
);
|
||||
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);
|
||||
print!("\nCompleted Gram matrix:{}", config.tr_mul(&*Q) * &config);
|
||||
print!("Configuration:{}", config);
|
||||
if success {
|
||||
println!("Target accuracy achieved!");
|
||||
} else {
|
||||
println!("Failed to reach target accuracy");
|
||||
}
|
||||
println!("Steps: {}", history.scaled_loss.len() - 1);
|
||||
println!("Loss: {}", history.scaled_loss.last().unwrap());
|
||||
println!("\nStep │ Loss\n─────┼────────────────────────────────");
|
||||
for (step, scaled_loss) in history.scaled_loss.into_iter().enumerate() {
|
||||
println!("{:<4} │ {}", step, scaled_loss);
|
||||
}
|
||||
}
|
||||
|
||||
/* 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]
|
||||
@ -559,7 +523,7 @@ mod tests {
|
||||
]);
|
||||
let frozen = [(3, 0), (3, 1)];
|
||||
println!();
|
||||
let (config, _, success, history) = realize_gram(
|
||||
let (config, success, history) = realize_gram(
|
||||
&gram, guess.clone(), &frozen,
|
||||
1.0e-12, 0.5, 0.9, 1.1, 200, 110
|
||||
);
|
||||
|
||||
@ -1 +0,0 @@
|
||||
pub mod engine;
|
||||
@ -15,6 +15,7 @@ use outline::Outline;
|
||||
#[derive(Clone)]
|
||||
struct AppState {
|
||||
assembly: Assembly,
|
||||
assembly_serial: Signal<u32>,
|
||||
selection: Signal<FxHashSet<ElementKey>>
|
||||
}
|
||||
|
||||
@ -22,34 +23,13 @@ impl AppState {
|
||||
fn new() -> AppState {
|
||||
AppState {
|
||||
assembly: Assembly::new(),
|
||||
assembly_serial: create_signal(0),
|
||||
selection: create_signal(FxHashSet::default())
|
||||
}
|
||||
}
|
||||
|
||||
// in single-selection mode, select the element with the given key. in
|
||||
// multiple-selection mode, toggle whether the element with the given key
|
||||
// is selected
|
||||
fn select(&self, key: ElementKey, multi: bool) {
|
||||
if multi {
|
||||
self.selection.update(|sel| {
|
||||
if !sel.remove(&key) {
|
||||
sel.insert(key);
|
||||
}
|
||||
});
|
||||
} else {
|
||||
self.selection.update(|sel| {
|
||||
sel.clear();
|
||||
sel.insert(key);
|
||||
});
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
fn main() {
|
||||
// set the console error panic hook
|
||||
#[cfg(feature = "console_error_panic_hook")]
|
||||
console_error_panic_hook::set_once();
|
||||
|
||||
sycamore::render(|| {
|
||||
provide_context(AppState::new());
|
||||
|
||||
|
||||
@ -83,7 +83,18 @@ fn ElementOutlineItem(key: ElementKey, element: assembly::Element) -> View {
|
||||
move |event: KeyboardEvent| {
|
||||
match event.key().as_str() {
|
||||
"Enter" => {
|
||||
state.select(key, event.shift_key());
|
||||
if event.shift_key() {
|
||||
state.selection.update(|sel| {
|
||||
if !sel.remove(&key) {
|
||||
sel.insert(key);
|
||||
}
|
||||
});
|
||||
} else {
|
||||
state.selection.update(|sel| {
|
||||
sel.clear();
|
||||
sel.insert(key);
|
||||
});
|
||||
}
|
||||
event.prevent_default();
|
||||
},
|
||||
"ArrowRight" if constrained.get() => {
|
||||
@ -169,11 +180,15 @@ pub fn Outline() -> View {
|
||||
|
||||
// list the elements alphabetically by ID
|
||||
let element_list = state.assembly.elements.map(
|
||||
|elts| elts
|
||||
.clone()
|
||||
.into_iter()
|
||||
.sorted_by_key(|(_, elt)| elt.id.clone())
|
||||
.collect()
|
||||
move |elts| {
|
||||
let asm_serial = state.assembly_serial.get_untracked();
|
||||
elts
|
||||
.clone()
|
||||
.into_iter()
|
||||
.sorted_by_key(|(_, elt)| elt.id.clone())
|
||||
.map(|(key, elt)| (asm_serial, key, elt))
|
||||
.collect()
|
||||
}
|
||||
);
|
||||
|
||||
view! {
|
||||
@ -186,10 +201,10 @@ pub fn Outline() -> View {
|
||||
) {
|
||||
Keyed(
|
||||
list=element_list,
|
||||
view=|(key, elt)| view! {
|
||||
view=|(_, key, elt)| view! {
|
||||
ElementOutlineItem(key=key, element=elt)
|
||||
},
|
||||
key=|(_, elt)| elt.serial
|
||||
key=|(asm_serial, key, _)| (asm_serial.clone(), key.clone())
|
||||
)
|
||||
}
|
||||
}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user