Lazy evaluation of a continuous field of signals by Austreelis · Pull Request #600 · Leafwing-Studios/Emergence

Austreelis · 2023-03-24T04:06:22Z

Fixes #511

These patches implement an algorithm solving the diffusion equation "on-demand". It is still only a draft, and needs optimization to perform better than the fixed-time step diffusion passes.

Closed form of the solution to diffusion

As pointed out by @RobWalt here, solving the simplified diffusion equation in our case boils down to computing a convolution of a so-called fundamental solution and initial conditions. The fundamental solution is the solution with no border condition and a Dirac delta as initial conditions.

The approach I took may be different than what was outlined in the comment linked above, so I'll develop a bit further. I'll ignore border conditions (obstruction by structures) for now.

It would make sense to model emission of signals as Dirac delta functions (at the instant they are emitted), which then get diffused. This means that solving for 1 emission is trivial, this is just the fundamental solution:

Which we can adjust by keeping track of when and where signals are emitted:

u(x, t) = ϕ(x-x0, t-T0)

For several emissions (potentially emitted at different times), we could solve for one emission while the second hasn't been emitted, then do the complete convolution, and so forth for the following ones. However, convolutions are distributive, so the solution is the sum of the contributions of each emissions, no matter when they were emitted.

So generally for n emissions:

u(x, t) = ϕ(x-x0, t-T0) + ϕ(x-x1, t-T1) + ... + ϕ(x-xn, t-Tn)

This allows us to model the field exactly with 2 parameters per signal emissions.

In practice, 3. Because emissions also have a strength, representing an amount of signal emitted.

Expressing decay

The fundamental solution of the simplified diffusion equation we use represents a "density" of the diffused material. We could see it as pheromone concentration. The field being a density means we can integrate over its surface to get a total "amount" of pheromones. This "volume" is constant for the fundamental solution previously described, which means no signal strength is lost, it only gets displaced.

We could interpret decay as a gradual loss of signal strength, which would need us to use a fundamental solution that keeps the same shape as without decay, but with a decreasing volume. This is actually quite easy, we just need to add a dividing factor to the fundamental solution:

ϕ(x, t) = 1/(1 + D t) exp(-x^2/(4 k t)) / sqrt(4 pi k t)
          ^^^^^^^^^^^

The algorithm

We use 3 types:

SignalEmission, representing a single emission pulse: its location, the time when it was emitted, and its strength.
DiffusionEquation, representing a single equation, parametrized by several SignalEmissions, and the current time.
Signals, mapping signal types to their equation.

Every time the update_signals system runs, we increment the current time of each DiffusionEquation. We then prune SignalEmissions which are considered negligible: any emission whose strength is smaller than Th (1 + D t) (where D is the decay rate, t the time since emission, and Th some threshold) is removed from the equation. This is equivalent to pruning emissions whose volume under the curve is less than Th.

To emit a signal, we simply add a SignalEmission to the DiffusionEquation, with emission time set to the current time.

