GameCult/CultLib
Watch
9
0
Fork
You've already forked CultLib
 0
Code Issues Pull requests 1 Projects Releases Packages Wiki Activity Actions

Merge hands/cultmath-lattice-weighted: lattice-weighted-cellular (4a)

Browse source 
Spec cultmath-tapes:cut_spec:cut-lattice-weighted-cellular.r1; reports h1, h2; Soul s2 holds.
No conflicts.

Co-Authored-By: Claude Sonnet 5.5 <noreply@anthropic.com>
This commit is contained in:
Meta Vens 2026-10-17 11:14:56 +02:00
commit b496131a6d
6 changed files with 723 additions and 5 deletions
Download patch file Download diff file Expand all files Collapse all files

4
packages/cultmath/README.md
View file

@ -34,6 +34,10 @@ shader semantics.
`soft_gabriel_gate` gates a pair of sites by the empty-disc margin. Each is a loop over the shader steps
`cultmath_power_score` / `cultmath_multiplicative_score`, `cultmath_site_cells_empty`, `cultmath_site_cells_insert`,
`cultmath_site_cells_border`, `cultmath_power_border_step` (the exact power cell border, `power_border`) and `cultmath_gabriel_margin`; a shader loops over its own buffer.
- Weighted lattice cells (C# only until the shader cut): `weighted_lattice_cells` returns `float4(F1, F2, smooth, id)` over a jittered,
drifting lattice whose cell weights come from the exact power border of the body span at each cell centre, with a stretched
metric and an exact ring search over the stated range (`LatticeMaxBoost`, `LatticeMaxStretch`, `LatticeMaxSmoothing`). It composes the
public steps `lattice_site`, `lattice_weight`, `lattice_distance` and `lattice_smooth_step`.
- `shaders/CultMath.hlsl`, a canonical HLSL mirror include for shader-side
parity. HLSL already owns `float2`, `float3`, and `float4`; the include
exposes `cultmath_*` functions for shared semantics HLSL does not provide.

28
packages/cultmath/docs/design.md
View file

@ -532,6 +532,34 @@ is `smoothstep(-softness, softness, least margin over every other site)` on posi
exists; it is C# only. `distance_to_segment(point, start, end)` has out-free 2D and 3D forms beside the `out`
form, all over one public `segment_closest` per dimension, each with its 2D and 3D shader overload.
## Weighted lattice cells
`weighted_lattice_cells(p, time, bodies, cellSize, drift, boost, falloff, smoothing, stretchNormal, stretch)` in
`math.Lattice.cs` returns `float4(F1, F2, smooth, id)` over a jittered, drifting lattice, in cell units times weight.
Four public steps compose it, and a shader or a tape composes the same four:
- `lattice_site(cell, time, drift)`: the site `cell + 0.5 + 0.5 sin(drift time + 2 pi u)` per axis (so it stays within
one cell of the corner) and the base weight `0.5 + 0.5 u`, in [0.5, 1), from its own hash lane.
- `lattice_weight(base, border, boost, falloff)`: `base (1 + boost exp(-(border / falloff)^2))`. `border` is
`power_border` at the cell centre over the body span (exact, so the weight is continuous as bodies move; the first
order `site_cells_border` jumps where the second-nearest body changes and is not used).
- `lattice_distance(q, site, weight, normal, stretch)`: `weight |(stretch a, c)|`, `a` and `c` the offset along and
across the unit normal; isotropic at stretch 1.
- `lattice_smooth_step(m, d, f1, smoothing)`: the smooth minimum fold, started at `F1 + smoothing` and clamped to it,
so a cell with `d >= F1 + 2 smoothing` changes nothing and the result does not depend on the box visited.
The search visits Chebyshev rings around the query's cell and stops before ring `k` once `0.5 (k - 1) >= max(F2,
F1 + 2 smoothing)`: every cell on ring `k` has `d >= 0.5 (k - 1)` (site within one cell of its corner, weight at
least 0.5, stretch at least 1). A second pass folds the visited box in absolute row-major order, so ties and the
smooth minimum do not depend on how far the search went. The proof is the comment above `weighted_lattice_cells`.
The accepted range is boost 0 to `LatticeMaxBoost` (4, a weight ratio of 10:1), stretch 1 to `LatticeMaxStretch` (2)
and smoothing up to `LatticeMaxSmoothing` (1); outside it the function throws `ArgumentException`. `LatticeMaxRing`
(37) is the proven greatest ring at the range edge. Measured over 24,000 points on orbiting bodies (six bodies,
sigma ratio 6:1, the ruled look and the range edge), the greatest ring visited has median 3, 99th percentile 7 and
maximum 8 (49, 225 and 289 cells); `RingsStayWellInsideTheCap` asserts these. The function is C# only; the shader
twins are the next cut.
## Intervals
`math.Interval.cs` and `shaders/CultMath.Interval.hlsl` (included from

212
packages/cultmath/src/CultMath/math.Lattice.cs Normal file
View file

@ -0,0 +1,212 @@
namespace CultMath;
public static partial class math
{
// Weighted lattice cells: a jittered, drifting cellular field whose cells are weighted by how close the cell's
// centre lies to a border of the body diagram (the power diagram of the zone's body span). The public steps are
// the four per-cell rules a tape or a shader composes; weighted_lattice_cells is the loop over them.
/// <summary>The greatest boost <see cref="weighted_lattice_cells"/> accepts: weights span 0.5 to 5, a ratio of 10:1.</summary>
public const float LatticeMaxBoost = 4.0f;
/// <summary>The greatest stretch <see cref="weighted_lattice_cells"/> accepts.</summary>
public const float LatticeMaxStretch = 2.0f;
/// <summary>The greatest smoothing <see cref="weighted_lattice_cells"/> accepts, in cell units.</summary>
public const float LatticeMaxSmoothing = 1.0f;
/// <summary>
/// The outermost ring the search can reach anywhere in the accepted range (see the proof above
/// <see cref="weighted_lattice_cells"/>).
/// </summary>
public const int LatticeMaxRing = 37;
// A query beyond this many cells from the origin has no integer cell id in float precision.
private const float LatticeMaxCoordinate = 8388608.0f;
/// <summary>
/// The site of a lattice cell: <c>float3(site.xy, base weight)</c>. The site is <c>cell + 0.5 + 0.5 sin(drift
/// time + 2 pi u)</c> per axis, so it stays within one cell of the cell's corner and drifts at rate
/// <paramref name="drift"/>. The base weight is <c>0.5 + 0.5 u</c>, in [0.5, 1), from its own hash lane.
/// </summary>
public static float3 lattice_site(float2 cell, float time, float drift)
{
var h = pcg3d(new int3((int)cell.x, (int)cell.y, 0));
var phase = new float2(cellular_unit(h.x), cellular_unit(h.y));
var site = cell + 0.5f + 0.5f * new float2(sin(drift * time + TAU * phase.x), sin(drift * time + TAU * phase.y));
return new float3(site, 0.5f + 0.5f * cellular_unit(h.z));
}
/// <summary>
/// The weight of a cell: <c>base (1 + boost exp(-(border / falloff)^2))</c>, where <paramref name="border"/> is
/// the distance from the cell's centre to the nearest body-diagram border.
/// </summary>
public static float lattice_weight(float @base, float border, float boost, float falloff)
{
var t = border / falloff;
return @base * (1.0f + boost * exp(-(t * t)));
}
/// <summary>
/// The weighted distance from <paramref name="q"/> to <paramref name="site"/>, in cell units:
/// <c>weight |(stretch a, c)|</c> with <c>a</c> and <c>c</c> the offset's components along and across the unit
/// <paramref name="normal"/>. Isotropic at stretch 1.
/// </summary>
public static float lattice_distance(float2 q, float2 site, float weight, float2 normal, float stretch)
{
var v = site - q;
var a = dot(v, normal);
var c = v.x * normal.y - v.y * normal.x;
return weight * length(new float2(stretch * a, c));
}
/// <summary>
/// One step of the smooth minimum fold: <paramref name="m"/> folded with distance <paramref name="d"/>, for a
/// nearest distance <paramref name="f1"/>. A distance of at least <c>f1 + 2 smoothing</c> leaves
/// <paramref name="m"/> unchanged, and the result never exceeds <c>f1 + smoothing</c>. Start the fold at
/// <c>f1 + smoothing</c> and visit every cell once, in an order that does not depend on the search box.
/// </summary>
public static float lattice_smooth_step(float m, float d, float f1, float smoothing)
{
if (d >= f1 + 2.0f * smoothing)
return m;
var h = smoothstep(0.0f, 1.0f, 0.5f + 0.5f * (m - d) / smoothing);
var next = m + (d - m) * h - h * (1.0f - h) * smoothing / (1.0f + 3.0f * smoothing);
return min(next, f1 + smoothing);
}
// Proof that the ring search finds every cell that can change an output (ruling lattice-ring-exit-sound).
//
// Let c0 = floor(q) and let a cell lie on Chebyshev ring k around c0 (an offset of k on some axis). Its site lies
// within [cell, cell + 1] per axis and q within [c0, c0 + 1), so on that axis |site - q| >= k - 1. The length of
// the offset is at least that, the stretch is at least 1 (it only lengthens), and the weight is at least
// 0.5 (base >= 0.5, boost >= 0). So every cell on ring k or beyond has d >= 0.5 (k - 1), and the same holds for
// every farther ring.
//
// A cell with d >= F2 changes neither F1 nor F2, and a cell with d >= F1 + 2 smoothing leaves the smooth minimum
// unchanged (lattice_smooth_step returns m before folding it). So the search stops before ring k once
// 0.5 (k - 1) >= max(F2, F1 + 2 smoothing), with F1 and F2 taken over the rings visited. Those values only fall
// as rings are added, so the running cutoff is at least the final one: stopping on it never stops early. Every
// unvisited cell then changes no output, for every weight, stretch and drift.
//
// The loop terminates inside the accepted range. Cell c0 gives a distance below sqrt(2) in cell units. The
// neighbour on the nearer side of the axis whose query fraction sits closest to a cell boundary (fraction m <= 1/2)
// lies within (1 + m, 1) per axis, below sqrt(3.25). Weights are below 1 + LatticeMaxBoost and stretch is at most
// LatticeMaxStretch, so F2 < 10 sqrt(3.25) = 18.03 and F1 + 2 smoothing < 10 sqrt(2) + 2 = 16.2. The exit holds at
// 0.5 (k - 1) >= 18.03, that is at k = 38, so rings 0 to 37 (LatticeMaxRing) are the most visited.
//
// The result does not depend on the box visited. The second pass folds F1, F2 and the winning id in absolute
// row-major order (strict <, so an exact tie keeps the earlier cell), and the smooth minimum starts at
// F1 + smoothing and is clamped to it, so a cell with d >= F1 + 2 smoothing, visited or not, changes nothing.
/// <summary>
/// Weighted cellular noise over a jittered lattice: <c>float4(F1, F2, smooth, id)</c>. F1 and F2 are the two
/// smallest weighted distances from <paramref name="p"/> in cell units times weight, smooth is the smooth minimum
/// (<paramref name="smoothing"/>, in the same units), and id is a hash of the winning cell in [0, 1).
/// </summary>
/// <remarks>
/// A cell's weight is <see cref="lattice_weight"/> of the exact power border distance (<see cref="power_border"/>)
/// of the cell's centre in the body span <paramref name="bodies"/> (<c>float4(position.xyz, sigma)</c>, z = 0), so it is
/// continuous as the bodies move. Distances are stretched by <paramref name="stretch"/> along
/// <paramref name="stretchNormal"/> (isotropic at 1). The search is exact in the accepted range: boost 0 to
/// <see cref="LatticeMaxBoost"/>, stretch 1 to <see cref="LatticeMaxStretch"/>, smoothing up to <see cref="LatticeMaxSmoothing"/>.
/// </remarks>
/// <exception cref="ArgumentException">An argument is not finite, a size is not above zero, a value is outside its range, the normal is zero while stretch exceeds 1, or the body span is refused as <see cref="power_cells"/> refuses it.</exception>
public static float4 weighted_lattice_cells(
float2 p, float time, ReadOnlySpan<float4> bodies, float cellSize, float drift, float boost, float falloff,
float smoothing, float2 stretchNormal, float stretch)
{
// Cell centres stay finite: the farthest one the search can reach is a ring past the query's cell.
var reachSpan = cellSize * (LatticeMaxRing + 2);
if (!float.IsFinite(cellSize) || !(cellSize > 0.0f) || !float.IsFinite(reachSpan))
throw new ArgumentException("cellSize must be finite and greater than zero, and small enough for its cell centres to be finite.", nameof(cellSize));
if (!float.IsFinite(p.x) || !float.IsFinite(p.y) || !(abs(p.x / cellSize) < LatticeMaxCoordinate) || !(abs(p.y / cellSize) < LatticeMaxCoordinate)
|| !float.IsFinite(abs(p.x) + reachSpan) || !float.IsFinite(abs(p.y) + reachSpan))
throw new ArgumentException("p must be finite and within 2^23 cells of the origin.", nameof(p));
if (!float.IsFinite(time))
throw new ArgumentException("time must be finite.", nameof(time));
if (!float.IsFinite(drift))
throw new ArgumentException("drift must be finite.", nameof(drift));
if (!float.IsFinite(drift * time))
throw new ArgumentException("drift times time must be finite.", nameof(drift));
if (!float.IsFinite(boost) || boost < 0.0f || boost > LatticeMaxBoost)
throw new ArgumentException("boost must be within [0, LatticeMaxBoost].", nameof(boost));
if (!float.IsFinite(falloff) || !(falloff > 0.0f))
throw new ArgumentException("falloff must be finite and greater than zero.", nameof(falloff));
if (!float.IsFinite(smoothing) || !(smoothing > 0.0f) || smoothing > LatticeMaxSmoothing)
throw new ArgumentException("smoothing must be within (0, LatticeMaxSmoothing].", nameof(smoothing));
if (!float.IsFinite(stretch) || stretch < 1.0f || stretch > LatticeMaxStretch)
throw new ArgumentException("stretch must be within [1, LatticeMaxStretch].", nameof(stretch));
if (!float.IsFinite(stretchNormal.x) || !float.IsFinite(stretchNormal.y))
throw new ArgumentException("stretchNormal must be finite.", nameof(stretchNormal));
// Scaled by its largest component first, so a tiny or huge normal normalises to unit length.
var scale = max(abs(stretchNormal.x), abs(stretchNormal.y));
if (scale > 0.0f)
{
stretchNormal /= scale;
stretchNormal /= length(stretchNormal);
}
else if (stretch > 1.0f)
throw new ArgumentException("stretchNormal must not be zero while stretch exceeds 1.", nameof(stretchNormal));
else
stretchNormal = new float2(1.0f, 0.0f);
require_cell_span(new float3(p, 0.0f), bodies, positiveWeights: false, nameof(bodies));
const int Width = 2 * LatticeMaxRing + 1;
Span<float> cache = stackalloc float[Width * Width];
var q = p / cellSize;
var origin = floor(q);
var cx0 = (int)origin.x;
var cz0 = (int)origin.y;
var f1 = 1.0e30f;
var f2 = 1.0e30f;
var reach = 0;
for (var k = 0; ; k++)
{
if (k > 1 && 0.5f * (k - 1) >= max(f2, f1 + 2.0f * smoothing))
break;
if (k > LatticeMaxRing)
throw new InvalidOperationException("The lattice ring search left its proven range.");
for (var dz = -k; dz <= k; dz++)
{
var step = dz == -k || dz == k || k == 0 ? 1 : 2 * k;
for (var dx = -k; dx <= k; dx += step)
{
var d = lattice_cell_distance(q, new float2(cx0 + dx, cz0 + dz), time, drift, bodies, cellSize, boost, falloff, stretchNormal, stretch);
cache[(dz + LatticeMaxRing) * Width + dx + LatticeMaxRing] = d;
if (d < f1) { f2 = f1; f1 = d; }
else if (d < f2) f2 = d;
}
}
reach = k;
}
var m = f1 + smoothing;
var best = 1.0e30f;
var second = 1.0e30f;
var winDx = 0;
var winDz = 0;
for (var dz = -reach; dz <= reach; dz++)
for (var dx = -reach; dx <= reach; dx++)
{
var d = cache[(dz + LatticeMaxRing) * Width + dx + LatticeMaxRing];
if (d < best) { second = best; best = d; winDx = dx; winDz = dz; }
else if (d < second) second = d;
m = lattice_smooth_step(m, d, f1, smoothing);
}
var id = pcg4d(new int4(cx0 + winDx, cz0 + winDz, 0, 0));
return new float4(best, second, m, cellular_unit(id.w));
}
private static float lattice_cell_distance(
float2 q, float2 cell, float time, float drift, ReadOnlySpan<float4> bodies, float cellSize, float boost, float falloff,
float2 normal, float stretch)
{
var site = lattice_site(cell, time, drift);
var centre = (cell + 0.5f) * cellSize;
var border = power_border_unchecked(new float3(centre, 0.0f), bodies).w;
return lattice_distance(q, site.xy, lattice_weight(site.z, border, boost, falloff), normal, stretch);
}
}

2
packages/cultmath/src/CultMath/math.Lattice.cs.meta Normal file
View file

@ -0,0 +1,2 @@
fileFormatVersion: 2
guid: 0a182e056c8c4258881c347ead6d571d

22
packages/cultmath/src/CultMath/math.Sites.cs
View file

@ -99,18 +99,18 @@ public static partial class math
// The one validation of a cell span: finite query, finite positions and weights, weights above zero when
// positiveWeights, at most 2^24 sites. Messages name the parameter, never a value.
private static void require_cell_span(float3 p, ReadOnlySpan<float4> sites, bool positiveWeights)
private static void require_cell_span(float3 p, ReadOnlySpan<float4> sites, bool positiveWeights, string spanName = "sites")
{
if (!all_finite(p))
throw new ArgumentException("The query point must be finite.", nameof(p));
if (sites.Length > MaxSiteSpan)
throw new ArgumentException("A site span holds at most 2^24 sites.", nameof(sites));
throw new ArgumentException("A site span holds at most 2^24 sites.", spanName);
for (var i = 0; i < sites.Length; i++)
{
if (!all_finite(sites[i].xyz))
throw new ArgumentException("Site positions must be finite.", nameof(sites));
throw new ArgumentException("Site positions must be finite.", spanName);
if (!float.IsFinite(sites[i].w) || (positiveWeights && !(sites[i].w > 0.0f)))
throw new ArgumentException(positiveWeights ? "Site weights must be finite and greater than zero." : "Site weights must be finite.", nameof(sites));
throw new ArgumentException(positiveWeights ? "Site weights must be finite and greater than zero." : "Site weights must be finite.", spanName);
}
}
@ -122,6 +122,11 @@ public static partial class math
public static CultSiteCells power_cells(float3 p, ReadOnlySpan<float4> sites)
{
require_cell_span(p, sites, positiveWeights: false);
return power_cells_unchecked(p, sites);
}
private static CultSiteCells power_cells_unchecked(float3 p, ReadOnlySpan<float4> sites)
{
var cells = site_cells_empty();
for (var i = 0; i < sites.Length; i++)
cells = site_cells_insert(cells, power_score(p, sites[i]), i);
@ -137,9 +142,16 @@ public static partial class math
/// </summary>
/// <exception cref="ArgumentException">As <see cref="power_cells"/>.</exception>
public static float4 power_border(float3 p, ReadOnlySpan<float4> sites)
{
require_cell_span(p, sites, positiveWeights: false);
return power_border_unchecked(p, sites);
}
// The border over an already validated span, which the lattice calls once per cell.
private static float4 power_border_unchecked(float3 p, ReadOnlySpan<float4> sites)
{
var border = new float4(0.0f, 0.0f, 0.0f, 1.0e30f);
var owner = (int)power_cells(p, sites).nearestSite;
var owner = (int)power_cells_unchecked(p, sites).nearestSite;
if (owner < 0)
return border;
for (var k = 0; k < sites.Length; k++)

460
packages/cultmath/tests/CultMath.Tests/LatticeCellsTests.cs Normal file
View file

@ -0,0 +1,460 @@
using Xunit;
using static CultMath.math;
namespace CultMath.Tests;
/// <summary>
/// The weighted lattice cell function and its public steps. The oracle and the double reference are written here
/// from the public steps and the spec formula; the ring search is never the reference for itself.
/// </summary>
public sealed class LatticeCellsTests
{
private const float Cell = 260.0f;
private const float Falloff = 120.0f;
private const float Drift = 0.7f;
private static readonly float2 Normal = new(0.6f, 0.8f);
private sealed record Look(float Boost, float Stretch, float Smoothing, float Falloff = LatticeCellsTests.Falloff);
private static readonly Look[] Looks =
{
new(1.6f, 1.0f, 0.35f), new(1.6f, 2.0f, 0.35f), new(LatticeMaxBoost, LatticeMaxStretch, LatticeMaxSmoothing), new(0.0f, 1.0f, 0.05f),
};
private static float4[] Bodies(float t)
{
var bodies = new float4[6];
for (var i = 0; i < bodies.Length; i++)
{
var angle = 0.001f * (i + 1) * t + i;
var radius = 300.0f + 220.0f * i;
bodies[i] = new float4(radius * MathF.Cos(angle), radius * MathF.Sin(angle), 0.0f, 40.0f + 40.0f * i);
}
return bodies;
}
private static float CellWeight(float2 cell, float time, float4[] bodies, Look look)
{
var site = lattice_site(cell, time, Drift);
var border = power_border(new float3((cell + 0.5f) * Cell, 0.0f), bodies).w;
return lattice_weight(site.z, border, look.Boost, look.Falloff);
}
// Folds F1, F2, id and smooth over the box of the given radius in absolute row-major order.
private static float4 BoxFold(float2 p, float time, float4[] bodies, Look look, int radius,
Dictionary<(int, int), float> weights, float2? unitNormal = null)
{
var q = p / Cell;
var c0 = floor(q);
var cells = new List<(float2 cell, float d)>();
for (var dz = -radius; dz <= radius; dz++)
for (var dx = -radius; dx <= radius; dx++)
{
var cell = c0 + new float2(dx, dz);
var key = ((int)cell.x, (int)cell.y);
if (!weights.TryGetValue(key, out var w))
weights[key] = w = CellWeight(cell, time, bodies, look);
var site = lattice_site(cell, time, Drift);
cells.Add((cell, lattice_distance(q, site.xy, w, unitNormal ?? Normal, look.Stretch)));
}
float f1 = 1.0e30f, f2 = 1.0e30f;
var win = default(float2);
foreach (var (cell, d) in cells)
{
if (d < f1) { f2 = f1; f1 = d; win = cell; }
else if (d < f2) f2 = d;
}
var m = f1 + look.Smoothing;
foreach (var (_, d) in cells)
m = lattice_smooth_step(m, d, f1, look.Smoothing);
var id = pcg4d(new int4((int)win.x, (int)win.y, 0, 0));
return new float4(f1, f2, m, cellular_unit(id.w));
}
private static float4 Call(float2 p, float time, float4[] bodies, Look look, float2? normal = null) =>
weighted_lattice_cells(p, time, bodies, Cell, Drift, look.Boost, look.Falloff, look.Smoothing, normal ?? Normal, look.Stretch);
private static float2 RandomPoint(System.Random r) =>
new((r.NextSingle() * 2.0f - 1.0f) * 1500.0f, (r.NextSingle() * 2.0f - 1.0f) * 1500.0f);
[Fact]
public void RingSearchEqualsTheFullBoxOracle()
{
var r = new System.Random(41);
foreach (var look in Looks)
for (var batch = 0; batch < 20; batch++)
{
var time = 13.0f * batch;
var bodies = Bodies(time);
var weights = new Dictionary<(int, int), float>();
for (var i = 0; i < 1000; i++)
{
var p = RandomPoint(r);
var expected = BoxFold(p, time, bodies, look, 12, weights);
var actual = Call(p, time, bodies, look);
Assert.Equal(expected, actual);
}
}
}
private delegate bool Stop(int k, float f1, float f2, Look look);
private static readonly (string name, Stop stop)[] Rules =
{
("correct", (k, f1, f2, l) => k > 1 && 0.5f * (k - 1) >= max(f2, f1 + 2.0f * l.Smoothing)),
("cutoff F1 + k", (k, f1, f2, l) => k > 1 && 0.5f * (k - 1) >= max(f2, f1 + l.Smoothing)),
("bound 0.5*k", (k, f1, f2, l) => k > 0 && 0.5f * k >= max(f2, f1 + 2.0f * l.Smoothing)),
("greatest weight", (k, f1, f2, l) => k > 1 && (1.0f + l.Boost) * (k - 1) >= max(f2, f1 + 2.0f * l.Smoothing)),
("F1 only", (k, f1, f2, l) => k > 1 && 0.5f * (k - 1) >= f1 + 2.0f * l.Smoothing),
};
// The radius of the box a search with the given exit rule visits (rings 0 to k - 1 when it stops before ring k).
private static int Radius(float2 p, float time, float4[] bodies, Look look, Stop stop,
Dictionary<(int, int), float> weights)
{
var q = p / Cell;
var c0 = floor(q);
float f1 = 1.0e30f, f2 = 1.0e30f;
for (var k = 0; k <= LatticeMaxRing; k++)
{
if (stop(k, f1, f2, look))
return k - 1;
for (var dz = -k; dz <= k; dz++)
for (var dx = -k; dx <= k; dx++)
{
if (Math.Max(Math.Abs(dx), Math.Abs(dz)) != k)
continue;
var cell = c0 + new float2(dx, dz);
var key = ((int)cell.x, (int)cell.y);
if (!weights.TryGetValue(key, out var w))
weights[key] = w = CellWeight(cell, time, bodies, look);
var d = lattice_distance(q, lattice_site(cell, time, Drift).xy, w, Normal, look.Stretch);
if (d < f1) { f2 = f1; f1 = d; }
else if (d < f2) f2 = d;
}
}
return LatticeMaxRing;
}
[Fact]
public void ConstructedPointsRealizeTheExitBound()
{
var failed = new Dictionary<string, int>();
var r = new System.Random(7);
foreach (var look in new[] { Looks[2], Looks[3] })
for (var batch = 0; batch < 20; batch++)
{
var time = 29.0f * batch;
var bodies = Bodies(time);
var weights = new Dictionary<(int, int), float>();
for (var i = 0; i < 1500; i++)
{
var p = RandomPoint(r);
var expected = BoxFold(p, time, bodies, look, 12, weights);
foreach (var (name, stop) in Rules)
{
var got = BoxFold(p, time, bodies, look, Math.Max(0, Radius(p, time, bodies, look, stop, weights)), weights);
if (!got.Equals(expected))
failed[name] = failed.GetValueOrDefault(name) + 1;
}
}
}
Assert.False(failed.ContainsKey("correct"));
foreach (var (name, _) in Rules.Skip(1))
Assert.True(failed.GetValueOrDefault(name) > 0, name + " is not caught by any point");
}
[Fact]
public void SmoothIgnoresCellsBeyondTheBand()
{
var r = new System.Random(3);
for (var i = 0; i < 20000; i++)
{
var k = 0.05f + 0.95f * r.NextSingle();
var f1 = 8.0f * r.NextSingle();
var m = f1 + k - (r.NextSingle() < 0.5f ? 0.0f : 2.0f * k * r.NextSingle());
var far = f1 + 2.0f * k + 5.0f * r.NextSingle();
Assert.Equal(m, lattice_smooth_step(m, far, f1, k));
var near = f1 + 2.0f * k * r.NextSingle();
Assert.True(lattice_smooth_step(m, near, f1, k) <= f1 + k);
}
Assert.NotEqual(1.0f, lattice_smooth_step(1.0f + 0.35f, 1.0f + 0.35f, 1.0f, 0.35f));
}
private static float2 Dir(System.Random r)
{
var a = r.NextSingle() * TAU;
return new float2(MathF.Cos(a), MathF.Sin(a));
}
[Fact]
public void ContinuousAcrossCellBoundariesAndBorders()
{
var look = Looks[1];
var bodies = Bodies(11.0f);
var r = new System.Random(5);
var worst = 0.0f;
for (var line = 0; line < 6; line++)
{
var a = RandomPoint(r);
var dir = Dir(r);
const int steps = 4000;
var length = 6.0f * Cell;
var prev = Call(a, 11.0f, bodies, look);
for (var s = 1; s <= steps; s++)
{
var cur = Call(a + dir * (length * s / steps), 11.0f, bodies, look);
var step = length / steps / Cell;
// A weighted distance changes at most weight * stretch per cell unit of travel.
Assert.True(abs(cur.z - prev.z) < 12.0f * step + 1e-3f);
Assert.True(abs(cur.x - prev.x) < 12.0f * step + 1e-5f * 12.0f);
worst = max(worst, abs((cur.x - cur.z) - (prev.x - prev.z)));
prev = cur;
}
}
Assert.True(worst < 12.0f * 6.0f / 4000 * 2.0f + 1e-3f);
// Straddling a lattice cell boundary by 1e-3 of a cell.
for (var i = 0; i < 300; i++)
{
var p = RandomPoint(r);
var edge = new float2(MathF.Round(p.x / Cell) * Cell, p.y);
var lo = Call(edge - new float2(1e-3f * Cell, 0.0f), 11.0f, bodies, look);
var hi = Call(edge + new float2(1e-3f * Cell, 0.0f), 11.0f, bodies, look);
Assert.True(abs(lo.x - hi.x) < 1e-3f * 24.0f);
Assert.True(abs(lo.y - hi.y) < 1e-3f * 24.0f);
Assert.True(abs(lo.z - hi.z) < 1e-3f * 24.0f);
Assert.True(abs((lo.x - lo.z) - (hi.x - hi.z)) < 1e-3f * 24.0f);
}
}
[Fact]
public void WeightsAreContinuousAsBodiesOrbit()
{
var look = Looks[2];
var cells = Enumerable.Range(0, 576).Select(i => new float2(i % 24 - 12, i / 24 - 12)).ToArray();
var prev = new float[cells.Length];
for (var s = 0; s <= 300 * 60; s++)
{
var bodies = Bodies(s / 60.0f);
for (var i = 0; i < cells.Length; i++)
{
var site = lattice_site(cells[i], 0.0f, Drift);
var border = power_border(new float3((cells[i] + 0.5f) * Cell, 0.0f), bodies).w;
var w = lattice_weight(site.z, border, look.Boost, Falloff) / site.z;
if (s > 0)
Assert.True(abs(w - prev[i]) <= 0.01f * look.Boost, $"cell {cells[i]} at {s}: {prev[i]} -> {w}");
prev[i] = w;
}
}
}
[Fact]
public void F1F2AndSmoothAreContinuousThroughWinnerSwaps()
{
var r = new System.Random(19);
var bodies = Bodies(11.0f);
var swaps = 0;
foreach (var look in new[] { Looks[1], Looks[2] })
for (var line = 0; line < 6; line++)
{
var a = RandomPoint(r);
var b = RandomPoint(r);
float2 At(float t) => a + (b - a) * t;
const int Steps = 4000;
for (var s = 1; s <= Steps; s++)
{
float lo = (s - 1) / (float)Steps, hi = s / (float)Steps;
if (Call(At(lo), 11.0f, bodies, look).w == Call(At(hi), 11.0f, bodies, look).w)
continue;
swaps++;
for (var i = 0; i < 40; i++)
{
var mid = 0.5f * (lo + hi);
if (Call(At(mid), 11.0f, bodies, look).w == Call(At(lo), 11.0f, bodies, look).w) lo = mid; else hi = mid;
}
var before = Call(At(lo), 11.0f, bodies, look);
var after = Call(At(hi), 11.0f, bodies, look);
var bound = 12.0f * length(At(hi) - At(lo)) / Cell + 1e-5f * 12.0f;
Assert.True(abs(before.x - after.x) < bound, $"F1 {before.x} -> {after.x}");
Assert.True(abs(before.y - after.y) < bound, $"F2 {before.y} -> {after.y}");
Assert.True(abs(before.z - after.z) < bound, $"smooth {before.z} -> {after.z}");
}
}
Assert.True(swaps >= 20, $"only {swaps} winner swaps crossed");
}
[Fact]
public void TheCapIsReachedAndMatchesTheFullBox()
{
// A cell whose every neighbour is boosted far and stretched needs ring 17, past the cap a smaller bound would set.
var look = new Look(LatticeMaxBoost, LatticeMaxStretch, LatticeMaxSmoothing, 1.0e6f);
var p = new float2(1446.7345f, -680.7637f);
const float Time = 97.69269f;
var bodies = Bodies(Time);
var weights = new Dictionary<(int, int), float>();
var ring = Radius(p, Time, bodies, look, Rules[0].stop, weights);
Assert.InRange(ring, 13, LatticeMaxRing);
Assert.Equal(BoxFold(p, Time, bodies, look, LatticeMaxRing, weights), Call(p, Time, bodies, look));
// The cap is the greatest ring the range's worst case can need: the exit holds once 0.5 (k - 1) reaches the
// larger of the F2 and F1 + 2 smoothing bounds, which are (1 + boost) stretch sqrt(3.25) and (1 + boost) stretch sqrt(2) + 2 smoothing.
var f2Bound = (1.0f + LatticeMaxBoost) * LatticeMaxStretch * MathF.Sqrt(3.25f);
var bandBound = (1.0f + LatticeMaxBoost) * LatticeMaxStretch * MathF.Sqrt(2.0f) + 2.0f * LatticeMaxSmoothing;
Assert.Equal((int)MathF.Ceiling(2.0f * max(f2Bound, bandBound)), LatticeMaxRing);
}
[Fact]
public void AnExactTieChangesOnlyTheIdToTheEarlierRowMajorCell()
{
// Cells sharing a base weight and, at these coordinates, a rounded site: a query on that site is at distance 0 from both.
foreach (var (earlier, time) in new[] { (new float2(6003486, 3996 + 6000000), 0.0f), (new float2(6007572, 2713 + 6000000), 4.07f) })
{
var later = earlier + new float2(0, 1);
var siteA = lattice_site(earlier, time, 1.0f);
Assert.Equal(siteA, lattice_site(later, time, 1.0f));
var got = weighted_lattice_cells(siteA.xy, time, Array.Empty<float4>(), 1.0f, 1.0f, 0.0f, 1.0f, 0.01f, Normal, 1.0f);
var idOf = (float2 c) => cellular_unit(pcg4d(new int4((int)c.x, (int)c.y, 0, 0)).w);
Assert.NotEqual(idOf(earlier), idOf(later));
Assert.Equal(new float4(0.0f, 0.0f, got.z, idOf(earlier)), got);
}
}
private static void Refuses(string parameter, Action call) =>
Assert.Equal(parameter, Assert.Throws<ArgumentException>(call).ParamName);
[Fact]
public void RefusesOutsideTheRange()
{
var b = Bodies(0.0f);
var p = new float2(10.0f, 20.0f);
var n = Normal;
float4 W(float2? pp = null, float time = 0, float cell = Cell, float drift = Drift, float boost = 1.6f, float falloff = Falloff,
float smoothing = 0.35f, float2? normal = null, float stretch = 1.5f, float4[]? bodies = null) =>
weighted_lattice_cells(pp ?? p, time, bodies ?? b, cell, drift, boost, falloff, smoothing, normal ?? n, stretch);
W(boost: 0.0f); W(boost: LatticeMaxBoost); W(stretch: 1.0f); W(stretch: LatticeMaxStretch);
W(smoothing: LatticeMaxSmoothing); W(normal: float2.zero, stretch: 1.0f); W(normal: new float2(1.0e-45f, 0.0f)); W(bodies: Array.Empty<float4>());
Refuses("p", () => W(pp: new float2(float.NaN, 0)));
Refuses("cellSize", () => W(cell: 0.0f));
Refuses("cellSize", () => W(cell: float.PositiveInfinity));
Refuses("cellSize", () => W(cell: 1.0e37f));
Refuses("p", () => W(pp: new float2(3.0e38f, 0)));
Refuses("p", () => W(pp: new float2(0, 1.0e10f)));
Refuses("p", () => W(pp: new float2(3.3e38f, 0), cell: 1.0e36f));
Refuses("time", () => W(time: float.PositiveInfinity));
Refuses("drift", () => W(drift: float.NaN));
Refuses("drift", () => W(time: 3.0e38f, drift: 2.0f));
W(time: 1.0e30f, drift: 1.0e-30f);
Refuses("boost", () => W(boost: -0.01f));
Refuses("boost", () => W(boost: LatticeMaxBoost + 0.01f));
Refuses("falloff", () => W(falloff: 0.0f));
Refuses("falloff", () => W(falloff: float.NaN));
Refuses("smoothing", () => W(smoothing: 0.0f));
Refuses("smoothing", () => W(smoothing: LatticeMaxSmoothing + 0.01f));
Refuses("stretch", () => W(stretch: 0.99f));
Refuses("stretch", () => W(stretch: LatticeMaxStretch + 0.01f));
Refuses("stretchNormal", () => W(normal: float2.zero));
Refuses("stretchNormal", () => W(normal: new float2(float.NaN, 1)));
Refuses("bodies", () => W(bodies: new[] { new float4(float.NaN, 0, 0, 1) }));
Refuses("bodies", () => W(bodies: new[] { new float4(0, 0, 0, float.PositiveInfinity) }));
}
[Fact]
public void RingsStayWellInsideTheCap()
{
var r = new System.Random(41);
var rings = new List<int>();
foreach (var look in Looks)
for (var batch = 0; batch < 20; batch++)
{
var time = 13.0f * batch;
var bodies = Bodies(time);
var weights = new Dictionary<(int, int), float>();
for (var i = 0; i < 300; i++)
rings.Add(Radius(RandomPoint(r), time, bodies, look, Rules[0].stop, weights));
}
Assert.True(rings.Max() <= LatticeMaxRing - 2, $"greatest ring {rings.Max()}");
rings.Sort();
Assert.Equal((3, 7, 8), (rings[rings.Count / 2], rings[rings.Count * 99 / 100], rings[^1]));
}
[Fact]
public void StretchNormalIsNormalised()
{
var bodies = Bodies(3.0f);
var r = new System.Random(8);
for (var i = 0; i < 50; i++)
{
var p = RandomPoint(r);
Assert.Equal(Call(p, 3.0f, bodies, Looks[1]), Call(p, 3.0f, bodies, Looks[1], Normal * 7.0f));
}
}
[Fact]
public void ANormalOfAnyMagnitudeActsAsItsDirection()
{
// A tiny normal once normalised to length 0.944 and a huge one to zero; both then missed a winning cell.
var none = Array.Empty<float4>();
var look = new Look(0.0f, 1.0f, 0.01f);
var invSqrt2 = new float2(0.70710677f, 0.70710677f);
var cases = new (float2 given, float2 unit)[]
{
(new float2(3.0e-23f, 4.0e-23f), Normal), (new float2(6.0e29f, 8.0e29f), Normal),
(new float2(1.0e30f, 1.0e30f), invSqrt2), (new float2(1.0e-44f, 1.0e-44f), invSqrt2), (new float2(3.0e38f, 3.0e38f), invSqrt2),
};
var weights = new Dictionary<(int, int), float>();
foreach (var (given, unit) in cases)
{
var miss = new float2(-1040.0f, -1481.7773f);
Assert.Equal(BoxFold(miss, 54.725235f, none, look, 12, weights, unit), Call(miss, 54.725235f, none, look, given));
var r = new System.Random(3);
for (var i = 0; i < 300; i++)
{
var p = RandomPoint(r);
Assert.Equal(BoxFold(p, 7.0f, none, Looks[1], 12, weights, unit), Call(p, 7.0f, none, Looks[1], given));
}
}
}
[Fact]
public void MatchesTheProbeFormula()
{
var look = Looks[1];
var r = new System.Random(2);
var bodies = Bodies(5.0f);
for (var i = 0; i < 2000; i++)
{
var p = RandomPoint(r);
var got = Call(p, 5.0f, bodies, look);
var nrm = (x: (double)Normal.x, y: (double)Normal.y);
var q = (x: (double)p.x / Cell, y: (double)p.y / Cell);
var ds = new List<double>();
for (var cz = (int)Math.Floor(q.y) - 8; cz <= (int)Math.Floor(q.y) + 8; cz++)
for (var cx = (int)Math.Floor(q.x) - 8; cx <= (int)Math.Floor(q.x) + 8; cx++)
{
var h = pcg3d(new int3(cx, cz, 0));
double u0 = cellular_unit(h.x), u1 = cellular_unit(h.y), u2 = cellular_unit(h.z);
var ph = 2 * Math.PI;
var sx = cx + 0.5 + 0.5 * Math.Sin(Drift * 5.0 + ph * u0);
var sy = cz + 0.5 + 0.5 * Math.Sin(Drift * 5.0 + ph * u1);
var border = power_border(new float3(((float)cx + 0.5f) * Cell, ((float)cz + 0.5f) * Cell, 0), bodies).w;
var weight = (0.5 + 0.5 * u2) * (1 + look.Boost * Math.Exp(-Math.Pow(border / (double)Falloff, 2)));
double vx = sx - q.x, vy = sy - q.y;
double a = vx * nrm.x + vy * nrm.y, c = vx * nrm.y - vy * nrm.x;
ds.Add(weight * Math.Sqrt(look.Stretch * look.Stretch * a * a + c * c));
}
var d = ds.OrderBy(x => x).ToArray();
double m = d[0] + look.Smoothing;
foreach (var dd in ds)
{
if (dd >= d[0] + 2 * look.Smoothing) continue;
var h = Math.Clamp(0.5 + 0.5 * (m - dd) / look.Smoothing, 0, 1);
h = h * h * (3 - 2 * h);
m = Math.Min(m + (dd - m) * h - h * (1 - h) * look.Smoothing / (1 + 3 * look.Smoothing), d[0] + look.Smoothing);
}
Assert.InRange(got.x - d[0], -1e-4, 1e-4);
Assert.InRange(got.y - d[1], -1e-4, 1e-4);
Assert.InRange(got.z - m, -1e-4, 1e-4);
}
}
}