Simplify the realization triggering system #105

Merged
glen merged 6 commits from Vectornaut/dyna3:reactive-realization-cleanup into main 2025-07-31 22:21:34 +00:00
4 changed files with 36 additions and 96 deletions
Showing only changes of commit c73008d702 - Show all commits

View file

@ -16,7 +16,6 @@ use crate::{
components::{display::DisplayItem, outline::OutlineItem},
engine::{
Q,
change_half_curvature,
local_unif_to_std,
point,
project_point_to_normalized,
@ -358,16 +357,6 @@ pub trait Regulator: Serial + ProblemPoser + OutlineItem {
fn subjects(&self) -> Vec<Rc<dyn Element>>;
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) -> bool {
self.set_point().with(|set_pt| set_pt.is_present())
}
}
impl Hash for dyn Regulator {
@ -488,18 +477,6 @@ impl Regulator for HalfCurvatureRegulator {
fn set_point(&self) -> Signal<SpecifiedValue> {
self.set_point
}
fn try_activate(&self) -> bool {
match self.set_point.with(|set_pt| set_pt.value) {
Some(half_curv) => {
self.subject.representation().update(
|rep| change_half_curvature(rep, half_curv)
);
true
}
None => false
}
}
}
impl Serial for HalfCurvatureRegulator {
@ -552,8 +529,7 @@ pub struct Assembly {
pub elements_by_id: Signal<BTreeMap<String, Rc<dyn Element>>>,
// realization control
pub keep_realized: Signal<bool>,
pub needs_realization: Signal<bool>,
pub realization_trigger: Signal<()>,
// realization diagnostics
pub realization_status: Signal<Result<(), String>>,
@ -568,21 +544,23 @@ impl Assembly {
regulators: create_signal(BTreeSet::new()),
tangent: create_signal(ConfigSubspace::zero(0)),
elements_by_id: create_signal(BTreeMap::default()),
keep_realized: create_signal(true),
needs_realization: create_signal(false),
realization_trigger: create_signal(()),
realization_status: create_signal(Ok(())),
descent_history: create_signal(DescentHistory::new())
};
// realize the assembly whenever it becomes simultaneously true that
// we're trying to keep it realized and it needs realization
// realize the assembly whenever the element list, the regulator list,
// a regulator's set point, or the realization trigger is updated
let assembly_for_effect = assembly.clone();
create_effect(move || {
let should_realize = assembly_for_effect.keep_realized.get()
&& assembly_for_effect.needs_realization.get();
if should_realize {
assembly_for_effect.realize();
}
assembly_for_effect.elements.track();
assembly_for_effect.regulators.with(
|regs| for reg in regs {
reg.set_point().track();
}
);
assembly_for_effect.realization_trigger.track();
assembly_for_effect.realize();
});
assembly
@ -646,19 +624,6 @@ impl Assembly {
regulators.update(|regs| regs.insert(regulator.clone()));
}
// request a 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!("Updated regulator with subjects {:?}", regulator.subjects());
if regulator.try_activate() {
self_for_effect.needs_realization.set(true);
}
});
/* DEBUG */
// print an updated list of regulators
console_log!("Regulators:");
@ -726,8 +691,10 @@ impl Assembly {
} else {
console_log!("✅️ Target accuracy achieved!");
}
console_log!("Steps: {}", history.scaled_loss.len() - 1);
console_log!("Loss: {}", history.scaled_loss.last().unwrap());
if history.scaled_loss.len() > 0 {
console_log!("Steps: {}", history.scaled_loss.len() - 1);
console_log!("Loss: {}", history.scaled_loss.last().unwrap());
}
// report the loss history
self.descent_history.set(history);
@ -750,9 +717,6 @@ impl Assembly {
// save the tangent space
self.tangent.set_silent(tangent);
// clear the realization request flag
self.needs_realization.set(false);
},
Err(message) => {
// report the realization status. the `Err(message)` we're
@ -848,10 +812,10 @@ impl Assembly {
});
}
// request a realization to bring the configuration back onto the
// trigger a realization to bring the configuration back onto the
// solution variety. this also gets the elements' column indices and the
// saved tangent space back in sync
self.needs_realization.set(true);
self.realization_trigger.set(());
}
}

View file