When a signal is accessed, we get the associated DiffusionEquation and iterate over its emissions. We compute the fundamental solution for each emission and accumulate them before returning the result. There is one subtlety: The distance is computed in tiles (using hexx' Hex::distance_to function), meaning that contour lines of equal signal strength are hexagons instead of circles.

Possible improvements

Some of the features the previous signals implementation had, and some that may be needed/nice later, are not handled yet by this algorithm:

Structures should block signal propagation.
Diffusivity may vary from place to place.

The amount of signal emissions stored seems to be the bottleneck, making this slower than the previous signal diffusion. But I'm pretty sure we could prune way more aggressively emissions, and that past some time their shape is so similar we can approximate a bunch with just one.

Approximating "old" emissions with only one.
Prune emissions more aggressively
Prune emissions based on the amount stored rather than on the time they were emitted.

There's also some other things to investigate:

When incrementing time, we iterate over each signal type, is this okay ?
Would caching field value per tile be useful ?
How to speed up the "debug" signal overlay ?

Todos

In no particular order:

I'm not really pleased with the fact I had to disable the `display_tile_overlay` system when using the signal overlay but this can be cleanup later. Color scheme has been updated to look better with the new terrain and shadows, and the overlay material has been made transparent a bit.

This only takes into account sources for now. This means that it only reaches part of the goal set in Leafwing-Studios#511: 1. It only tracks signal emitters. 2. It does determine an equation for each emitter. The equation simply is the solution of the diffusion equation for that emitter: `1/(1+D*t) * 1/sqrt(4*pi*k*t) * exp(d^2 / (4*k*t))` where: - `D` is a decay rate constant, and `1/(1+D*t)` is a decay term. - `k` is a diffusivity constant. - `d` is the distance between the emitter and a tile. - `t` is the time elapsed since the signal was emitted. 3. Contributions of each emitter are then integrated to compute the actual field value. This works because of the linearity of this subset of the problem, but may not hold when tackling e.g. barriers or signal modifiers (I haven't yet investigated it). 4. It lazyly looks up the field value when needed, but I suspect this is pretty costly when displaying the signal overlay. I don't know if and what we'd want to do about it yet.

Our application is in 2D. The 1D fundamental solution was used previously.

Austreelis · 2023-03-24T04:10:39Z

I did not have the time to push this draft as soon as I wanted, but here it is.

There's 2 issues right now that are on top of my priority list:

Emissions pruning logic is broken
The debug signal overlay is flickering. I realized the issue after rebasing on main, but it's getting too late for me to investigate.

Then next up probably is a bit of clean up and design talk, before I attempt tackling obstructions.

JasmineLowen · 2023-03-24T07:22:12Z

Very well done! ⭐ I love that you implemented the debug overlay aswell! ✨ Would you be interested in a review?

Austreelis · 2023-03-24T12:38:40Z

Would you be interested in a review?

Yes please ! There's a lot to improve, so I'd gladly get feedback early to go in the right direction 😄

alice-i-cecile · 2023-03-24T13:35:54Z

Super cool! I'm going to start by splitting out the signal-visualization work: this is valuable regardless of the approach used and I think we can handle it more robustly in its own PR :)

Co-authored-by: Austreelis <dev@austreelis.net>

* Extract initial signal viz code from #600 Co-authored-by: Austreelis <dev@austreelis.net> * Move tile selection code into infoviz module * Store an Option for visualized signal type * Use a single resource to store overlay info * Refactor material fetching into a method * Log error, don't panic * Reparameterize color ramp scaling for clarity * Fix math * Tweak max signal strength for better visualization * Use N_COLORS consistently * Display color ramp in legend * Organize color palette constants * Don't screw up the base color * Get tile interactions and signal viz to play nice together * Create color ramp via HSLA color interpolation * Store the legend in the right orientation * Tweak overlay colors for now * Fix math for normalizing signal strength * Attempt to do UI layout * Display currently visualized signal * Cycle signal types randomly * Clippy --------- Co-authored-by: Austreelis <dev@austreelis.net>

alice-i-cecile · 2023-03-24T23:42:15Z

+use super::{SignalStrength, DECAY_RATE, DECAY_THRESHOLD, DIFFUSIVITY};
+
+#[derive(Clone, Copy, Debug)]
+pub struct SignalEmission {


Doc strings for the type as a whole please :)

alice-i-cecile · 2023-03-24T23:42:28Z

+///
+/// where:
+/// - ϕ is the density of the diffusing material. In emergence's case, ϕ is the signal strength.
+/// - D is the diffusivity of the environment. For emergence, this is a constant for now, making


Suggested change

/// - D is the diffusivity of the environment. For emergence, this is a constant for now, making

/// - D is the diffusivity of the environment. For Emergence, this is a constant for now, making

alice-i-cecile · 2023-03-24T23:43:03Z

+///
+/// The fundamental solution solves a problem with no boundary conditions and
+/// an initial density given by the [Dirac delta function](https://en.wikipedia.org/wiki/Dirac_delta_function).
+/// The general solution, for any initial condition is a convolution with the fundamental solution:


Suggested change

/// The general solution, for any initial condition is a convolution with the fundamental solution:

/// The general solution for any initial condition is a convolution with the fundamental solution:

alice-i-cecile · 2023-03-24T23:44:20Z

+/// - y is a surface element in the same space as x.
+///
+/// Practically, we don't want to deal with costly convolutions, so we'll try to avoid them.
+/// Because signals are emitted in pulses, and convolution is distributive, we can separate each


Do we want signals to be implemented in pulses? I think ideally these would be better modeled as plateaus, to avoid situations where wave fronts propagate away from the actual source.

alice-i-cecile · 2023-03-24T23:44:49Z

+/// s(x, t) = ∫ Φ(x−y, t-T0) e0(y) dy + ∫ Φ(x−y, t-T1) e1(y) dy + ... + ∫ Φ(x−y, t-Tn) en(y) dy
+///
+/// where:
+/// - en is the initial condition of the nth emission.


I think that in comments, this is clearer as e_n

alice-i-cecile · 2023-03-24T23:45:27Z

+/// s(x, t) = Φ(x-y0, t-T0) + Φ(x-y1, t-T1) + ... + Φ(x-yn, t-Tn)
+///
+/// where yn is the location of the nth signal emission.
+// NOTE: One equation per signal type implies they don't interact together.


This should be a doc comment on emissions

alice-i-cecile · 2023-03-24T23:48:07Z

+        });
+    }
+
+    /// Sets the time at which the equation will be evaluated to the given `current_time`. This


