Curvature regulators #80

Merged
glen merged 21 commits from Vectornaut/dyna3:curvature-regulators into main 2025-04-21 23:40:43 +00:00
6 changed files with 640 additions and 331 deletions

View file

@ -1,26 +1,19 @@
use nalgebra::DMatrix;
use dyna3::engine::{Q, point, realize_gram, sphere, PartialMatrix};
use dyna3::engine::{Q, point, realize_gram, sphere, ConstraintProblem};
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(&[
let mut problem = ConstraintProblem::from_guess(&[
point(0.0, 0.0, 2.0),
sphere(0.0, 0.0, 0.0, 1.0)
]);
let frozen = [(3, 0)];
for j in 0..2 {
for k in j..2 {
problem.gram.push_sym(j, k, if (j, k) == (1, 1) { 1.0 } else { 0.0 });
}
}
problem.frozen.push(3, 0, problem.guess[(3, 0)]);
println!();
let (config, _, success, history) = realize_gram(
&gram, guess, &frozen,
1.0e-12, 0.5, 0.9, 1.1, 200, 110
&problem, 1.0e-12, 0.5, 0.9, 1.1, 200, 110
);
print!("\nCompleted Gram matrix:{}", config.tr_mul(&*Q) * &config);
print!("Configuration:{}", config);

View file

@ -1,29 +1,22 @@
use nalgebra::DMatrix;
use dyna3::engine::{Q, realize_gram, sphere, PartialMatrix};
use dyna3::engine::{Q, realize_gram, sphere, ConstraintProblem};
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 mut problem = ConstraintProblem::from_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)
])
};
]
});
for j in 0..3 {
for k in j..3 {
problem.gram.push_sym(j, k, if j == k { 1.0 } else { -1.0 });
}
}
println!();
let (config, _, success, history) = realize_gram(
&gram, guess, &[],
1.0e-12, 0.5, 0.9, 1.1, 200, 110
&problem, 1.0e-12, 0.5, 0.9, 1.1, 200, 110
);
print!("\nCompleted Gram matrix:{}", config.tr_mul(&*Q) * &config);
if success {

View file

@ -4,14 +4,14 @@ use web_sys::{console, wasm_bindgen::JsValue};
use crate::{
engine,
AppState,
assembly::{Assembly, Element}
assembly::{Assembly, Element, InversiveDistanceRegulator}
};
/* DEBUG */
// load an example assembly for testing. this code will be removed once we've
// built a more formal test assembly system
fn load_gen_assemb(assembly: &Assembly) {
let _ = assembly.try_insert_element(
let _ = assembly.try_insert_sphere(
Element::new(
String::from("gemini_a"),
String::from("Castor"),
@ -19,7 +19,7 @@ fn load_gen_assemb(assembly: &Assembly) {
engine::sphere(0.5, 0.5, 0.0, 1.0)
)
);
let _ = assembly.try_insert_element(
let _ = assembly.try_insert_sphere(
Element::new(
String::from("gemini_b"),
String::from("Pollux"),
@ -27,7 +27,7 @@ fn load_gen_assemb(assembly: &Assembly) {
engine::sphere(-0.5, -0.5, 0.0, 1.0)
)
);
let _ = assembly.try_insert_element(
let _ = assembly.try_insert_sphere(
Element::new(
String::from("ursa_major"),
String::from("Ursa major"),
@ -35,7 +35,7 @@ fn load_gen_assemb(assembly: &Assembly) {
engine::sphere(-0.5, 0.5, 0.0, 0.75)
)
);
let _ = assembly.try_insert_element(
let _ = assembly.try_insert_sphere(
Element::new(
String::from("ursa_minor"),
String::from("Ursa minor"),
@ -43,7 +43,7 @@ fn load_gen_assemb(assembly: &Assembly) {
engine::sphere(0.5, -0.5, 0.0, 0.5)
)
);
let _ = assembly.try_insert_element(
let _ = assembly.try_insert_sphere(
Element::new(
String::from("moon_deimos"),
String::from("Deimos"),
@ -51,7 +51,7 @@ fn load_gen_assemb(assembly: &Assembly) {
engine::sphere(0.0, 0.15, 1.0, 0.25)
)
);
let _ = assembly.try_insert_element(
let _ = assembly.try_insert_sphere(
Element::new(
String::from("moon_phobos"),
String::from("Phobos"),
@ -66,7 +66,7 @@ fn load_gen_assemb(assembly: &Assembly) {
// built a more formal test assembly system
fn load_low_curv_assemb(assembly: &Assembly) {
let a = 0.75_f64.sqrt();
let _ = assembly.try_insert_element(
let _ = assembly.try_insert_sphere(
Element::new(
"central".to_string(),
"Central".to_string(),
@ -74,7 +74,7 @@ fn load_low_curv_assemb(assembly: &Assembly) {
engine::sphere(0.0, 0.0, 0.0, 1.0)
)
);
let _ = assembly.try_insert_element(
let _ = assembly.try_insert_sphere(
Element::new(
"assemb_plane".to_string(),
"Assembly plane".to_string(),
@ -82,7 +82,7 @@ fn load_low_curv_assemb(assembly: &Assembly) {
engine::sphere_with_offset(0.0, 0.0, 1.0, 0.0, 0.0)
)
);
let _ = assembly.try_insert_element(
let _ = assembly.try_insert_sphere(
Element::new(
"side1".to_string(),
"Side 1".to_string(),
@ -90,7 +90,7 @@ fn load_low_curv_assemb(assembly: &Assembly) {
engine::sphere_with_offset(1.0, 0.0, 0.0, 1.0, 0.0)
)
);
let _ = assembly.try_insert_element(
let _ = assembly.try_insert_sphere(
Element::new(
"side2".to_string(),
"Side 2".to_string(),
@ -98,7 +98,7 @@ fn load_low_curv_assemb(assembly: &Assembly) {
engine::sphere_with_offset(-0.5, a, 0.0, 1.0, 0.0)
)
);
let _ = assembly.try_insert_element(
let _ = assembly.try_insert_sphere(
Element::new(
"side3".to_string(),
"Side 3".to_string(),
@ -106,7 +106,7 @@ fn load_low_curv_assemb(assembly: &Assembly) {
engine::sphere_with_offset(-0.5, -a, 0.0, 1.0, 0.0)
)
);
let _ = assembly.try_insert_element(
let _ = assembly.try_insert_sphere(
Element::new(
"corner1".to_string(),
"Corner 1".to_string(),
@ -114,7 +114,7 @@ fn load_low_curv_assemb(assembly: &Assembly) {
engine::sphere(-4.0/3.0, 0.0, 0.0, 1.0/3.0)
)
);
let _ = assembly.try_insert_element(
let _ = assembly.try_insert_sphere(
Element::new(
"corner2".to_string(),
"Corner 2".to_string(),
@ -122,7 +122,7 @@ fn load_low_curv_assemb(assembly: &Assembly) {
engine::sphere(2.0/3.0, -4.0/3.0 * a, 0.0, 1.0/3.0)
)
);
let _ = assembly.try_insert_element(
let _ = assembly.try_insert_sphere(
Element::new(
String::from("corner3"),
String::from("Corner 3"),
@ -148,6 +148,7 @@ pub fn AddRemove() -> View {
let assembly = &state.assembly;
// clear state
assembly.regulators.update(|regs| regs.clear());
assembly.elements.update(|elts| elts.clear());
assembly.elements_by_id.update(|elts_by_id| elts_by_id.clear());
state.selection.update(|sel| sel.clear());
@ -166,18 +167,7 @@ pub fn AddRemove() -> View {
button(
on:click=|_| {
let state = use_context::<AppState>();
state.assembly.insert_new_element();
/* DEBUG */
// print updated list of elements by identifier
console::log_1(&JsValue::from("elements by identifier:"));
for (id, key) in state.assembly.elements_by_id.get_clone().iter() {
console::log_3(
&JsValue::from(" "),
&JsValue::from(id),
&JsValue::from(*key)
);
}
state.assembly.insert_new_sphere();
}
) { "+" }
button(
@ -188,13 +178,20 @@ pub fn AddRemove() -> View {
},
on:click=|_| {
let state = use_context::<AppState>();
let subjects = state.selection.with(
|sel| {
let subject_vec: Vec<_> = sel.into_iter().collect();
(subject_vec[0].clone(), subject_vec[1].clone())
}
let subjects: [_; 2] = state.selection.with(
// the button is only enabled when two elements are
// selected, so we know the cast to a two-element array
// will succeed
|sel| sel
.clone()
.into_iter()
.collect::<Vec<_>>()
.try_into()
.unwrap()
);
state.assembly.insert_regulator(
InversiveDistanceRegulator::new(subjects, &state.assembly)
);
state.assembly.insert_new_regulator(subjects);
state.selection.update(|sel| sel.clear());
}
) { "🔗" }

View file

@ -1,12 +1,21 @@
use nalgebra::{DMatrix, DVector, DVectorView, Vector3};
use rustc_hash::FxHashMap;
use slab::Slab;
use std::{collections::BTreeSet, sync::atomic::{AtomicU64, Ordering}};
use std::{collections::BTreeSet, rc::Rc, sync::atomic::{AtomicU64, Ordering}};
use sycamore::prelude::*;
use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */
use crate::{
engine::{Q, local_unif_to_std, realize_gram, ConfigSubspace, PartialMatrix},
engine::{
Q,
change_half_curvature,
local_unif_to_std,
realize_gram,
sphere,
ConfigSubspace,
ConstraintProblem
},
outline::OutlineItem,
specified::SpecifiedValue
};
@ -23,6 +32,10 @@ pub type ElementColor = [f32; 3];
// each assembly has a key that identifies it within the sesssion
static NEXT_ELEMENT_SERIAL: AtomicU64 = AtomicU64::new(0);
pub trait ProblemPoser {
fn pose(&self, problem: &mut ConstraintProblem, elts: &Slab<Element>);
}
#[derive(Clone, PartialEq)]
pub struct Element {
pub id: String,
@ -30,8 +43,8 @@ pub struct Element {
pub color: ElementColor,
pub representation: Signal<DVector<f64>>,
// All regulators with this element as a subject. The assembly owning
// this element is responsible for keeping this set up to date.
// the regulators this element is subject to. the assembly that owns the
// element is responsible for keeping this set up to date
pub regulators: Signal<BTreeSet<RegulatorKey>>,
// a serial number, assigned by `Element::new`, that uniquely identifies
@ -45,6 +58,8 @@ pub struct Element {
}
impl Element {
const CURVATURE_COMPONENT: usize = 3;
pub fn new(
id: String,
label: String,
@ -117,13 +132,148 @@ impl Element {
}
}
#[derive(Clone, Copy)]
pub struct Regulator {
pub subjects: (ElementKey, ElementKey),
impl ProblemPoser for Element {
fn pose(&self, problem: &mut ConstraintProblem, _elts: &Slab<Element>) {
let index = self.column_index.expect(
glen marked this conversation as resolved Outdated

The check here and the potential panic if it fails (and analogous code in two other ProblemPoser implementations) is evidence of misfactoring: in fact, the ProblemPosers only ever get to execute when every element has just been assigned a column index. So there is no need to have such a check.

I am open to any reasonable way to factor these checks out. One is to insist that all elements in an assembly always have a column index. As far as I can see, the only thing that None column indices are currently used for is a transient period between an Element being added to an Assembly and the next call to that Assembly being realized. That period could be eliminated by simply assigning each Element the next available column index and immediately calling realize. Or there could be a method (likely in engine?) that takes a realized assembly and a new completely unconstrained element and updates the realization to include the unconstrained element -- that should presumably be far less work than a full call to realize. As far as I can tell, the only consistency bit that needs to be kept up is the tangent space bases in the realization. But I suppose it should be possible to calculate the new tangent space basis matriceso directly in the case of a new unconstrained element added to a realized assembly.

Or we can insist that elements always have an index by assigning indices on creation, and dropping the invariant that an element with an index has to have tangent motions in the corresponding column of the tangent space basis matrices. That might not be so terrible -- if an element's column is out of range for the tangent space basis matrices, that must mean it is totally unconstrained, and so we can actually deform it freely, and so just switch to simpler deformation code for that case.

Another possibility is to initially use column index -1 for a new element, and not worry about any -1s in the problem posing because the code is structured such that they never occur when a problem is being posed; we only need to worry about the -1s when deforming, and handle the case there.

Another possibility (if I understand Rust and the Option type properly) is to simply use unwrap in these posers. It nominally could panic, but we know it won't because problems are only posed just after every element has been assigned an index. This route seems the least satisfying because we still have all of the baggage of Option even though the other possibilities make it clear that Option isn't filling a particularly critical need for the column indices of elements. But this route also might be the path of least resistance, so I'm OK if that's the way you want to go here.

Or there might be some other mechanism whereby the data structure passed to a problem poser is something of a slightly different type that definitely has a column for each element. (Perhaps this mild redundant-checking issue is even evidence that elements shouldn't store their column indices; conceptually, a column index is more of a property of an Element as a part of an Assembly, rather than intrinsically of the Element itself. So maybe instead the Assembly has some record of the column index for each Element, where the indices are plain ol' integers, and that record is provided to problem posers. But I am definitely not mandating a refactor of that scope as part of this PR.)

What doesn't seem good is writing out checks and manual panic messages in three places, for panics that can literally never occur because of the structure of the code. That redundant checking is what I mean by "evidence of misfactoring".

The check here and the potential panic if it fails (and analogous code in two other ProblemPoser implementations) is evidence of misfactoring: in fact, the ProblemPosers only ever get to execute when every element has just been assigned a column index. So there is no need to have such a check. I am open to any reasonable way to factor these checks out. One is to insist that all elements in an assembly always have a column index. As far as I can see, the only thing that None column indices are currently used for is a transient period between an Element being added to an Assembly and the next call to that Assembly being realized. That period could be eliminated by simply assigning each Element the next available column index and immediately calling realize. Or there could be a method (likely in engine?) that takes a realized assembly and a new completely unconstrained element and updates the realization to include the unconstrained element -- that should presumably be far less work than a full call to realize. As far as I can tell, the only consistency bit that needs to be kept up is the tangent space bases in the realization. But I suppose it should be possible to calculate the new tangent space basis matriceso directly in the case of a new unconstrained element added to a realized assembly. Or we can insist that elements always have an index by assigning indices on creation, and dropping the invariant that an element with an index has to have tangent motions in the corresponding column of the tangent space basis matrices. That might not be so terrible -- if an element's column is out of range for the tangent space basis matrices, that must mean it is totally unconstrained, and so we can actually deform it freely, and so just switch to simpler deformation code for that case. Another possibility is to initially use column index -1 for a new element, and not worry about any -1s in the problem posing because the code is structured such that they never occur when a problem is being posed; we only need to worry about the -1s when deforming, and handle the case there. Another possibility (if I understand Rust and the Option type properly) is to simply use unwrap in these posers. It nominally could panic, but we know it won't because problems are only posed just after every element has been assigned an index. This route seems the least satisfying because we still have all of the baggage of Option even though the other possibilities make it clear that Option isn't filling a particularly critical need for the column indices of elements. But this route also might be the path of least resistance, so I'm OK if that's the way you want to go here. Or there might be some other mechanism whereby the data structure passed to a problem poser is something of a slightly different type that definitely has a column for each element. (Perhaps this mild redundant-checking issue is even evidence that elements shouldn't store their column indices; conceptually, a column index is more of a property of an Element as a part of an Assembly, rather than intrinsically of the Element itself. So maybe instead the Assembly has some record of the column index for each Element, where the indices are plain ol' integers, and that record is provided to problem posers. But I am definitely not mandating a refactor of that scope as part of this PR.) What doesn't seem good is writing out checks and manual panic messages in three places, for panics that can literally never occur because of the structure of the code. That redundant checking is what I mean by "evidence of misfactoring".

I think there are two problems here—a narrow one that falls within the scope of this pull request, and a broader one that deserves its own pull request.

Narrow: semantics issue

The conditional panic doesn't do a great job of communicating that we expect every element to have a column index when pose is called. Unwrapping the Option, as you point out, has the semantics we want, and that's what we do on the main branch (commit b86f176). However, unwrapping doesn't provide an informative panic message, which would make it harder to debug an invariant violation.

I've switched to the expect method, which has the same semantics, and allows us to provide our own panic message (commit 5eeb093). I've also rewritten the panic messages in the style recommended by the expect documentation.

Broad: factoring issue

I agree that we should eventually rewrite the column index system to avoid these checks entirely. I think that deserves its own pull request, since this pull request didn't create and doesn't modify the column index system.

I think there are two problems here—a narrow one that falls within the scope of this pull request, and a broader one that deserves its own pull request. #### Narrow: semantics issue The conditional panic doesn't do a great job of communicating that we expect every element to have a column index when `pose` is called. Unwrapping the `Option`, as you point out, has the semantics we want, and that's what we do on the main branch (commit b86f176). However, unwrapping doesn't provide an informative panic message, which would make it harder to debug an invariant violation. I've switched to the [`expect`](https://doc.rust-lang.org/std/option/enum.Option.html#method.expect) method, which has the same semantics, and allows us to provide our own panic message (commit 5eeb093). I've also rewritten the panic messages in the [style recommended](https://doc.rust-lang.org/std/option/enum.Option.html#recommended-message-style) by the `expect` documentation. #### Broad: factoring issue I agree that we should eventually rewrite the column index system to avoid these checks entirely. I think that deserves its own pull request, since this pull request didn't create and doesn't modify the column index system.
format!("Element \"{}\" should be indexed before writing problem data", self.id).as_str()
);
problem.gram.push_sym(index, index, 1.0);
problem.guess.set_column(index, &self.representation.get_clone_untracked());
}
}
pub trait Regulator: ProblemPoser + OutlineItem {
fn subjects(&self) -> Vec<ElementKey>;
fn measurement(&self) -> ReadSignal<f64>;
fn set_point(&self) -> Signal<SpecifiedValue>;
// this method is used to responsively precondition the assembly for
// realization when the regulator becomes a constraint, or is edited while
// acting as a constraint. it should track the set point, do any desired
// preconditioning when the set point is present, and use its return value
// to report whether the set is present. the default implementation does no
// preconditioning
fn try_activate(&self, _assembly: &Assembly) -> bool {
self.set_point().with(|set_pt| set_pt.is_present())
}
}
pub struct InversiveDistanceRegulator {
pub subjects: [ElementKey; 2],
pub measurement: ReadSignal<f64>,
pub set_point: Signal<SpecifiedValue>
}
impl InversiveDistanceRegulator {
pub fn new(subjects: [ElementKey; 2], assembly: &Assembly) -> InversiveDistanceRegulator {
let measurement = assembly.elements.map(
move |elts| {
let representations = subjects.map(|subj| elts[subj].representation);
representations[0].with(|rep_0|
representations[1].with(|rep_1|
rep_0.dot(&(&*Q * rep_1))
)
)
}
);
let set_point = create_signal(SpecifiedValue::from_empty_spec());
InversiveDistanceRegulator { subjects, measurement, set_point }
}
}
impl Regulator for InversiveDistanceRegulator {
fn subjects(&self) -> Vec<ElementKey> {
self.subjects.into()
}
fn measurement(&self) -> ReadSignal<f64> {
self.measurement
}
fn set_point(&self) -> Signal<SpecifiedValue> {
self.set_point
}
}
impl ProblemPoser for InversiveDistanceRegulator {
fn pose(&self, problem: &mut ConstraintProblem, elts: &Slab<Element>) {
self.set_point.with_untracked(|set_pt| {
if let Some(val) = set_pt.value {
let [row, col] = self.subjects.map(
|subj| elts[subj].column_index.expect(
"Subjects should be indexed before inversive distance regulator writes problem data"
)
);
problem.gram.push_sym(row, col, val);
}
});
}
glen marked this conversation as resolved Outdated

At some point, it would be nice to have symbolic names for all of the components in some standard place and use them anywhere needed. I am not saying it has to be part of this PR. But it could be. Thoughts?

At some point, it would be nice to have symbolic names for all of the components in some standard place and use them anywhere needed. I am not saying it has to be part of this PR. But it could be. Thoughts?

Since you have to make the other 3 the same as this 3, likely you may as well do centralized symbolic names for all the components right at the moment, although anything that makes the two 3s the same thing and not magic numbers is OK for the moment if you prefer not to do a full central enumeration of components.

Since you have to make the _other_ 3 the same as this 3, likely you may as well do centralized symbolic names for all the components right at the moment, although anything that makes the two 3s the same thing and not magic numbers is OK for the moment if you prefer not to do a full central enumeration of components.

Done, by making CURVATURE_COMPONENT an associated constant of the Element structure (commit 7f21e7e). When we introduce separate structures for different kinds of elements, this constant will go in the Sphere structure.

Done, by making `CURVATURE_COMPONENT` an associated constant of the `Element` structure (commit 7f21e7e). When we introduce separate structures for different kinds of elements, this constant will go in the `Sphere` structure.

There's still pretty tight coupling between Elements and the Engine, because Elements are storing their own representations in Engine coordinates. So the Elements are not really abstract geometric elements, divorced from Engine details. Or put another way, an Engine change would entail deeper changes to the innards of Elements than might be ideal. But there's no real need to pre-worry about fixing that until we actually want to try another Engine. And the revised code does eliminate unlinked occurrences of the same magic number, so resolving this conversation.

There's still pretty tight coupling between Elements and the Engine, because Elements are storing their own representations in Engine coordinates. So the Elements are not really abstract geometric elements, divorced from Engine details. Or put another way, an Engine change would entail deeper changes to the innards of Elements than might be ideal. But there's no real need to pre-worry about fixing that until we actually want to try another Engine. And the revised code does eliminate unlinked occurrences of the same magic number, so resolving this conversation.
}
pub struct HalfCurvatureRegulator {
pub subject: ElementKey,
pub measurement: ReadSignal<f64>,
pub set_point: Signal<SpecifiedValue>
}
impl HalfCurvatureRegulator {
pub fn new(subject: ElementKey, assembly: &Assembly) -> HalfCurvatureRegulator {
let measurement = assembly.elements.map(
move |elts| elts[subject].representation.with(
|rep| rep[Element::CURVATURE_COMPONENT]
)
);
let set_point = create_signal(SpecifiedValue::from_empty_spec());
HalfCurvatureRegulator { subject, measurement, set_point }
}
}
impl Regulator for HalfCurvatureRegulator {
fn subjects(&self) -> Vec<ElementKey> {
vec![self.subject]
}
fn measurement(&self) -> ReadSignal<f64> {
self.measurement
}
fn set_point(&self) -> Signal<SpecifiedValue> {
self.set_point
}
fn try_activate(&self, assembly: &Assembly) -> bool {
match self.set_point.with(|set_pt| set_pt.value) {
Some(half_curv) => {
let representation = assembly.elements.with_untracked(
|elts| elts[self.subject].representation
);
representation.update(
|rep| change_half_curvature(rep, half_curv)
);
true
}
None => false
}
}
}
impl ProblemPoser for HalfCurvatureRegulator {
fn pose(&self, problem: &mut ConstraintProblem, elts: &Slab<Element>) {
self.set_point.with_untracked(|set_pt| {
if let Some(val) = set_pt.value {
let col = elts[self.subject].column_index.expect(
"Subject should be indexed before half-curvature regulator writes problem data"
);
problem.frozen.push(Element::CURVATURE_COMPONENT, col, val);
}
});
}
}
// the velocity is expressed in uniform coordinates
pub struct ElementMotion<'a> {
pub key: ElementKey,
@ -137,7 +287,7 @@ type AssemblyMotion<'a> = Vec<ElementMotion<'a>>;
pub struct Assembly {
// elements and regulators
pub elements: Signal<Slab<Element>>,
pub regulators: Signal<Slab<Regulator>>,
pub regulators: Signal<Slab<Rc<dyn Regulator>>>,
// solution variety tangent space. the basis vectors are stored in
// configuration matrix format, ordered according to the elements' column
@ -167,26 +317,33 @@ impl Assembly {
// --- inserting elements and regulators ---
// insert an element into the assembly without checking whether we already
// insert a sphere 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
fn insert_element_unchecked(&self, elt: Element) {
fn insert_sphere_unchecked(&self, elt: Element) -> ElementKey {
// insert the sphere
let id = elt.id.clone();
let key = self.elements.update(|elts| elts.insert(elt));
self.elements_by_id.update(|elts_by_id| elts_by_id.insert(id, key));
// regulate the sphere's curvature
self.insert_regulator(HalfCurvatureRegulator::new(key, &self));
key
}
pub fn try_insert_element(&self, elt: Element) -> bool {
pub fn try_insert_sphere(&self, elt: Element) -> Option<ElementKey> {
let can_insert = self.elements_by_id.with_untracked(
|elts_by_id| !elts_by_id.contains_key(&elt.id)
);
if can_insert {
self.insert_element_unchecked(elt);
Some(self.insert_sphere_unchecked(elt))
} else {
None
}
can_insert
}
pub fn insert_new_element(&self) {
pub fn insert_new_sphere(&self) {
// find the next unused identifier in the default sequence
let mut id_num = 1;
let mut id = format!("sphere{}", id_num);
@ -197,70 +354,69 @@ impl Assembly {
id = format!("sphere{}", id_num);
glen marked this conversation as resolved Outdated

I figured out what's bothering me about this code. It's connected with my feeling that these two methods (insert_new_product_regulator and insert_new_half_curvature_regulator) should either be the same or be more highly factored than they are now. It's that there's a de-localization of concerns at the moment. It appears that the definition of the measurement of a product regulator or a half-curvature regulator currently "lives" in these insertion methods. But such code should be part of the code for that flavor of regulator, not in an "outside" method like this insertion gadget. The current code structure suggests that it would be perfectly legitimate to insert a half-curvature regulator that measures the mean of the first three coordinates. It would not be.

In other words, we don't have ideal encapsulation of regulators at the moment. Please refactor; and hopefully, after refactoring, these two functions will be identical or very thin veneer over something identical. Thanks!

I figured out what's bothering me about this code. It's connected with my feeling that these two methods (insert_new_product_regulator and insert_new_half_curvature_regulator) should either be the same or be more highly factored than they are now. It's that there's a de-localization of concerns at the moment. It appears that the definition of the measurement of a product regulator or a half-curvature regulator currently "lives" in these insertion methods. But such code should be part of the code for that flavor of regulator, not in an "outside" method like this insertion gadget. The current code structure suggests that it would be perfectly legitimate to insert a half-curvature regulator that measures the mean of the first three coordinates. It would not be. In other words, we don't have ideal encapsulation of regulators at the moment. Please refactor; and hopefully, after refactoring, these two functions will be identical or very thin veneer over something identical. Thanks!

But such code should be part of the code for that flavor of regulator, not in an "outside" method like this insertion gadget.

I think that's a great idea. I propose moving the insertion code into a method of the Regulator trait, analogous to how the code that writes constraint problem data is kept in the pose method of the ProblemPoser trait. If you'd like, we could try to make a supertrait for anything that can be inserted into an assembly, which would be even more analogous to the ProblemPoser trait, but I don't think that's called for yet.

> But such code should be part of the code for that flavor of regulator, not in an "outside" method like this insertion gadget. I think that's a great idea. I propose moving the insertion code into a method of the `Regulator` trait, analogous to how the code that writes constraint problem data is kept in the `pose` method of the `ProblemPoser` trait. If you'd like, we could try to make a supertrait for anything that can be inserted into an assembly, which would be even more analogous to the `ProblemPoser` trait, but I don't think that's called for yet.

yes a method of regulator trait sounds good, no not necessary to generalize into an assembly-inserting supertrait, at least not until we have something else that would also use that supertrait, i.e. no need for premature abstraction. Thanks!

yes a method of regulator trait sounds good, no not necessary to generalize into an assembly-inserting supertrait, at least not until we have something else that would also use that supertrait, i.e. no need for premature abstraction. Thanks!

Done (commits ee8a01b52d9975). The code that creates and inserts regulators looks much better-organized now!

Done (commits ee8a01b – 52d9975). The code that creates and inserts regulators looks much better-organized now!
}
// create and insert a new element
self.insert_element_unchecked(
// create and insert a sphere
let _ = self.insert_sphere_unchecked(
Element::new(
id,
format!("Sphere {}", id_num),
[0.75_f32, 0.75_f32, 0.75_f32],
DVector::<f64>::from_column_slice(&[0.0, 0.0, 0.0, 0.5, -0.5])
sphere(0.0, 0.0, 0.0, 1.0)
)
);
}
fn insert_regulator(&self, regulator: Regulator) {
let subjects = regulator.subjects;
let key = self.regulators.update(|regs| regs.insert(regulator));
let subject_regulators = self.elements.with(
|elts| (elts[subjects.0].regulators, elts[subjects.1].regulators)
pub fn insert_regulator<T: Regulator + 'static>(&self, regulator: T) {
// add the regulator to the assembly's regulator list
let regulator_rc = Rc::new(regulator);
let key = self.regulators.update(
|regs| regs.insert(regulator_rc.clone())
);
subject_regulators.0.update(|regs| regs.insert(key));
subject_regulators.1.update(|regs| regs.insert(key));
}
pub fn insert_new_regulator(self, subjects: (ElementKey, ElementKey)) {
// create and insert a new regulator
let measurement = self.elements.map(
move |elts| {
let reps = (
elts[subjects.0].representation.get_clone(),
elts[subjects.1].representation.get_clone()
);
reps.0.dot(&(&*Q * reps.1))
// add the regulator to each subject's regulator list
let subjects = regulator_rc.subjects();
let subject_regulators: Vec<_> = self.elements.with_untracked(
|elts| subjects.into_iter().map(
|subj| elts[subj].regulators
).collect()
);
for regulators in subject_regulators {
regulators.update(|regs| regs.insert(key));
}
glen marked this conversation as resolved Outdated

Aha, this 3 is the same 3 as that other one I recently commented about. They need to be instances of the same symbol. Otherwise, they're just magic numbers that may leave a maintainer wondering whether and why they need to be the same. Thanks for fixing.

Aha, this 3 is the _same_ 3 as that other one I recently commented about. They need to be instances of the same symbol. Otherwise, they're just magic numbers that may leave a maintainer wondering whether and why they need to be the same. Thanks for fixing.
// update the realization when the regulator becomes a constraint, or is
// edited while acting as a constraint
let self_for_effect = self.clone();
create_effect(move || {
/* DEBUG */
// log the regulator update
console::log_1(&JsValue::from(
format!("Updated regulator with subjects {:?}", regulator_rc.subjects())
));
if regulator_rc.try_activate(&self_for_effect) {
glen marked this conversation as resolved Outdated

This function evokes a few questions:

A) Seems like there is some duplication at least of structure/behavior here with the insert_new_product_regulator; is there anything that can be profitable factored out (some common insert_regulator common functionality that both of these can use)?

B) Why do you have to do so much work updating the guess (I think that's what's going on) when you start regulating curvature, but it seems like you don't do much of anything when you start regulating an inversive distance?

C) This very much has the look of engine code that has leaked into the assembly: nitty gritty dealing with the internal coordinates of an element. Any chance it could be implemented in the engine and just called here, or perhaps even just deferred until a pre-processing step when the realization is called? Either way, it also might make the code between insert curvature regulator and insert product regulator more similar and easier to share between the two.

This function evokes a few questions: A) Seems like there is some duplication at least of structure/behavior here with the insert_new_product_regulator; is there anything that can be profitable factored out (some common insert_regulator common functionality that both of these can use)? B) Why do you have to do so much work updating the guess (I think that's what's going on) when you start regulating curvature, but it seems like you don't do much of anything when you start regulating an inversive distance? C) This very much has the look of engine code that has leaked into the assembly: nitty gritty dealing with the internal coordinates of an element. Any chance it could be implemented in the engine and just called here, or perhaps even just deferred until a pre-processing step when the realization is called? Either way, it also might make the code between insert curvature regulator and insert product regulator more similar and easier to share between the two.

A) Seems like there is some duplication at least of structure/behavior here with the insert_new_product_regulator; is there anything that can be profitable factored out (some common insert_regulator common functionality that both of these can use)?

I agree that the code is organized similarly, but I haven't found anything we can usefully factor out. For example, I don't see a way to put the shared organization in a shared constructor, because ProductRegulator and HalfCurvatureRegulator are different structures. The only thing that seems really straightforwardly shared is initializing the set point to

create_signal(SpecifiedValue::from_empty_spec()),

which seems like too little code to factor out profitably.

> A) Seems like there is some duplication at least of structure/behavior here with the insert_new_product_regulator; is there anything that can be profitable factored out (some common insert_regulator common functionality that both of these can use)? I agree that the code is organized similarly, but I haven't found anything we can usefully factor out. For example, I don't see a way to put the shared organization in a shared constructor, because `ProductRegulator` and `HalfCurvatureRegulator` are different structures. The only thing that seems really straightforwardly shared is initializing the set point to ```create_signal(SpecifiedValue::from_empty_spec())```, which seems like too little code to factor out profitably.

B) Why do you have to do so much work updating the guess (I think that's what's going on) when you start regulating curvature, but it seems like you don't do much of anything when you start regulating an inversive distance?

We don't have to be so careful about updating the guess. My first prototype of the curvature regulator didn't do anything special when it started regulating; the curvature component of the representation vector would just get overwritten with the desired value at the beginning of the engine's realization routine. That seemed to work fine in the examples I played with.

I felt compelled to be careful about updating the guess for three related reasons:

  1. Just overwriting the curvature component of the representation vector messes up the vector's normalization. I worried that starting from a non-normalized guess could make the engine more likely to stall, and could make the assembly change especially erratically when the curvature regulator is turned on.
  2. I could imagine using a curvature regulator to make extreme curvature changes. For example, someone who's doing conformal geometry on the boundary of hyperbolic 3-space might switch between the Poincaré model and the upper half-space model by switching the curvature of the boundary sphere between 1 and 0. I worried that if we just took whatever realization the engine happened to find, extreme curvature changes might throw the sphere far outside the field of view.
  3. For a moderately curved sphere with nothing constrained except its curvature, one might expect that changing the curvature shouldn't affect the center much. I wanted to demonstrate that we could meet this expectation, in case you thought it was important.

I'd be open to also choosing the initial guess more carefully when we start regulating an inversive distance. That seems tricker to me, though, and I haven't felt a need for it.

> B) Why do you have to do so much work updating the guess (I think that's what's going on) when you start regulating curvature, but it seems like you don't do much of anything when you start regulating an inversive distance? We don't *have* to be so careful about updating the guess. My first prototype of the curvature regulator didn't do anything special when it started regulating; the curvature component of the representation vector would just get overwritten with the desired value at the beginning of the engine's realization routine. That seemed to work fine in the examples I played with. I felt compelled to be careful about updating the guess for three related reasons: 1. Just overwriting the curvature component of the representation vector messes up the vector's normalization. I worried that starting from a non-normalized guess could make the engine more likely to stall, and could make the assembly change especially erratically when the curvature regulator is turned on. 2. I could imagine using a curvature regulator to make extreme curvature changes. For example, someone who's doing conformal geometry on the boundary of hyperbolic 3-space might switch between the Poincaré model and the upper half-space model by switching the curvature of the boundary sphere between 1 and 0. I worried that if we just took whatever realization the engine happened to find, extreme curvature changes might throw the sphere far outside the field of view. 3. For a moderately curved sphere with nothing constrained except its curvature, one might expect that changing the curvature shouldn't affect the center much. I wanted to demonstrate that we could meet this expectation, in case you thought it was important. I'd be open to also choosing the initial guess more carefully when we start regulating an inversive distance. That seems tricker to me, though, and I haven't felt a need for it.

C) This very much has the look of engine code that has leaked into the assembly: nitty gritty dealing with the internal coordinates of an element.

Yes, I think this code is on the border between the assembly module's responsibilities and the engine's responsibilities. I decided to put it in the assembly module because its core purpose is to decide how an otherwise unconstrained sphere should behave when you change its curvature. That seems more like user interaction than like constraint solving to me.

On the other hand, the engine does provide representation-specific implementations of other user-facing, mostly-representation-agnostic tasks, like creating a sphere from a center and a radius (sphere) or a direction, an offset, and a curvature (sphere_with_offset). If we can find a mostly-representation-agnostic way to describe how we update the curvature, maybe it would feel similar in spirit to those tasks.

> C) This very much has the look of engine code that has leaked into the assembly: nitty gritty dealing with the internal coordinates of an element. Yes, I think this code is on the border between the assembly module's responsibilities and the engine's responsibilities. I decided to put it in the assembly module because its core purpose is to decide how an otherwise unconstrained sphere should behave when you change its curvature. That seems more like user interaction than like constraint solving to me. On the other hand, the engine does provide representation-specific implementations of other user-facing, mostly-representation-agnostic tasks, like creating a sphere from a center and a radius (`sphere`) or a direction, an offset, and a curvature (`sphere_with_offset`). If we can find a mostly-representation-agnostic way to describe how we update the curvature, maybe it would feel similar in spirit to those tasks.

On the other hand, the engine does provide representation-specific implementations of other user-facing, mostly-representation-agnostic tasks […]. If we can find a mostly-representation-agnostic way to describe how we update the curvature, maybe it would feel similar in spirit to those tasks.

I've adopted this approach in commit 4654bf0 by moving the half-curvature change routine into the engine module. I haven't come up with a representation-agnostic description of how the sphere is supposed to change, so I'm just calling the routine change_half_curvature for now.

> On the other hand, the engine does provide representation-specific implementations of other user-facing, mostly-representation-agnostic tasks […]. If we can find a mostly-representation-agnostic way to describe how we update the curvature, maybe it would feel similar in spirit to those tasks. I've adopted this approach in commit 4654bf0 by moving the half-curvature change routine into the engine module. I haven't come up with a representation-agnostic description of how the sphere is supposed to change, so I'm just calling the routine `change_half_curvature` for now.
self_for_effect.realize();
}
);
let set_point = create_signal(SpecifiedValue::from_empty_spec());
self.insert_regulator(Regulator {
subjects: subjects,
measurement: measurement,
set_point: set_point
});
/* DEBUG */
// print an updated list of regulators
console::log_1(&JsValue::from("Regulators:"));
self.regulators.with(|regs| {
self.regulators.with_untracked(|regs| {
for (_, reg) in regs.into_iter() {
console::log_5(
&JsValue::from(" "),
&JsValue::from(reg.subjects.0),
&JsValue::from(reg.subjects.1),
&JsValue::from(":"),
&reg.set_point.with_untracked(
|set_pt| JsValue::from(set_pt.spec.as_str())
console::log_1(&JsValue::from(format!(
" {:?}: {}",
reg.subjects(),
reg.set_point().with_untracked(
|set_pt| {
let spec = &set_pt.spec;
if spec.is_empty() {
"__".to_string()
} else {
spec.clone()
}
}
)
);
}
});
// update the realization when the regulator becomes a constraint, or is
// edited while acting as a constraint
create_effect(move || {
console::log_1(&JsValue::from(
format!("Updated constraint with subjects ({}, {})", subjects.0, subjects.1)
));
if set_point.with(|set_pt| set_pt.is_present()) {
self.realize();
)));
}
});
}
@ -275,55 +431,39 @@ impl Assembly {
}
});
// 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();
// set up the constraint problem
let problem = self.elements.with_untracked(|elts| {
let mut problem = ConstraintProblem::new(elts.len());
for (_, elt) in elts {
elt.pose(&mut problem, elts);
}
self.regulators.with_untracked(|regs| {
for (_, reg) in regs {
reg.set_point.with_untracked(|set_pt| {
if let Some(val) = set_pt.value {
let subjects = reg.subjects;
let row = elts[subjects.0].column_index.unwrap();
let col = elts[subjects.1].column_index.unwrap();
gram_to_be.push_sym(row, col, val);
}
});
reg.pose(&mut problem, elts);
}
});
// set up the initial configuration matrix and the diagonal of the
// Gram matrix
let mut guess_to_be = DMatrix::<f64>::zeros(5, elts.len());
for (_, elt) in elts {
let index = elt.column_index.unwrap();
gram_to_be.push_sym(index, index, 1.0);
guess_to_be.set_column(index, &elt.representation.get_clone_untracked());
}
(gram_to_be, guess_to_be)
problem
});
/* DEBUG */
// log the Gram matrix
console::log_1(&JsValue::from("Gram matrix:"));
gram.log_to_console();
problem.gram.log_to_console();
glen marked this conversation as resolved Outdated

It's OK with me if you just call the interior variable problem even though it's being returned out to an exterior variable named problem -- I don't think it would reduce clarity. Not necessary to change, just wanted to make sure you knew that I don't have any issue with using the same name in a case like this. And I do think that concise variable names are definitely helpful.

It's OK with me if you just call the interior variable `problem` even though it's being returned out to an exterior variable named `problem` -- I don't think it would reduce clarity. Not necessary to change, just wanted to make sure you knew that I don't have any issue with using the same name in a case like this. And I do think that concise variable names are definitely helpful.

Good to know. For now, I think it's worth clearly distinguishing the interior "builder" variable and the exterior "storage" variable, because they have some differences in usage: for example, the builder variable is mutable, and the storage variable isn't. However, I'd be fine with switching to the convention of using the same name for both.

Good to know. For now, I think it's worth clearly distinguishing the interior "builder" variable and the exterior "storage" variable, because they have some differences in usage: for example, the builder variable is mutable, and the storage variable isn't. However, I'd be fine with switching to the convention of using the same name for both.

Well, we are either going to have the convention of using the same variable name inside and outside, or always distinguishing the names. We're not going to do it one way in some places and the other way in other places. I like the "same name" convention better -- more concise with no loss in clarity. So if you are not against that change in convention, please just go ahead and change this instance to conform. We can discuss (A,B,C) from the other thread in person.

Well, we are either going to have the convention of using the same variable name inside and outside, or always distinguishing the names. We're not going to do it one way in some places and the other way in other places. I like the "same name" convention better -- more concise with no loss in clarity. So if you are not against that change in convention, please just go ahead and change this instance to conform. We can discuss (A,B,C) from the other thread in person.

I've switched to the same-name convention in commit 955220c.

I've switched to the same-name convention in commit 955220c.

Great, I see that, and this code is fine. It occurs to me that the organization that would be used for this code in some languages would be to have the problem = ConstraintProblem::new(...) outside the closure in self.elements.with_untracked(...) and then have some sort of reference to the problem be visible inside that closure. I suppose that is possible in Rust as well. If you think it's cleaner/clearer to organize it that way, feel free to refactor this, but also feel free to leave it as is if you prefer. In the latter case, please just resolve this thread.

Great, I see that, and this code is fine. It occurs to me that the organization that would be used for this code in some languages would be to have the `problem = ConstraintProblem::new(...)` outside the closure in self.elements.with_untracked(...) and then have some sort of reference to the problem be visible inside that closure. I suppose that is possible in Rust as well. If you think it's cleaner/clearer to organize it that way, feel free to refactor this, but also feel free to leave it as is if you prefer. In the latter case, please just resolve this thread.

I chose this organization so that problem could be immutable in the scope where it's used, and only mutable in the scope where it's created.

I chose this organization so that `problem` could be immutable in the scope where it's used, and only mutable in the scope where it's created.
/* DEBUG */
// log the initial configuration matrix
console::log_1(&JsValue::from("Old configuration:"));
for j in 0..guess.nrows() {
for j in 0..problem.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());
for k in 0..problem.guess.ncols() {
row_str.push_str(format!(" {:>8.3}", problem.guess[(j, k)]).as_str());
}
console::log_1(&JsValue::from(row_str));
}
// look for a configuration with the given Gram matrix
let (config, tangent, success, history) = realize_gram(
&gram, guess, &[],
1.0e-12, 0.5, 0.9, 1.1, 200, 110
&problem, 1.0e-12, 0.5, 0.9, 1.1, 200, 110
);
/* DEBUG */
@ -458,4 +598,48 @@ impl Assembly {
// sync
self.realize();
}
}
#[cfg(test)]
mod tests {
use crate::engine;
use super::*;
#[test]
#[should_panic(expected = "Element \"sphere\" should be indexed before writing problem data")]
fn unindexed_element_test() {
let _ = create_root(|| {
Element::new(
"sphere".to_string(),
"Sphere".to_string(),
[1.0_f32, 1.0_f32, 1.0_f32],
engine::sphere(0.0, 0.0, 0.0, 1.0)
).pose(&mut ConstraintProblem::new(1), &Slab::new());
});
}
#[test]
#[should_panic(expected = "Subjects should be indexed before inversive distance regulator writes problem data")]
fn unindexed_subject_test_inversive_distance() {
let _ = create_root(|| {
let mut elts = Slab::new();
let subjects = [0, 1].map(|k| {
elts.insert(
Element::new(
format!("sphere{k}"),
format!("Sphere {k}"),
[1.0_f32, 1.0_f32, 1.0_f32],
engine::sphere(0.0, 0.0, 0.0, 1.0)
)
)
});
elts[subjects[0]].column_index = Some(0);
InversiveDistanceRegulator {
subjects: subjects,
measurement: create_memo(|| 0.0),
set_point: create_signal(SpecifiedValue::try_from("0.0".to_string()).unwrap())
}.pose(&mut ConstraintProblem::new(2), &elts);
});
}
}

View file

@ -35,9 +35,43 @@ pub fn sphere_with_offset(dir_x: f64, dir_y: f64, dir_z: f64, off: f64, curv: f6
])
}
// given a sphere's representation vector, change the sphere's half-curvature to
// `half-curv` and then restore normalization by contracting the representation
// vector toward the curvature axis
pub fn change_half_curvature(rep: &mut DVector<f64>, half_curv: f64) {
// set the sphere's half-curvature to the desired value
rep[3] = half_curv;
// restore normalization by contracting toward the curvature axis
const SIZE_THRESHOLD: f64 = 1e-9;
let half_q_lt = -2.0 * half_curv * rep[4];
let half_q_lt_sq = half_q_lt * half_q_lt;
let mut spatial = rep.fixed_rows_mut::<3>(0);
let q_sp = spatial.norm_squared();
if q_sp < SIZE_THRESHOLD && half_q_lt_sq < SIZE_THRESHOLD {
spatial.copy_from_slice(
&[0.0, 0.0, (1.0 - 2.0 * half_q_lt).sqrt()]
);
} else {
let scaling = half_q_lt + (q_sp + half_q_lt_sq).sqrt();
spatial.scale_mut(1.0 / scaling);
rep[4] /= scaling;
}
/* DEBUG */
// verify normalization
let rep_for_debug = rep.clone();
console::log_1(&JsValue::from(
format!(
"Sphere self-product after curvature change: {}",
rep_for_debug.dot(&(&*Q * &rep_for_debug))
)
));
}
// --- partial matrices ---
struct MatrixEntry {
pub struct MatrixEntry {
index: (usize, usize),
value: f64
}
@ -49,42 +83,72 @@ impl PartialMatrix {
PartialMatrix(Vec::<MatrixEntry>::new())
}
pub fn push_sym(&mut self, row: usize, col: usize, value: f64) {
pub fn push(&mut self, row: usize, col: usize, value: f64) {
let PartialMatrix(entries) = self;
entries.push(MatrixEntry { index: (row, col), value: value });
}
pub fn push_sym(&mut self, row: usize, col: usize, value: f64) {
self.push(row, col, value);
if row != col {
entries.push(MatrixEntry { index: (col, row), value: value });
self.push(col, row, 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()));
for &MatrixEntry { index: (row, col), value } in self {
console::log_1(&JsValue::from(
format!(" {} {} {}", row, col, value)
));
}
}
fn freeze(&self, a: &DMatrix<f64>) -> DMatrix<f64> {
let mut result = a.clone();
for &MatrixEntry { index, value } in self {
result[index] = value;
}
result
}
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];
for &MatrixEntry { index, .. } in self {
result[index] = a[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];
for &MatrixEntry { index, value } in self {
result[index] = value - rhs[index];
}
result
}
}
impl IntoIterator for PartialMatrix {
type Item = MatrixEntry;
type IntoIter = std::vec::IntoIter<Self::Item>;
fn into_iter(self) -> Self::IntoIter {
let PartialMatrix(entries) = self;
entries.into_iter()
}
}
impl<'a> IntoIterator for &'a PartialMatrix {
type Item = &'a MatrixEntry;
type IntoIter = std::slice::Iter<'a, MatrixEntry>;
fn into_iter(self) -> Self::IntoIter {
let PartialMatrix(entries) = self;
entries.into_iter()
}
}
// --- configuration subspaces ---
#[derive(Clone)]
@ -195,6 +259,34 @@ impl DescentHistory {
}
}
// --- constraint problems ---
pub struct ConstraintProblem {
pub gram: PartialMatrix,
pub frozen: PartialMatrix,
pub guess: DMatrix<f64>,
}
impl ConstraintProblem {
pub fn new(element_count: usize) -> ConstraintProblem {
const ELEMENT_DIM: usize = 5;
ConstraintProblem {
gram: PartialMatrix::new(),
frozen: PartialMatrix::new(),
guess: DMatrix::<f64>::zeros(ELEMENT_DIM, element_count)
}
}
#[cfg(feature = "dev")]
pub fn from_guess(guess_columns: &[DVector<f64>]) -> ConstraintProblem {
ConstraintProblem {
gram: PartialMatrix::new(),
frozen: PartialMatrix::new(),
guess: DMatrix::from_columns(guess_columns)
}
}
}
// --- gram matrix realization ---
// the Lorentz form
@ -286,12 +378,12 @@ fn seek_better_config(
None
}
// seek a matrix `config` for which `config' * Q * config` matches the partial
// matrix `gram`. use gradient descent starting from `guess`
// seek a matrix `config` that matches the partial matrix `problem.frozen` and
// has `config' * Q * config` matching the partial matrix `problem.gram`. start
// at `problem.guess`, set the frozen entries to their desired values, and then
// use a regularized Newton's method to seek the desired Gram matrix
pub fn realize_gram(
gram: &PartialMatrix,
guess: DMatrix<f64>,
frozen: &[(usize, usize)],
problem: &ConstraintProblem,
scaled_tol: f64,
min_efficiency: f64,
backoff: f64,
@ -299,6 +391,11 @@ pub fn realize_gram(
max_descent_steps: i32,
max_backoff_steps: i32
) -> (DMatrix<f64>, ConfigSubspace, bool, DescentHistory) {
// destructure the problem data
let ConstraintProblem {
gram, guess, frozen
} = problem;
// start the descent history
let mut history = DescentHistory::new();
@ -313,11 +410,11 @@ pub fn realize_gram(
// convert the frozen indices to stacked format
let frozen_stacked: Vec<usize> = frozen.into_iter().map(
|index| index.1*element_dim + index.0
|MatrixEntry { index: (row, col), .. }| col*element_dim + row
).collect();
// use Newton's method with backtracking and gradient descent backup
let mut state = SearchState::from_config(gram, guess);
// use a regularized Newton's method with backtracking
let mut state = SearchState::from_config(gram, frozen.freeze(guess));
let mut hess = DMatrix::zeros(element_dim, assembly_dim);
for _ in 0..max_descent_steps {
// find the negative gradient of the loss function
@ -415,7 +512,7 @@ pub fn realize_gram(
#[cfg(feature = "dev")]
pub mod examples {
use std::{array, f64::consts::PI};
use std::f64::consts::PI;
use super::*;
@ -428,35 +525,7 @@ pub mod examples {
// https://www.nippon.com/en/japan-topics/c12801/
//
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(
let mut problem = ConstraintProblem::from_guess(
[
sphere(0.0, 0.0, 0.0, 15.0),
sphere(0.0, 0.0, -9.0, 5.0),
@ -471,42 +540,45 @@ pub mod examples {
).collect::<Vec<_>>().as_slice()
);
for s in 0..9 {
// each sphere is represented by a spacelike vector
problem.gram.push_sym(s, s, 1.0);
// the circumscribing sphere is tangent to all of the other
// spheres, with matching orientation
if s > 0 {
problem.gram.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 {
problem.gram.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;
problem.gram.push_sym(s, s_next, -1.0);
}
}
// 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));
for k in 0..4 {
problem.frozen.push(3, k, problem.guess[(3, k)]);
}
realize_gram(
&gram, guess, &frozen,
scaled_tol, 0.5, 0.9, 1.1, 200, 110
)
realize_gram(&problem, scaled_tol, 0.5, 0.9, 1.1, 200, 110)
}
// set up a kaleidocycle, made of points with fixed distances between them,
// and find its tangent space
pub fn realize_kaleidocycle(scaled_tol: f64) -> (DMatrix<f64>, ConfigSubspace, bool, DescentHistory) {
const N_POINTS: usize = 12;
let gram = {
let mut gram_to_be = PartialMatrix::new();
for block in (0..N_POINTS).step_by(2) {
let block_next = (block + 2) % N_POINTS;
for j in 0..2 {
// diagonal and hinge edges
for k in j..2 {
gram_to_be.push_sym(block + j, block + k, if j == k { 0.0 } else { -0.5 });
}
// non-hinge edges
for k in 0..2 {
gram_to_be.push_sym(block + j, block_next + k, -0.625);
}
}
}
gram_to_be
};
let guess = {
const N_HINGES: usize = 6;
let guess_elts = (0..N_HINGES).step_by(2).flat_map(
const N_HINGES: usize = 6;
let mut problem = ConstraintProblem::from_guess(
(0..N_HINGES).step_by(2).flat_map(
|n| {
let ang_hor = (n as f64) * PI/3.0;
let ang_vert = ((n + 1) as f64) * PI/3.0;
@ -519,16 +591,30 @@ pub mod examples {
point(x_vert, y_vert, 0.5)
]
}
).collect::<Vec<_>>();
DMatrix::from_columns(&guess_elts)
};
).collect::<Vec<_>>().as_slice()
);
let frozen: [_; N_POINTS] = array::from_fn(|k| (3, k));
const N_POINTS: usize = 2 * N_HINGES;
for block in (0..N_POINTS).step_by(2) {
let block_next = (block + 2) % N_POINTS;
for j in 0..2 {
// diagonal and hinge edges
for k in j..2 {
problem.gram.push_sym(block + j, block + k, if j == k { 0.0 } else { -0.5 });
}
// non-hinge edges
for k in 0..2 {
problem.gram.push_sym(block + j, block_next + k, -0.625);
}
}
}
realize_gram(
&gram, guess, &frozen,
scaled_tol, 0.5, 0.9, 1.1, 200, 110
)
for k in 0..N_POINTS {
problem.frozen.push(3, k, problem.guess[(3, k)])
}
realize_gram(&problem, scaled_tol, 0.5, 0.9, 1.1, 200, 110)
}
}
@ -539,6 +625,25 @@ mod tests {
use super::{*, examples::*};
#[test]
fn freeze_test() {
let frozen = PartialMatrix(vec![
MatrixEntry { index: (0, 0), value: 14.0 },
MatrixEntry { index: (0, 2), value: 28.0 },
MatrixEntry { index: (1, 1), value: 42.0 },
MatrixEntry { index: (1, 2), value: 49.0 }
]);
let config = 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, &[
14.0, 2.0, 28.0,
4.0, 42.0, 49.0
]);
assert_eq!(frozen.freeze(&config), expected_result);
}
#[test]
fn sub_proj_test() {
let target = PartialMatrix(vec![
@ -560,18 +665,12 @@ mod tests {
#[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 }
});
}
let mut gram = PartialMatrix::new();
for j in 0..3 {
for k in 0..3 {
gram.push(j, k, if j == k { 1.0 } else { -1.0 });
}
entries
});
}
let config = {
let a = 0.75_f64.sqrt();
DMatrix::from_columns(&[
@ -584,37 +683,33 @@ mod tests {
assert!(state.loss.abs() < f64::EPSILON);
}
/* TO DO */
// at the frozen indices, the optimization steps should have exact zeros,
// and the realized configuration should match the initial guess
// and the realized configuration should have the desired values
#[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(&[
let mut problem = ConstraintProblem::from_guess(&[
point(0.0, 0.0, 2.0),
sphere(0.0, 0.0, 0.0, 1.0)
sphere(0.0, 0.0, 0.0, 0.95)
]);
let frozen = [(3, 0), (3, 1)];
println!();
for j in 0..2 {
for k in j..2 {
problem.gram.push_sym(j, k, if (j, k) == (1, 1) { 1.0 } else { 0.0 });
}
}
problem.frozen.push(3, 0, problem.guess[(3, 0)]);
problem.frozen.push(3, 1, 0.5);
let (config, _, success, history) = realize_gram(
&gram, guess.clone(), &frozen,
1.0e-12, 0.5, 0.9, 1.1, 200, 110
&problem, 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 {
for &MatrixEntry { index, .. } in &problem.frozen {
assert_eq!(base_step[index], 0.0);
}
}
for index in frozen {
assert_eq!(config[index], guess[index]);
for MatrixEntry { index, value } in problem.frozen {
assert_eq!(config[index], value);
}
}
@ -635,34 +730,32 @@ mod tests {
#[test]
fn tangent_test_three_spheres() {
const SCALED_TOL: f64 = 1.0e-12;
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(&[
const ELEMENT_DIM: usize = 5;
let mut problem = ConstraintProblem::from_guess(&[
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));
for j in 0..3 {
for k in j..3 {
problem.gram.push_sym(j, k, if j == k { 1.0 } else { -1.0 });
}
}
for n in 0..ELEMENT_DIM {
problem.frozen.push(n, 0, problem.guess[(n, 0)]);
}
let (config, tangent, success, history) = realize_gram(
&gram, guess.clone(), &frozen,
SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
&problem, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
);
assert_eq!(config, guess);
assert_eq!(config, problem.guess);
assert_eq!(success, true);
assert_eq!(history.scaled_loss.len(), 1);
// list some motions that should form a basis for the tangent space of
// the solution variety
const UNIFORM_DIM: usize = 4;
let element_dim = guess.nrows();
let assembly_dim = guess.ncols();
let element_dim = problem.guess.nrows();
let assembly_dim = problem.guess.ncols();
let tangent_motions_unif = vec![
basis_matrix((0, 1), UNIFORM_DIM, assembly_dim),
basis_matrix((1, 1), UNIFORM_DIM, assembly_dim),
@ -805,22 +898,17 @@ mod tests {
fn proj_equivar_test() {
// find a pair of spheres that meet at 120°
const SCALED_TOL: f64 = 1.0e-12;
let gram = {
let mut gram_to_be = PartialMatrix::new();
gram_to_be.push_sym(0, 0, 1.0);
gram_to_be.push_sym(1, 1, 1.0);
gram_to_be.push_sym(0, 1, 0.5);
gram_to_be
};
let guess_orig = DMatrix::from_columns(&[
let mut problem_orig = ConstraintProblem::from_guess(&[
sphere(0.0, 0.0, 0.5, 1.0),
sphere(0.0, 0.0, -0.5, 1.0)
]);
problem_orig.gram.push_sym(0, 0, 1.0);
problem_orig.gram.push_sym(1, 1, 1.0);
problem_orig.gram.push_sym(0, 1, 0.5);
let (config_orig, tangent_orig, success_orig, history_orig) = realize_gram(
&gram, guess_orig.clone(), &[],
SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
&problem_orig, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
);
assert_eq!(config_orig, guess_orig);
assert_eq!(config_orig, problem_orig.guess);
assert_eq!(success_orig, true);
assert_eq!(history_orig.scaled_loss.len(), 1);
@ -833,11 +921,15 @@ mod tests {
sphere(-a, 0.0, 7.0 - a, 1.0)
])
};
let problem_tfm = ConstraintProblem {
gram: problem_orig.gram,
guess: guess_tfm,
frozen: problem_orig.frozen
};
let (config_tfm, tangent_tfm, success_tfm, history_tfm) = realize_gram(
&gram, guess_tfm.clone(), &[],
SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
&problem_tfm, SCALED_TOL, 0.5, 0.9, 1.1, 200, 110
);
assert_eq!(config_tfm, guess_tfm);
assert_eq!(config_tfm, problem_tfm.guess);
assert_eq!(success_tfm, true);
assert_eq!(history_tfm.scaled_loss.len(), 1);
@ -869,7 +961,7 @@ mod tests {
// the comparison tolerance because the transformation seems to
// introduce some numerical error
const SCALED_TOL_TFM: f64 = 1.0e-9;
let tol_sq = ((guess_orig.nrows() * guess_orig.ncols()) as f64) * SCALED_TOL_TFM * SCALED_TOL_TFM;
let tol_sq = ((problem_orig.guess.nrows() * problem_orig.guess.ncols()) as f64) * SCALED_TOL_TFM * SCALED_TOL_TFM;
assert!((motion_proj_tfm - motion_tfm_proj).norm_squared() < tol_sq);
}
}

View file

@ -1,4 +1,5 @@
use itertools::Itertools;
use std::rc::Rc;
use sycamore::prelude::*;
use web_sys::{
KeyboardEvent,
@ -9,24 +10,38 @@ use web_sys::{
use crate::{
AppState,
assembly,
assembly::{ElementKey, Regulator, RegulatorKey},
assembly::{
ElementKey,
HalfCurvatureRegulator,
InversiveDistanceRegulator,
Regulator,
RegulatorKey
},
specified::SpecifiedValue
};
// an editable view of a regulator
#[component(inline_props)]
fn RegulatorInput(regulator: Regulator) -> View {
fn RegulatorInput(regulator: Rc<dyn Regulator>) -> View {
// get the regulator's measurement and set point signals
let measurement = regulator.measurement();
let set_point = regulator.set_point();
// the `valid` signal tracks whether the last entered value is a valid set
// point specification
let valid = create_signal(true);
// the `value` signal holds the current set point specification
let value = create_signal(
regulator.set_point.with_untracked(|set_pt| set_pt.spec.clone())
set_point.with_untracked(|set_pt| set_pt.spec.clone())
);
// this closure resets the input value to the regulator's set point
// specification
// this `reset_value` closure resets the input value to the regulator's set
// point specification
let reset_value = move || {
batch(|| {
valid.set(true);
value.set(regulator.set_point.with(|set_pt| set_pt.spec.clone()));
value.set(set_point.with(|set_pt| set_pt.spec.clone()));
})
};
@ -39,7 +54,7 @@ fn RegulatorInput(regulator: Regulator) -> View {
r#type="text",
class=move || {
if valid.get() {
regulator.set_point.with(|set_pt| {
set_point.with(|set_pt| {
if set_pt.is_present() {
"regulator-input constraint"
} else {
@ -50,13 +65,13 @@ fn RegulatorInput(regulator: Regulator) -> View {
"regulator-input invalid"
}
},
placeholder=regulator.measurement.with(|result| result.to_string()),
placeholder=measurement.with(|result| result.to_string()),
bind:value=value,
on:change=move |_| {
valid.set(
match SpecifiedValue::try_from(value.get_clone_untracked()) {
Ok(set_pt) => {
regulator.set_point.set(set_pt);
set_point.set(set_pt);
true
}
Err(_) => false
@ -75,26 +90,53 @@ fn RegulatorInput(regulator: Regulator) -> View {
}
}
pub trait OutlineItem {
fn outline_item(self: Rc<Self>, element_key: ElementKey) -> View;
}
impl OutlineItem for InversiveDistanceRegulator {
fn outline_item(self: Rc<Self>, element_key: ElementKey) -> View {
let state = use_context::<AppState>();
let other_subject = if self.subjects[0] == element_key {
self.subjects[1]
} else {
self.subjects[0]
};
let other_subject_label = state.assembly.elements.with(
|elts| elts[other_subject].label.clone()
);
view! {
li(class="regulator") {
div(class="regulator-label") { (other_subject_label) }
div(class="regulator-type") { "Inversive distance" }
RegulatorInput(regulator=self)
div(class="status")
}
}
}
}
impl OutlineItem for HalfCurvatureRegulator {
fn outline_item(self: Rc<Self>, _element_key: ElementKey) -> View {
view! {
li(class="regulator") {
div(class="regulator-label") // for spacing
div(class="regulator-type") { "Half-curvature" }
RegulatorInput(regulator=self)
div(class="status")
}
}
}
}
// a list item that shows a regulator in an outline view of an element
#[component(inline_props)]
fn RegulatorOutlineItem(regulator_key: RegulatorKey, element_key: ElementKey) -> View {
let state = use_context::<AppState>();
let assembly = &state.assembly;
let regulator = assembly.regulators.with(|regs| regs[regulator_key]);
let other_subject = if regulator.subjects.0 == element_key {
regulator.subjects.1
} else {
regulator.subjects.0
};
let other_subject_label = assembly.elements.with(|elts| elts[other_subject].label.clone());
view! {
li(class="regulator") {
div(class="regulator-label") { (other_subject_label) }
div(class="regulator-type") { "Inversive distance" }
RegulatorInput(regulator=regulator)
div(class="status")
}
}
let regulator = state.assembly.regulators.with(
|regs| regs[regulator_key].clone()
);
regulator.outline_item(element_key)
}
// a list item that shows an element in an outline view of an assembly
@ -117,7 +159,15 @@ fn ElementOutlineItem(key: ElementKey, element: assembly::Element) -> View {
};
let regulated = element.regulators.map(|regs| regs.len() > 0);
let regulator_list = element.regulators.map(
|regs| regs.clone().into_iter().collect()
move |elt_reg_keys| elt_reg_keys
.clone()
.into_iter()
.sorted_by_key(
|&reg_key| state.assembly.regulators.with(
|regs| regs[reg_key].subjects().len()
)
)
.collect()
);
let details_node = create_node_ref();
view! {