dyna3/app-proto/src/assembly.rs

use nalgebra::{DMatrix, DVector, DVectorView};
use std::{
    cell::Cell,
    cmp::Ordering,
    collections::{BTreeMap, BTreeSet},
    fmt,
    fmt::{Debug, Formatter},
    hash::{Hash, Hasher},
    rc::Rc,
    sync::{atomic, atomic::AtomicU64},
};
use sycamore::prelude::*;
use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */

use crate::{
    components::{display::DisplayItem, outline::OutlineItem},
    engine::{
        Q,
        local_unif_to_std,
        point,
        project_point_to_normalized,
        project_sphere_to_normalized,
        realize_gram,
        sphere,
        ConfigNeighborhood,
        ConfigSubspace,
        ConstraintProblem,
        DescentHistory,
        Realization,
    },
    specified::SpecifiedValue,
};

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 item has a key that identifies it within its assembly and
// each assembly has a key that identifies it within the sesssion
static NEXT_SERIAL: AtomicU64 = AtomicU64::new(0);

pub trait Serial {
    // a serial number that uniquely identifies this element
    fn serial(&self) -> u64;
    
    // take the next serial number, panicking if that was the last one left
    fn next_serial() -> u64 where Self: Sized {
        // 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
        //
        NEXT_SERIAL.fetch_update(
            atomic::Ordering::SeqCst, atomic::Ordering::SeqCst,
            |serial| serial.checked_add(1)
        ).expect("Out of serial numbers for elements")
    }
}

impl Hash for dyn Serial {
    fn hash<H: Hasher>(&self, state: &mut H) {
        self.serial().hash(state)
    }
}

impl PartialEq for dyn Serial {
    fn eq(&self, other: &Self) -> bool {
        self.serial() == other.serial()
    }
}

impl Eq for dyn Serial {}

impl PartialOrd for dyn Serial {
    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
        Some(self.cmp(other))
    }
}

impl Ord for dyn Serial {
    fn cmp(&self, other: &Self) -> Ordering {
        self.serial().cmp(&other.serial())
    }
}

pub trait ProblemPoser {
    fn pose(&self, problem: &mut ConstraintProblem);
}

pub trait Element: Serial + ProblemPoser + DisplayItem {
    // the default identifier for an element of this type
    fn default_id() -> String where Self: Sized;
    
    // the default example of an element of this type
    fn default(id: String, id_num: u64) -> Self where Self: Sized;
    
    // the default regulators that come with this element
    fn default_regulators(self: Rc<Self>) -> Vec<Rc<dyn Regulator>> {
        Vec::new()
    }
    
    fn id(&self) -> &String;
    fn label(&self) -> &String;
    fn representation(&self) -> Signal<DVector<f64>>;
    fn ghost(&self) -> Signal<bool>;
    
    // the regulators the element is subject to. the assembly that owns the
    // element is responsible for keeping this set up to date
    fn regulators(&self) -> Signal<BTreeSet<Rc<dyn Regulator>>>;
    
    // project a representation vector for this kind of element onto its
    // normalization variety
    fn project_to_normalized(&self, rep: &mut DVector<f64>);
    
    // the configuration matrix column index that was assigned to the element
    // last time the assembly was realized, or `None` if the element has never
    // been through a realization
    fn column_index(&self) -> Option<usize>;
    
    // assign the element a configuration matrix column index. this method must
    // be used carefully to preserve invariant (1), described in the comment on
    // the `tangent` field of the `Assembly` structure
    fn set_column_index(&self, index: usize);
}

impl Debug for dyn Element {
    fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), fmt::Error> {
        self.id().fmt(f)
    }
}

impl Hash for dyn Element {
    fn hash<H: Hasher>(&self, state: &mut H) {
        <dyn Serial>::hash(self, state)
    }
}

impl PartialEq for dyn Element {
    fn eq(&self, other: &Self) -> bool {
        <dyn Serial>::eq(self, other)
    }
}

impl Eq for dyn Element {}

impl PartialOrd for dyn Element {
    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
        <dyn Serial>::partial_cmp(self, other)
    }
}

impl Ord for dyn Element {
    fn cmp(&self, other: &Self) -> Ordering {
        <dyn Serial>::cmp(self, other)
    }
}

pub struct Sphere {
    pub id: String,
    pub label: String,
    pub color: ElementColor,
    pub representation: Signal<DVector<f64>>,
    pub ghost: Signal<bool>,
    pub regulators: Signal<BTreeSet<Rc<dyn Regulator>>>,
    serial: u64,
    column_index: Cell<Option<usize>>,
}