Nit: the explanatory text belongs in a new paragraph; this plays nicer with IDE tooltip hovers.

alice-i-cecile · 2023-03-24T23:49:58Z

+                // negligible being more for fast decaying signals (we'll trim it earlier), and
+                // less for slowly decaying ones (we'll trim it later)
+                Some(SignalEmission { time, strength, .. })
+                    if DECAY_THRESHOLD * (1.0 + DECAY_RATE * (self.time - time)) > strength => {}


I'm having trouble following the math here when skimming; I think intermediate variables will make this clearer.

alice-i-cecile · 2023-03-24T23:58:02Z

Very well-made, and fascinating math.

Main conclusions:

We can't replace the existing solution until a) this handles obstructions (pathfinding around obstacles feels far too important) and b) this is actually a performance improvement. I think a feature flag is the way to go.
Ultimately I think that modelling emitters as pulsing Dirac deltas is the wrong approach conceptually, and will substantially increase costs. Instead, we should be treating this as piece-wise linear, and using (careful) change detection to determine the knots.
Very clever to extinguish sufficiently old events! That's something that we're going to want to build on. I think we'll eventually want to do a two pass system here. Periodically, we'll want to clean up old events that won't matter even if you're on top of them. And then, when computing the signal at a specific location, you'll only want to fetch nearby events. Not needed for an MVP though.
Fetching a "patch" of 7 signals of the same type is a super common operation, and is probably worth abstracting out here in order to unlock dedicated performance optimizations.
I've extracted the signal visualization code from here, so you can now delete it :)

Austreelis · 2023-03-25T07:09:34Z

About 4: I'm guessing this is done by ants to follow gradients of signal, right ? If so caching results for the time of one signal pass can probably be an easy win, given how ants tends to be close together.

alice-i-cecile · 2023-03-25T19:34:36Z

About 4: I'm guessing this is done by ants to follow gradients of signal, right ? If so caching results for the time of one signal pass can probably be an easy win, given how ants tends to be close together.

Correct :)

Austreelis · 2023-04-02T22:40:19Z

I'm unlikely to have time to work on this PR before after the bevy jam, but I'm still looking at it. If anyone has suggestions/changes to add, I can probably take a quick look. And I'll continue after the Jam :)

Austreelis added 6 commits March 24, 2023 02:44

Add a "debug" visualization of signal strengths

9823d9f

Move signal.rs -> signals/mod.rs

2c89ef6

Fix the formula of the fundamental solution of the diffusion equation

e9cf281

Our application is in 2D. The 1D fundamental solution was used previously.

Document maths of the diffusion equation solver

07bc38f

alice-i-cecile added performance Code go brrr architectural Changes to the organization and patterns of the code labels Mar 24, 2023

alice-i-cecile added this to the Foundations milestone Mar 24, 2023

alice-i-cecile added a commit that referenced this pull request Mar 24, 2023

Extract initial signal viz code from #600

37fa8c2

Co-authored-by: Austreelis <dev@austreelis.net>

alice-i-cecile mentioned this pull request Mar 24, 2023

Visualize signal strengths by location #601

Merged

5 tasks

alice-i-cecile added a commit that referenced this pull request Mar 24, 2023

Extract initial signal viz code from #600

e820c78

Co-authored-by: Austreelis <dev@austreelis.net>

alice-i-cecile reviewed Mar 24, 2023

View reviewed changes

alice-i-cecile added the controversial Needs careful analysis and testing before proceeding label Mar 25, 2023

alice-i-cecile removed this from the Foundations milestone May 10, 2023

	/// - D is the diffusivity of the environment. For emergence, this is a constant for now, making
	/// - D is the diffusivity of the environment. For Emergence, this is a constant for now, making

	/// The general solution, for any initial condition is a convolution with the fundamental solution:
	/// The general solution for any initial condition is a convolution with the fundamental solution:

Uh oh!

Conversation

Austreelis commented Mar 24, 2023 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Closed form of the solution to diffusion

Expressing decay

The algorithm

Possible improvements

Todos

Uh oh!

Austreelis commented Mar 24, 2023

Uh oh!

JasmineLowen commented Mar 24, 2023

Uh oh!

Austreelis commented Mar 24, 2023

Uh oh!

alice-i-cecile commented Mar 24, 2023

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

alice-i-cecile commented Mar 24, 2023 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Austreelis commented Mar 25, 2023

Uh oh!

alice-i-cecile commented Mar 25, 2023

Uh oh!

Austreelis commented Apr 2, 2023

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

Austreelis commented Mar 24, 2023 •

edited

Loading

alice-i-cecile commented Mar 24, 2023 •

edited

Loading