feat: Find tangent space of solution variety, use for perturbations
### Tangent space
#### Implementation
The structure `engine::ConfigSubspace` represents a subspace of the configuration vector space $\operatorname{Hom}(\mathbb{R}^n, \mathbb{R}^5)$. It holds a basis for the subspace which is orthonormal with respect to the Euclidean inner product. The method `ConfigSubspace::symmetric_kernel` takes an endomorphism of the configuration vector space, which must be symmetric with respect to the Euclidean inner product, and returns its approximate kernel in the form of a `ConfigSubspace`.
At the end of `engine::realize_gram`, we use the computed Hessian to find the tangent space of the solution variety, and we return it alongside the realization. Since altering the constraints can change the tangent space without changing the solution, we compute the tangent space even when the guess passed to the realization routine is already a solution.
After `Assembly::realize` calls `engine::realize_gram`, it saves the returned tangent space in the assembly's `tangent` signal. The basis vectors are stored in configuration matrix format, ordered according to the elements' column indices. To help maintain consistency between the storage layout of the tangent space and the elements' column indices, we switch the column index data type from `usize` to `Option<usize>` and enforce the following invariants:
1. If an element has a column index, its tangent motions can be found in that column of the tangent space basis matrices.
2. If an element is affected by a constraint, it has a column index.
The comments in `assembly.rs` state the invariants and describe how they're enforced.
#### Automated testing
The test `engine::tests::tangent_test` builds a simple assembly with a known tangent space, runs the realization routine, and checks the returned tangent space against a hand-computed basis.
#### Limitations
The method `ConfigSubspace::symmetric_kernel` approximates the kernel by taking all the eigenspaces whose eigenvalues are smaller than a hard-coded threshold size. We may need a more flexible system eventually.
### Deformation
#### Implementation
The main purpose of this implementation is to confirm that deformation works as we'd hoped. The code is messy, and the deformation routine has at least one numerical quirk.
For simplicity, the keyboard commands that manipulate the assembly are handled by the display, just like the keyboard commands that control the camera. Deformation happens at the beginning of the animation loop.
The function `Assembly::deform` works like this:
1. Take a list of element motions
2. Project them onto the tangent space of the solution variety
3. Sum them to get a deformation $v$ of the whole assembly
4. Step the assembly along the "mass shell" geodesic tangent to $v$
* This step stays on the solution variety to first order
5. Call `realize` to bring the assembly back onto the solution variety
#### Manual testing
To manipulate the assembly:
1. Select a sphere
2. Make sure the display has focus
3. Hold the following keys:
* **A**/**D** for $x$ translation
* **W**/**S** for $y$ translation
* **shift**+**W**/**S** for $z$ translation
#### Limitations
Because the manipulation commands are handled by the display, you can only manipulate the assembly when the display has focus.
Since our test assemblies only include spheres, we assume in `Assembly::deform` that every element is a sphere.
When the tangent space is zero, `Assembly::deform` does nothing except print "The assembly is rigid" to the console.
During a deformation, the curvature and co-curvature components of a sphere's vector representation can exhibit weird discontinuous "swaps" that don't visibly affect how the sphere is drawn. *[I'll write more about this in an issue.]*
Co-authored-by: Aaron Fenyes <aaron.fenyes@fareycircles.ooo>
Reviewed-on: #29
Co-authored-by: Vectornaut <vectornaut@nobody@nowhere.net>
Co-committed-by: Vectornaut <vectornaut@nobody@nowhere.net>
On the incoming branch, you can select a sphere by clicking it in the display. Holding *shift* while clicking enables multiple selection. These controls match the ones already implemented in the outline view.
Since the selection routine is now used in multiple places, the incoming branch factors it out into the `AppState::select` method.
Co-authored-by: Aaron Fenyes <aaron.fenyes@fareycircles.ooo>
Reviewed-on: #25
Co-authored-by: Vectornaut <vectornaut@nobody@nowhere.net>
Co-committed-by: Vectornaut <vectornaut@nobody@nowhere.net>
Some of the Cargo tests on the main branch are designed to print output for human inspection, not to verify computations automatically. The incoming branch turns these tests into Cargo examples. It also makes two organizational changes in pursuit of this goal:
- It introduces a dyna3 library target, which the examples use as a dependency. In the future, this target could grow into an officially maintained dyna3 library.
- It puts the code for realizing the Irisawa hexlet into a new conditionally compiled `engine::irisawa` module. This code is shared by a test and an example. Compilation is controlled by the `dev` feature, which is turned on by default in development mode.
I've verified that printed output of the examples hasn't changed between the head (848f7d6) and base (e917272) of the incoming branch.
Co-authored-by: Aaron Fenyes <aaron.fenyes@fareycircles.ooo>
Co-authored-by: Glen Whitney <glen@studioinfinity.org>
Reviewed-on: #24
Co-authored-by: Vectornaut <vectornaut@nobody@nowhere.net>
Co-committed-by: Vectornaut <vectornaut@nobody@nowhere.net>
Give each `Element` a serial number, which identifies it uniquely. The serial number is assigned by the `Element::new` constructor.
Because disallows potentially unsafe global state (at least without explicit `unsafe` blocks), the next serial number is stored in a thread-safe static atomic variable (`assembly::NEXT_ELEMENT_SERIAL`), as suggested in [this StackOverflow answer](https://stackoverflow.com/a/32936288). Since the overhead for keeping track of memory ordering should be minimal, we're using the strongest available ordering: [sequentially consistent](https://marabos.nl/atomics/memory-ordering.html#seqcst).
Resolves#20.
Co-authored-by: Aaron Fenyes <aaron.fenyes@fareycircles.ooo>
Reviewed-on: #22
Co-authored-by: Vectornaut <vectornaut@nobody@nowhere.net>
Co-committed-by: Vectornaut <vectornaut@nobody@nowhere.net>
Clean up the source code and interface of the outline view. In addition, [fix a bug](commit/6e42681b719d7ec97c4225ca321225979bf87b56) that could cause `Assembly::realize` to react to itself under certain circumstances. Those circumstances arose, making the bug noticeable, while this branch was being written.
#### Source code
- Modularize the `Outline` component into smaller components.
- Switch from static iteration to dynamic Sycamore lists. This reduces the amount of re-rendering that happens when an element or constraint changes. It also allows constraint details to stay open or closed during constraint updates, rather than resetting to closed.
- Make `Element::index` private, as discussed [here](pulls/15#issuecomment-1816).
#### Interface
- Make constraints editable, updating the assembly realization on input. Flag constraints where the Lorentz product value doesn't parse.
- Round element vector coordinates to prevent the displayed strings from overlapping.
Note that issue #20 was created by this PR, but it will be addressed shortly.
Co-authored-by: Aaron Fenyes <aaron.fenyes@fareycircles.ooo>
Reviewed-on: #19
Co-authored-by: Vectornaut <vectornaut@nobody@nowhere.net>
Co-committed-by: Vectornaut <vectornaut@nobody@nowhere.net>
Port the engine prototype to Rust, integrate it into the application prototype, and use it to enforce the constraints.
### Features
To see the engine in action:
1. Add a constraint by shift-clicking to select two spheres in the outline view and then hitting the 🔗 button
2. Click a summary arrow to see the outline item for the new constraint
2. Set the constraint's Lorentz product by entering a value in the text field at the right end of the outline item
* *The display should update as soon as you press* Enter *or focus away from the text field*
The checkbox at the left end of a constraint outline item controls whether the constraint is active. Activating a constraint triggers a solution update. (Deactivating a constraint doesn't, since the remaining active constraints are still satisfied.)
### Precision
The Julia prototype of the engine uses a generic scalar type, so you can pass in any type the linear algebra functions are implemented for. The examples use the [adjustable-precision](https://docs.julialang.org/en/v1/base/numbers/#Base.MPFR.setprecision) `BigFloat` type.
In the Rust port of the engine, the scalar type is currently fixed at `f64`. Switching to generic scalars shouldn't be too hard, but I haven't looked into [which other types](https://www.nalgebra.org/docs/user_guide/generic_programming) the linear algebra functions are implemented for.
### Testing
To confirm quantitatively that the Rust port of the engine is working, you can go to the `app-proto` folder and:
* Run some automated tests by calling `cargo test`.
* Inspect the optimization process in a few examples calling the `run-examples` script. The first example that prints is the same as the Irisawa hexlet example from the engine prototype. If you go into `engine-proto/gram-test`, launch Julia, and then
```
include("irisawa-hexlet.jl")
for (step, scaled_loss) in enumerate(history_alt.scaled_loss)
println(rpad(step-1, 4), " | ", scaled_loss)
end
```
you should see that it prints basically the same loss history until the last few steps, when the lower default precision of the Rust engine really starts to show.
### A small engine revision
The Rust port of the engine improves on the Julia prototype in one part of the constraint-solving routine: projecting the Hessian onto the subspace where the frozen entries stay constant. The Julia prototype does this by removing the rows and columns of the Hessian that correspond to the frozen entries, finding the Newton step from the resulting "compressed" Hessian, and then adding zero entries to the Newton step in the appropriate places. The Rust port instead replaces each frozen row and column with its corresponding standard unit vector, avoiding the finicky compressing and decompressing steps.
To confirm that this version of the constraint-solving routine works the same as the original, I implemented it in Julia as `realize_gram_alt_proj`. The solutions we get from this routine match the ones we get from the original `realize_gram` to very high precision, and in the simplest examples (`sphere-in-tetrahedron.jl` and `tetrahedron-radius-ratio.jl`), the descent paths also match to very high precision. In a more complicated example (`irisawa-hexlet.jl`), the descent paths diverge about a quarter of the way into the search, even though they end up in the same place.
Co-authored-by: Aaron Fenyes <aaron.fenyes@fareycircles.ooo>
Reviewed-on: #15
Co-authored-by: Vectornaut <vectornaut@nobody@nowhere.net>
Co-committed-by: Vectornaut <vectornaut@nobody@nowhere.net>
Creates a prototype user interface for dyna3 in the `app-proto` folder. The interface is dynamically constructed using [Sycamore](https://sycamore.dev).
The prototype includes:
* An application state model (the `AppState` type)
* A constraint problem model (the `Assembly` type), used in the application state
* Two views
* A 3D rendering of the assembly (the `Display` component)
* A list of elements and constraints (the `Outline` component)
The following features confirm that the views can reflect and send input to the model:
* You can select elements by clicking and shift-clicking them in the outline. The selected elements are highlighted in the display.
* You can add elements using a button above the outline. The new elements appear in the display.
Co-authored-by: Aaron Fenyes <aaron.fenyes@fareycircles.ooo>
Reviewed-on: #14
Co-authored-by: Vectornaut <vectornaut@nobody@nowhere.net>
Co-committed-by: Vectornaut <vectornaut@nobody@nowhere.net>
This PR adds code for a Julia-language prototype of a configuration solver, in the `engine-proto` folder. It uses Julia version 1.10.0.
### Approaches
Development of this PR tried two broad approaches to the constraint geometry problem. Each one suggested various solution techniques. The Gram matrix approach, with the low-rank factorization technique, seems the most promising.
- **Algebraic** *(In the `alg-test` subfolder).* Write the constraints as polynomials in the inversive coordinates of the elements, and use computational algebraic geometry techniques to solve the resulting system. We tried the following techniques.
- **Gröbner bases** *(`Engine.Algebraic.jl`).* Symbolic. Find a Gröbner basis for the ideal generated by the constraint equations. Information about the solution variety, like its codimension, is then relatively easy to extract.
- **Homotopy continuation** *(`Engine.Numerical.jl`).* Numerical. Cut the solution set along a random hyperplane to get a generic zero-dimensional slice, and then use a fancy homotopy technique to approximate the points in that slice.
A few notes about our experiences can be found on the [engine prototype](wiki/Engine-prototype) wiki page.
- **Gram matrix** *(in the `gram-test` subfolder).* A construction is described completely, up to conformal transformations, by the Gram matrix of the vectors representing its elements. Express the constraints as fixed entries of the Gram matrix, and use numerical linear algebra techniques to find a list of vectors whose Gram matrix fits the bill. We tried the following techniques.
- **LDL decomposition** *(`gram-test.sage`, `gram-test.jl`, `overlap-test.jl`).* Find a cluster of up to five elements whose Gram matrix is completely filled in by the constraints. Use LDL decomposition to find a list of vectors with that Gram matrix. This technique can be made algebraic, as seen in `overlap-test.jl`.
- **Low-rank factorization** *(source files listed in findings section).* Write down a quadratic loss function that says how far a set of vectors is from meeting the Gram matrix constraints. Use a smooth optimization technique like Newton's method or gradient descent to find a zero of the loss function. In addition to the polished prototype described in the results section, we have an early prototype using an off-the-shelf factorization package (`low-rank-test.jl`) and an visualization of the loss function landscape near global minima (`basin-shapes.jl`).
The [Gram matrix parameterization](wiki/Gram-matrix-parameterization) wiki page contains detailed notes on this approach.
### Findings
With the algebraic approach, we hit a performance wall pretty quickly as our constructions grew. It was often hard to find real solutions of the polynomial system, since the techniques we use work most naturally in the complex world.
With the Gram matrix approach, on the other hand, we could solve interesting problems in acceptably short times using the low-rank factorization technique. We put the optimization routine in its own module (`Engine.jl`) and used it to solve five example problems:
- `overlapping-pyramids.jl`
- `circles-in-triangle.jl`
- `sphere-in-tetrahedron.jl`
- `tetrahedron-radius-ratio.jl`
- `irisawa-hexlet.jl`
We plan to use low-rank factorization of the Gram matrix in our first app prototype.
### Visualizations
We used the visualizer in the `ganja-test` folder to visually check our low-rank factorization results. The visualizer runs [Ganja.js](https://enkimute.github.io/ganja.js/) in an Electron app, made with [Blink](https://github.com/JuliaGizmos/Blink.jl). Although Ganja.js makes beautiful pictures under most circumstances, we found two obstacles to using it in production.
- It seems to have precision problems with low-curvature spheres.
- We couldn't figure out how to customize its clipping and transparency settings, and the default settings often obscure construction details.
Co-authored-by: Aaron Fenyes <aaron.fenyes@fareycircles.ooo>
Co-authored-by: Glen Whitney <glen@studioinfinity.org>
Reviewed-on: #13
Co-authored-by: Vectornaut <vectornaut@nobody@nowhere.net>
Co-committed-by: Vectornaut <vectornaut@nobody@nowhere.net>
This commit adds a utility to parse package-lock.json and write the proper
contents of externals.js to standard output. In addition, if the utility
(src/helpers/pkglock_to_externals.litcoffee) is invoked with a --doc option,
it instead emits a Markdown bulleted list of all of the external dependencies.