impl Sphere {
    const CURVATURE_COMPONENT: usize = 3;
    
    pub fn new(
        id: String,
        label: String,
        color: ElementColor,
        representation: DVector<f64>,
    ) -> Self {
        Self {
            id,
            label,
            color,
            representation: create_signal(representation),
            ghost: create_signal(false),
            regulators: create_signal(BTreeSet::new()),
            serial: Self::next_serial(),
            column_index: None.into(),
        }
    }
}

impl Element for Sphere {
    fn default_id() -> String {
        "sphere".to_string()
    }
    
    fn default(id: String, id_num: u64) -> Self {
        Self::new(
            id,
            format!("Sphere {id_num}"),
            [0.75_f32, 0.75_f32, 0.75_f32],
            sphere(0.0, 0.0, 0.0, 1.0),
        )
    }
    
    fn default_regulators(self: Rc<Self>) -> Vec<Rc<dyn Regulator>> {
        vec![Rc::new(HalfCurvatureRegulator::new(self))]
    }
    
    fn id(&self) -> &String {
        &self.id
    }
    
    fn label(&self) -> &String {
        &self.label
    }
    
    fn representation(&self) -> Signal<DVector<f64>> {
        self.representation
    }
    
    fn ghost(&self) -> Signal<bool> {
        self.ghost
    }
    
    fn regulators(&self) -> Signal<BTreeSet<Rc<dyn Regulator>>> {
        self.regulators
    }
    
    fn project_to_normalized(&self, rep: &mut DVector<f64>) {
        project_sphere_to_normalized(rep);
    }
    
    fn column_index(&self) -> Option<usize> {
        self.column_index.get()
    }
    
    fn set_column_index(&self, index: usize) {
        self.column_index.set(Some(index));
    }
}

impl Serial for Sphere {
    fn serial(&self) -> u64 {
        self.serial
    }
}

