diff --git a/packages/cultmath/README.md b/packages/cultmath/README.md index 552df549..10881fdb 100644 --- a/packages/cultmath/README.md +++ b/packages/cultmath/README.md @@ -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. diff --git a/packages/cultmath/docs/design.md b/packages/cultmath/docs/design.md index ecd66d5a..d6abafd8 100644 --- a/packages/cultmath/docs/design.md +++ b/packages/cultmath/docs/design.md @@ -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 diff --git a/packages/cultmath/src/CultMath/math.Lattice.cs b/packages/cultmath/src/CultMath/math.Lattice.cs new file mode 100644 index 00000000..42ae963b --- /dev/null +++ b/packages/cultmath/src/CultMath/math.Lattice.cs @@ -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. + + /// The greatest boost accepts: weights span 0.5 to 5, a ratio of 10:1. + public const float LatticeMaxBoost = 4.0f; + + /// The greatest stretch accepts. + public const float LatticeMaxStretch = 2.0f; + + /// The greatest smoothing accepts, in cell units. + public const float LatticeMaxSmoothing = 1.0f; + + /// + /// The outermost ring the search can reach anywhere in the accepted range (see the proof above + /// ). + /// + 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; + + /// + /// The site of a lattice cell: float3(site.xy, base weight). The site is cell + 0.5 + 0.5 sin(drift + /// time + 2 pi u) per axis, so it stays within one cell of the cell's corner and drifts at rate + /// . The base weight is 0.5 + 0.5 u, in [0.5, 1), from its own hash lane. + /// + 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)); + } + + /// + /// The weight of a cell: base (1 + boost exp(-(border / falloff)^2)), where is + /// the distance from the cell's centre to the nearest body-diagram border. + /// + 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))); + } + + /// + /// The weighted distance from to , in cell units: + /// weight |(stretch a, c)| with a and c the offset's components along and across the unit + /// . Isotropic at stretch 1. + /// + 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)); + } + + /// + /// One step of the smooth minimum fold: folded with distance , for a + /// nearest distance . A distance of at least f1 + 2 smoothing leaves + /// unchanged, and the result never exceeds f1 + smoothing. Start the fold at + /// f1 + smoothing and visit every cell once, in an order that does not depend on the search box. + /// + 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. + + /// + /// Weighted cellular noise over a jittered lattice: float4(F1, F2, smooth, id). F1 and F2 are the two + /// smallest weighted distances from in cell units times weight, smooth is the smooth minimum + /// (, in the same units), and id is a hash of the winning cell in [0, 1). + /// + /// + /// A cell's weight is of the exact power border distance () + /// of the cell's centre in the body span (float4(position.xyz, sigma), z = 0), so it is + /// continuous as the bodies move. Distances are stretched by along + /// (isotropic at 1). The search is exact in the accepted range: boost 0 to + /// , stretch 1 to , smoothing up to . + /// + /// 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 refuses it. + public static float4 weighted_lattice_cells( + float2 p, float time, ReadOnlySpan 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 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 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); + } +} diff --git a/packages/cultmath/src/CultMath/math.Lattice.cs.meta b/packages/cultmath/src/CultMath/math.Lattice.cs.meta new file mode 100644 index 00000000..a4f24ecd --- /dev/null +++ b/packages/cultmath/src/CultMath/math.Lattice.cs.meta @@ -0,0 +1,2 @@ +fileFormatVersion: 2 +guid: 0a182e056c8c4258881c347ead6d571d diff --git a/packages/cultmath/src/CultMath/math.Sites.cs b/packages/cultmath/src/CultMath/math.Sites.cs index feab825c..a5a4937b 100644 --- a/packages/cultmath/src/CultMath/math.Sites.cs +++ b/packages/cultmath/src/CultMath/math.Sites.cs @@ -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 sites, bool positiveWeights) + private static void require_cell_span(float3 p, ReadOnlySpan 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 sites) { require_cell_span(p, sites, positiveWeights: false); + return power_cells_unchecked(p, sites); + } + + private static CultSiteCells power_cells_unchecked(float3 p, ReadOnlySpan 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 /// /// As . public static float4 power_border(float3 p, ReadOnlySpan 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 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++) diff --git a/packages/cultmath/tests/CultMath.Tests/LatticeCellsTests.cs b/packages/cultmath/tests/CultMath.Tests/LatticeCellsTests.cs new file mode 100644 index 00000000..741c3d58 --- /dev/null +++ b/packages/cultmath/tests/CultMath.Tests/LatticeCellsTests.cs @@ -0,0 +1,460 @@ +using Xunit; +using static CultMath.math; + +namespace CultMath.Tests; + +/// +/// 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. +/// +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(); + 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(), 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(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()); + 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(); + 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(); + 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(); + 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); + } + } +}