@ -14,7 +14,12 @@ pub fn AddRemove() -> View {
button(
on:click=|_| {
let state = use_context::<AppState>();
state.assembly.insert_element_default::<Sphere>();
batch(|| {
// this call is batched to avoid redundant realizations.
// it updates the element list and the regulator list,
// which are both tracked by the realization effect
state.assembly.insert_element_default::<Sphere>();
});
}
) { "Add sphere" }
button(

View file

@ -900,9 +900,6 @@ pub fn TestAssemblyChooser() -> View {
let state = use_context::<AppState>();
let assembly = &state.assembly;
// pause realization
assembly.keep_realized.set(false);
// clear state
assembly.regulators.update(|regs| regs.clear());
assembly.elements.update(|elts| elts.clear());
@ -923,9 +920,6 @@ pub fn TestAssemblyChooser() -> View {
"irisawa-hexlet" => load_irisawa_hexlet_assemb(assembly),
_ => ()
};
// resume realization
assembly.keep_realized.set(true);
});
});

View file

@ -50,40 +50,6 @@ pub fn project_point_to_normalized(rep: &mut DVector<f64>) {
rep.scale_mut(0.5 / rep[3]);
}
// 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 ---
pub struct MatrixEntry {
@ -425,9 +391,20 @@ pub fn realize_gram(
// start the descent history
let mut history = DescentHistory::new();
// handle the empty-assembly case
let assembly_dim = guess.ncols();
if assembly_dim == 0 {
let result = Ok(
ConfigNeighborhood {
config: guess.clone(),
nbhd: ConfigSubspace::zero(0)
}
);
return Realization { result, history }
}
// find the dimension of the search space
glen marked this conversation as resolved Outdated

(a) Why do we now need an empty-assembly special case, when we didn't before? Or is it just for efficiency?

(b) Could this be expanded relatively easily to the more common no-set-regulators case?

(a) Why do we now need an empty-assembly special case, when we didn't before? Or is it just for efficiency? (b) Could this be expanded relatively easily to the more common no-set-regulators case?

(a) Why do we now need an empty-assembly special case, when we didn't before? Or is it just for efficiency?

The realization triggering system on the main branch never tries to realize an empty assembly, so we've never had to worry about the empty-assembly case.

As I mentioned in the pull request description, the realization triggering system on the incoming branch allows some unnecessary realizations in exchange for simplicity. One of those unnecessary realizations happens when we construct an empty assembly using the constructor Assembly::new.

As it turns out, the realization function on the main branch panics when you try to realize an empty assembly. This is because an empty assembly has an empty (zero by zero) Hessian, and the DMatrix::from_columns function, which we use to build the Hessian, can't handle an empty list of columns.

(b) Could this be expanded relatively easily to the more common no-set-regulators case?

Sure. We could give PartialMatrix an is_empty method and then expand the early return condition to:

assembly_dim == 0 || (gram.is_empty() && frozen.is_empty())
> (a) Why do we now need an empty-assembly special case, when we didn't before? Or is it just for efficiency? The realization triggering system on the main branch never tries to realize an empty assembly, so we've never had to worry about the empty-assembly case. As I mentioned in the pull request description, the realization triggering system on the incoming branch allows some unnecessary realizations in exchange for simplicity. One of those unnecessary realizations happens when we construct an empty assembly using the constructor `Assembly::new`. As it turns out, the realization function on the main branch panics when you try to realize an empty assembly. This is because an empty assembly has an empty (zero by zero) Hessian, and the `DMatrix::from_columns` function, which we use to build the Hessian, can't handle an empty list of columns. > (b) Could this be expanded relatively easily to the more common no-set-regulators case? Sure. We could give `PartialMatrix` an `is_empty` method and then expand the early return condition to: ```rust assembly_dim == 0 || (gram.is_empty() && frozen.is_empty()) ```

Any reason not to go ahead and do that in this PR? Seems like more bang for not many more bucks -- as long as we are checking for one kind of trivial case, why not exclude many more trivial cases? Seems to me like there's no downside...

Any reason not to go ahead and do that in this PR? Seems like more bang for not many more bucks -- as long as we are checking for one kind of trivial case, why not exclude many more trivial cases? Seems to me like there's no downside...

Any reason not to go ahead and do that in this PR?

On second thought, we'd actually have to use a separate conditional for the unconstrained case, because the tangent space would be different. When there are no elements, we return the zero subspace as the tangent space; when there are elements but no constraints, we want to return the whole configuration vector space as the tangent space.

I'd recommend leaving realize_gram mostly as it is. When there are no constraints, realize_gram still exits before doing a Cholesky decomposition, which is the most costly part of the optimization step: the loop breaks halfway through the first iteration, in the block with the comment "stop if the loss is tolerably low." After the break, we use ConfigSubspace::symmetric_kernel to find the tangent space. This costs more computationally than hard-coding the tangent space for the unconstrained case, but I think it makes maintenance easier.

Note that to get the spectrum for the diagnostics panel, realize_gram finds the eigenvalues of the Hessian on every step. I'm not factoring this into the cost considerations above, because I'd only want to do this in development builds.

> Any reason not to go ahead and do that in this PR? On second thought, we'd actually have to use a separate conditional for the unconstrained case, because the tangent space would be different. When there are no elements, we return the zero subspace as the tangent space; when there are elements but no constraints, we want to return the whole configuration vector space as the tangent space. I'd recommend leaving `realize_gram` mostly as it is. When there are no constraints, `realize_gram` still exits before doing a Cholesky decomposition, which is the most costly part of the optimization step: the loop breaks halfway through the first iteration, in the block with the comment "stop if the loss is tolerably low." After the break, we use `ConfigSubspace::symmetric_kernel` to find the tangent space. This costs more computationally than hard-coding the tangent space for the unconstrained case, but I think it makes maintenance easier. Note that to get the spectrum for the diagnostics panel, `realize_gram` finds the eigenvalues of the Hessian on every step. I'm not factoring this into the cost considerations above, because I'd only want to do this in development builds.

@Vectornaut wrote in #105 (comment):

On second thought, we'd actually have to use a separate conditional for the unconstrained case, because the tangent space would be different. When there are no elements, we return the zero subspace as the tangent space; when there are elements but no constraints, we want to return the whole configuration vector space as the tangent space.

How do those differ? When there are no elements, the whole configuration vector space is the zero space...

@Vectornaut wrote in https://code.studioinfinity.org/StudioInfinity/dyna3/pulls/105#issuecomment-3114: > On second thought, we'd actually have to use a separate conditional for the unconstrained case, because the tangent space would be different. When there are no elements, we return the zero subspace as the tangent space; when there are elements but no constraints, we want to return the whole configuration vector space as the tangent space. How do those differ? When there are no elements, the whole configuration vector space is the zero space...

How do those differ? When there are no elements, the whole configuration vector space is the zero space...

Good point: in principle, a ConfigSubspace variant representing a full subspace could be used for both empty assemblies and unconstrained assemblies. I've prototyped a full subspace representation and saved the code locally in case we want it later. In the process, however, I realized that we'll never actually encounter an unconstrained assembly, because every element comes with a normalization constraint. I'm now strongly against adding special handling for unconstrained assemblies, because during normal operation, this code will never run.

> How do those differ? When there are no elements, the whole configuration vector space is the zero space... Good point: in principle, a `ConfigSubspace` variant representing a full subspace could be used for both empty assemblies and unconstrained assemblies. I've prototyped a full subspace representation and saved the code locally in case we want it later. In the process, however, I realized that we'll never actually encounter an unconstrained assembly, because every element comes with a normalization constraint. I'm now strongly against adding special handling for unconstrained assemblies, because during normal operation, this code will never run.

OK, yes, I had forgotten that our coordinates are themselves underdetermined and that we use essentially the same mechanism to keep the coordinates for a given item in the form we want as to place an actual Euclidean-geometric constraint. So I agree, there's no point in any code for unconstrained assemblies other than the totally empty assembly. So yes go ahead and leave the special case for empty assemblies; please just add a comment as to why it must be handled specially. Thanks!

OK, yes, I had forgotten that our coordinates are themselves underdetermined and that we use essentially the same mechanism to keep the coordinates for a given item in the form we want as to place an actual Euclidean-geometric constraint. So I agree, there's no point in any code for unconstrained assemblies other than the totally empty assembly. So yes go ahead and leave the special case for empty assemblies; please just add a comment as to why it must be handled specially. Thanks!

So yes go ahead and leave the special case for empty assemblies; please just add a comment as to why it must be handled specially. Thanks!

Done in commit 779c026!

> So yes go ahead and leave the special case for empty assemblies; please just add a comment as to why it must be handled specially. Thanks! Done in commit 779c026!
let element_dim = guess.nrows();
let assembly_dim = guess.ncols();
let total_dim = element_dim * assembly_dim;
// scale the tolerance