impl ProblemPoser for Sphere {
    fn pose(&self, problem: &mut ConstraintProblem) {
        let index = self.column_index().expect(
            format!("Sphere \"{}\" 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 struct Point {
    pub id: String,
    pub label: String,
    pub color: ElementColor,
    pub representation: Signal<DVector<f64>>,
    pub ghost: Signal<bool>,
    pub regulators: Signal<BTreeSet<Rc<dyn Regulator>>>,
    serial: u64,
    column_index: Cell<Option<usize>>,
}

impl Point {
    const WEIGHT_COMPONENT: usize = 3;
    
    pub fn new(
        id: String,
        label: String,
        color: ElementColor,
        representation: DVector<f64>,
    ) -> Self {
        Self {
            id,
            label,
            color,
            representation: create_signal(representation),
            ghost: create_signal(false),
            regulators: create_signal(BTreeSet::new()),
            serial: Self::next_serial(),
            column_index: None.into(),
        }
    }
}

impl Element for Point {
    fn default_id() -> String {
        "point".to_string()
    }
    
    fn default(id: String, id_num: u64) -> Self {
        Self::new(
            id,
            format!("Point {id_num}"),
            [0.75_f32, 0.75_f32, 0.75_f32],
            point(0.0, 0.0, 0.0),
        )
    }
    
    fn id(&self) -> &String {
        &self.id
    }
    
    fn label(&self) -> &String {
        &self.label
    }
    
    fn representation(&self) -> Signal<DVector<f64>> {
        self.representation
    }
    
    fn ghost(&self) -> Signal<bool> {
        self.ghost
    }
    
    fn regulators(&self) -> Signal<BTreeSet<Rc<dyn Regulator>>> {
        self.regulators
    }
    
    fn project_to_normalized(&self, rep: &mut DVector<f64>) {
        project_point_to_normalized(rep);
    }
    
    fn column_index(&self) -> Option<usize> {
        self.column_index.get()
    }
    
    fn set_column_index(&self, index: usize) {
        self.column_index.set(Some(index));
    }
}

impl Serial for Point {
    fn serial(&self) -> u64 {
        self.serial
    }
}

impl ProblemPoser for Point {
    fn pose(&self, problem: &mut ConstraintProblem) {
        let index = self.column_index().expect(
            format!("Point \"{}\" should be indexed before writing problem data", self.id).as_str()
        );
        problem.gram.push_sym(index, index, 0.0);
        problem.frozen.push(Self::WEIGHT_COMPONENT, index, 0.5);
        problem.guess.set_column(index, &self.representation.get_clone_untracked());
    }
}

pub trait Regulator: Serial + ProblemPoser + OutlineItem {
    fn subjects(&self) -> Vec<Rc<dyn Element>>;
    fn measurement(&self) -> ReadSignal<f64>;
    fn set_point(&self) -> Signal<SpecifiedValue>;
}

impl Hash for dyn Regulator {
    fn hash<H: Hasher>(&self, state: &mut H) {
        <dyn Serial>::hash(self, state)
    }
}

impl PartialEq for dyn Regulator {
    fn eq(&self, other: &Self) -> bool {
        <dyn Serial>::eq(self, other)
    }
}

impl Eq for dyn Regulator {}

impl PartialOrd for dyn Regulator {
    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
        <dyn Serial>::partial_cmp(self, other)
    }
}

impl Ord for dyn Regulator {
    fn cmp(&self, other: &Self) -> Ordering {
        <dyn Serial>::cmp(self, other)
    }
}

pub struct InversiveDistanceRegulator {
    pub subjects: [Rc<dyn Element>; 2],
    pub measurement: ReadSignal<f64>,
    pub set_point: Signal<SpecifiedValue>,
    serial: u64,
}

impl InversiveDistanceRegulator {
    pub fn new(subjects: [Rc<dyn Element>; 2]) -> Self {
        let representations = subjects.each_ref().map(|subj| subj.representation());
        let measurement = create_memo(move || {
            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());
        let serial = Self::next_serial();
        
        Self { subjects, measurement, set_point, serial }
    }
}

impl Regulator for InversiveDistanceRegulator {
    fn subjects(&self) -> Vec<Rc<dyn Element>> {
        self.subjects.clone().into()
    }
    
    fn measurement(&self) -> ReadSignal<f64> {
        self.measurement
    }
    
    fn set_point(&self) -> Signal<SpecifiedValue> {
        self.set_point
    }
}

impl Serial for InversiveDistanceRegulator {
    fn serial(&self) -> u64 {
        self.serial
    }
}

impl ProblemPoser for InversiveDistanceRegulator {
    fn pose(&self, problem: &mut ConstraintProblem) {
        self.set_point.with_untracked(|set_pt| {
            if let Some(val) = set_pt.value {
                let [row, col] = self.subjects.each_ref().map(
                    |subj| subj.column_index().expect(
                        "Subjects should be indexed before inversive distance regulator writes problem data"
                    )
                );
                problem.gram.push_sym(row, col, val);
            }
        });
    }
}

pub struct HalfCurvatureRegulator {
    pub subject: Rc<dyn Element>,
    pub measurement: ReadSignal<f64>,
    pub set_point: Signal<SpecifiedValue>,
    serial: u64,
}

impl HalfCurvatureRegulator {
    pub fn new(subject: Rc<dyn Element>) -> Self {
        let measurement = subject.representation().map(
            |rep| rep[Sphere::CURVATURE_COMPONENT]
        );
        
        let set_point = create_signal(SpecifiedValue::from_empty_spec());
        let serial = Self::next_serial();
        
        Self { subject, measurement, set_point, serial }
    }
}

impl Regulator for HalfCurvatureRegulator {
    fn subjects(&self) -> Vec<Rc<dyn Element>> {
        vec![self.subject.clone()]
    }
    
    fn measurement(&self) -> ReadSignal<f64> {
        self.measurement
    }
    
    fn set_point(&self) -> Signal<SpecifiedValue> {
        self.set_point
    }
}

impl Serial for HalfCurvatureRegulator {
    fn serial(&self) -> u64 {
        self.serial
    }
}

impl ProblemPoser for HalfCurvatureRegulator {
    fn pose(&self, problem: &mut ConstraintProblem) {
        self.set_point.with_untracked(|set_pt| {
            if let Some(val) = set_pt.value {
                let col = self.subject.column_index().expect(
                    "Subject should be indexed before half-curvature regulator writes problem data"
                );
                problem.frozen.push(Sphere::CURVATURE_COMPONENT, col, val);
            }
        });
    }
}

// the velocity is expressed in uniform coordinates
pub struct ElementMotion<'a> {
    pub element: Rc<dyn Element>,
    pub velocity: DVectorView<'a, f64>,
}

type AssemblyMotion<'a> = Vec<ElementMotion<'a>>;

// a complete, view-independent description of an assembly
#[derive(Clone)]
pub struct Assembly {
    // elements and regulators
    pub elements: Signal<BTreeSet<Rc<dyn Element>>>,
    pub regulators: Signal<BTreeSet<Rc<dyn Regulator>>>,
    
    // 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<BTreeMap<String, Rc<dyn Element>>>,
    
    // realization control
    pub realization_trigger: Signal<()>,
    
    // realization diagnostics
    pub realization_status: Signal<Result<(), String>>,
    pub descent_history: Signal<DescentHistory>,
}

impl Assembly {
    pub fn new() -> Assembly {
        // create an assembly
        let assembly = Assembly {
            elements: create_signal(BTreeSet::new()),
            regulators: create_signal(BTreeSet::new()),
            tangent: create_signal(ConfigSubspace::zero(0)),
            elements_by_id: create_signal(BTreeMap::default()),
            realization_trigger: create_signal(()),
            realization_status: create_signal(Ok(())),
            descent_history: create_signal(DescentHistory::new()),
        };
        
        // 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 || {
            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
    }
    
    // --- inserting elements and regulators ---
    
    // insert an element 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: impl Element + 'static) {
        // insert the element
        let id = elt.id().clone();
        let elt_rc = Rc::new(elt);
        self.elements.update(|elts| elts.insert(elt_rc.clone()));
        self.elements_by_id.update(|elts_by_id| elts_by_id.insert(id, elt_rc.clone()));
        
        // create and insert the element's default regulators
        for reg in elt_rc.default_regulators() {
            self.insert_regulator(reg);
        }
    }
    
    pub fn try_insert_element(&self, elt: impl Element + 'static) -> bool {
        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);
        }
        can_insert
    }
    
    pub fn insert_element_default<T: Element + 'static>(&self) {
        // find the next unused identifier in the default sequence
        let default_id = T::default_id();
        let mut id_num = 1;
        let mut id = format!("{default_id}{id_num}");
        while self.elements_by_id.with_untracked(
            |elts_by_id| elts_by_id.contains_key(&id)
        ) {
            id_num += 1;
            id = format!("{default_id}{id_num}");
        }
        
        // create and insert the default example of `T`
        let _ = self.insert_element_unchecked(T::default(id, id_num));
    }
    
    pub fn insert_regulator(&self, regulator: Rc<dyn Regulator>) {
        // add the regulator to the assembly's regulator list
        self.regulators.update(
            |regs| regs.insert(regulator.clone())
        );
        
        // add the regulator to each subject's regulator list
        let subject_regulators: Vec<_> = regulator.subjects().into_iter().map(
            |subj| subj.regulators()
        ).collect();
        for regulators in subject_regulators {
            regulators.update(|regs| regs.insert(regulator.clone()));
        }
        
        /* DEBUG */
        // print an updated list of regulators
        console_log!("Regulators:");
        self.regulators.with_untracked(|regs| {
            for reg in regs.into_iter() {
                console_log!(
                    "  {:?}: {}",
                    reg.subjects(),
                    reg.set_point().with_untracked(
                        |set_pt| {
                            let spec = &set_pt.spec;
                            if spec.is_empty() {
                                "__".to_string()
                            } else {
                                spec.clone()
                            }
                        }
                    )
                );
            }
        });
    }
    
    // --- realization ---
    
    pub fn realize(&self) {
        // index the elements
        self.elements.update_silent(|elts| {
            for (index, elt) in elts.iter().enumerate() {
                elt.set_column_index(index);
            }
        });
        
        // 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);
            }
            self.regulators.with_untracked(|regs| {
                for reg in regs {
                    reg.pose(&mut problem);
                }
            });
            problem
        });
        
        /* DEBUG */
        // log the Gram matrix
        console_log!("Gram matrix:\n{}", problem.gram);
        
        /* DEBUG */
        // log the initial configuration matrix
        console_log!("Old configuration:{:>8.3}", problem.guess);
        
        // look for a configuration with the given Gram matrix
        let Realization { result, history } = realize_gram(
            &problem, 1.0e-12, 0.5, 0.9, 1.1, 200, 110
        );
        
        /* DEBUG */
        // report the outcome of the search in the browser console
        if let Err(ref message) = result {
            console_log!("❌️ {message}");
        } else {
            console_log!("✅️ Target accuracy achieved!");
        }
        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);
        
        match result {
            Ok(ConfigNeighborhood { config, nbhd: tangent }) => {
                /* DEBUG */
                // report the tangent dimension
                console_log!("Tangent dimension: {}", tangent.dim());
                
                // report the realization status
                self.realization_status.set(Ok(()));
                
                // 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()))
                    );
                }
                
                // save the tangent space
                self.tangent.set_silent(tangent);
            },
            Err(message) => {
                // report the realization status. the `Err(message)` we're
                // setting the status to has a different type than the
                // `Err(message)` we received from the match: we're changing the
                // `Ok` type from `Realization` to `()`
                self.realization_status.set(Err(message))
            },
        }
    }
    
    // --- 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 = {
            let mut next_column_index = realized_dim;
            for elt_motion in motion.iter() {
                let moving_elt = &elt_motion.element;
                if moving_elt.column_index().is_none() {
                    moving_elt.set_column_index(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 = elt_motion.element.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);
                let unif_to_std = elt_motion.element.representation().with_untracked(
                    |rep| local_unif_to_std(rep.as_view())
                );
                target_column += unif_to_std * elt_motion.velocity;
            }
        }
        
        // step the assembly along the deformation. this changes the elements'
        // normalizations, so we restore those afterward
        for elt in self.elements.get_clone_untracked() {
            elt.representation().update_silent(|rep| {
                match elt.column_index() {
                    Some(column_index) => {
                        // step the element along the deformation and then
                        // restore its normalization
                        *rep += motion_proj.column(column_index);
                        elt.project_to_normalized(rep);
                    },
                    None => {
                        console_log!("No velocity to unpack for fresh element \"{}\"", elt.id())
                    },
                };
            });
        }
        
        // 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.realization_trigger.set(());
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    
    use crate::engine;
    
    #[test]
    #[should_panic(expected = "Sphere \"sphere\" should be indexed before writing problem data")]
    fn unindexed_element_test() {
        let _ = create_root(|| {
            let elt = Sphere::default("sphere".to_string(), 0);
            elt.pose(&mut ConstraintProblem::new(1));
        });
    }
    
    #[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 subjects = [0, 1].map(
                |k| Rc::new(Sphere::default(format!("sphere{k}"), k)) as Rc<dyn Element>
            );
            subjects[0].set_column_index(0);
            InversiveDistanceRegulator {
                subjects: subjects,
                measurement: create_memo(|| 0.0),
                set_point: create_signal(SpecifiedValue::try_from("0.0".to_string()).unwrap()),
                serial: InversiveDistanceRegulator::next_serial()
            }.pose(&mut ConstraintProblem::new(2));
        });
    }
    
    #[test]
    fn curvature_drift_test() {
        const INITIAL_RADIUS: f64 = 0.25;
        let _ = create_root(|| {
            // set up an assembly containing a single sphere centered at the
            // origin
            let assembly = Assembly::new();
            let sphere_id = "sphere0";
            let _ = assembly.try_insert_element(
                // we create the sphere by hand for two reasons: to choose the
                // curvature (which can affect drift rate) and to make the test
                // independent of `Sphere::default`
                Sphere::new(
                    String::from(sphere_id),
                    String::from("Sphere 0"),
                    [0.75_f32, 0.75_f32, 0.75_f32],
                    engine::sphere(0.0, 0.0, 0.0, INITIAL_RADIUS),
                )
            );
            
            // nudge the sphere repeatedly along the `z` axis
            const STEP_SIZE: f64 = 0.0025;
            const STEP_CNT: usize = 400;
            let sphere = assembly.elements_by_id.with(|elts_by_id| elts_by_id[sphere_id].clone());
            let velocity = DVector::from_column_slice(&[0.0, 0.0, STEP_SIZE, 0.0]);
            for _ in 0..STEP_CNT {
                assembly.deform(
                    vec![
                        ElementMotion {
                            element: sphere.clone(),
                            velocity: velocity.as_view(),
                        }
                    ]
                );
            }
            
            // check how much the sphere's curvature has drifted
            const INITIAL_HALF_CURV: f64 = 0.5 / INITIAL_RADIUS;
            const DRIFT_TOL: f64 = 0.015;
            let final_half_curv = sphere.representation().with_untracked(
                |rep| rep[Sphere::CURVATURE_COMPONENT]
            );
            assert!((final_half_curv / INITIAL_HALF_CURV - 1.0).abs() < DRIFT_TOL);
        });
    }
}