diff --git a/packages/cultmath/tests/CultMath.Tests/CellularAndSminGradTests.cs b/packages/cultmath/tests/CultMath.Tests/CellularAndSminGradTests.cs index c401fbfa..465f6961 100644 --- a/packages/cultmath/tests/CultMath.Tests/CellularAndSminGradTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/CellularAndSminGradTests.cs @@ -6,11 +6,32 @@ namespace CultMath.Tests; public sealed class CellularAndSminGradTests { - // Shared central-difference tolerance for smin_grad and cellular (F5). Soul measured float noise - // at the committed eps values (0.001f here, 0.0005f in cellular) at about 6e-5; this leaves - // roughly 5x headroom for legitimate finite-difference truncation while still killing a - // 1%-scale gradient mutant (about 1e-2 error against these fields). - private const float GradientTolerance = 3.0e-4f; + // Central-difference tolerances for smin_grad and cellular (F5), sized separately: the two + // functions have different curvature, so the same eps does not carry the same truncation + // error, and an earlier shared 3e-4 tolerance was fragile against both (Soul, cut 2a-i: measured + // max error smin 2.74e-4, cellular nearest 1.25e-4, cellular edge 1.88e-4 on a small sample — + // and re-measured here at 5.3e-4, 2.0e-4, 2.6e-4 on a larger one — against that 3e-4 tolerance, + // sometimes under 1.1x headroom, not the "~6e-5 noise" an earlier version of this comment + // claimed). Each tolerance below is set from its own measured worst case with real headroom. + + // smin_grad's value is exactly quadratic in position within a smoothing band (math.cs comment on + // smin_grad: h is affine in b.w - a.w, which is itself affine in position, so value(q) has zero + // third derivative there): central differences have no truncation error from curvature, only + // float32 rounding noise, so a larger eps directly dilutes that noise with no accuracy cost. + // Measured max error at eps=0.01f over 10,000 samples: 4.37e-5 (up from eps=0.001f's 5.3e-4, + // confirming the noise is eps-limited, not curvature-limited); this tolerance keeps about 5.7x + // headroom over that. A uniform-scale mutant on the returned gradient (a 1.01x M9/M16-style + // mutant applied here) still dies down to a 1.0005x scale (0.05%); it survives at 1.0002x. + private const float SminGradientTolerance = 2.5e-4f; + + // cellular's F1/F2 are Euclidean distance fields: real curvature away from a feature point means + // central differences DO carry O(eps^2) truncation error, so (unlike smin_grad) a larger eps + // makes things worse, not better; a sweep at eps in {0.0002f, 0.0003f, 0.0005f, 0.0007f, 0.001f} + // over 8,000 samples found the committed eps=0.0005f already close to the sweet spot. Measured + // max error at that eps: nearest 2.04e-4, edge 2.57e-4; this tolerance keeps just over 5x + // headroom over the larger of the two. M9 (grad1 * 1.01) and M16 (grad2 * 1.01) both still die + // down to a 1.002x scale (0.2%); both survive at 1.001x. + private const float CellularGradientTolerance = 1.3e-3f; // ---- smin_grad ---- @@ -81,24 +102,31 @@ public sealed class CellularAndSminGradTests // Stay clear of the |a.w - b.w| = k seam: central differences straddling a kink would // disagree with either side's analytic gradient (the same reason cellular excludes a - // band around its own F1 = F2 set). + // band around its own F1 = F2 set). The band has an absolute floor, not just a + // fraction of k: eps grew to 0.01f below (SminGradientTolerance's comment), and a + // fraction-of-k-only band could be narrower than the diff can shift by eps*|cf-cg| + // (up to about 0.035 here) when k is small, letting p +/- eps step across the seam even + // though p itself was clear of it. var diff = F(p) - G(p); - if (Math.Abs(Math.Abs(diff) - k) < 0.05f * k) + if (Math.Abs(Math.Abs(diff) - k) < Math.Max(0.05f * k, 0.05f)) continue; tested++; var analytic = math.smin_grad(new float4(cf, F(p)), new float4(cg, G(p)), k); - const float eps = 0.001f; + + // eps=0.01f, not 0.001f: smin is exactly quadratic within a band (no truncation error + // from curvature; see SminGradientTolerance's comment), so the larger eps only dilutes + // float32 rounding noise, at no accuracy cost, while still comfortably inside the band + // exclusion above. + const float eps = 0.01f; var numeric = new float3( (Smin(p + new float3(eps, 0.0f, 0.0f)) - Smin(p - new float3(eps, 0.0f, 0.0f))) / (2.0f * eps), (Smin(p + new float3(0.0f, eps, 0.0f)) - Smin(p - new float3(0.0f, eps, 0.0f))) / (2.0f * eps), (Smin(p + new float3(0.0f, 0.0f, eps)) - Smin(p - new float3(0.0f, 0.0f, eps))) / (2.0f * eps)); - // Measured float noise at this eps is about 6e-5 (Soul, cut 2a-i); 3e-4 leaves roughly - // 5x headroom while still killing a 1%-scale gradient mutant (about 1e-2 error here). - AssertWithinTolerance(numeric.x, analytic.x, GradientTolerance); - AssertWithinTolerance(numeric.y, analytic.y, GradientTolerance); - AssertWithinTolerance(numeric.z, analytic.z, GradientTolerance); + AssertWithinTolerance(numeric.x, analytic.x, SminGradientTolerance); + AssertWithinTolerance(numeric.y, analytic.y, SminGradientTolerance); + AssertWithinTolerance(numeric.z, analytic.z, SminGradientTolerance); } } @@ -145,11 +173,12 @@ public sealed class CellularAndSminGradTests { // math.cellular's own search order visits neighbour offsets with dz outermost and dx // innermost over a 5x5x5 block (F6), and it also prunes cells whose position-only lower - // bound already exceeds the running F2. Offsets with any |component| = 2 are provably never - // the true nearest (their minimum possible distance, sqrt(3), sits right at the maximum - // possible distance of the query's own cell, per the exactness proof in math.cs, so it takes - // a measure-zero coincidence for one to win — confirmed empirically at 0 wins in 3,000,000 - // samples), which only leaves the 3x3x3 corners as ever-reachable "last" cases. (dx, dy, dz) + // bound already exceeds the running F2. Offsets with any |component| = 2 are not provably + // unreachable as F1: their minimum possible distance is 1 (e.g. (2, 0, 0), whose lower bound + // per the exactness proof in math.cs is max(0, 2-1) = 1), well inside the query's own cell's + // sqrt(3) ceiling, so a |component| = 2 offset winning F1 is rare rather than impossible + // (Soul, cut 2a-i: measured 3 F1 wins in 2,000,000 samples). It is still rare enough that the + // 3x3x3 corners remain the practically-reachable "last" cases for this test. (dx, dy, dz) // = (1, 1, 1) is the last of those the fixed iteration order visits, and if pruning or search // order were broken, this is the offset most likely to be skipped. Rather than hand-picking a // hash value, search random points for one whose true nearest feature (by independent brute @@ -202,16 +231,23 @@ public sealed class CellularAndSminGradTests { attempts++; var p = new float3(random.NextSingle() * 6.0f - 3.0f, random.NextSingle() * 6.0f - 3.0f, random.NextSingle() * 6.0f - 3.0f); - var (f1, f2) = FindF1F2(p); + var (f1, f2, f3) = FindF1F2F3(p); // Exclude a band around the F1 = F2 set: the identity of the nearest and second-nearest // feature swap there, a genuine kink in F1 and F2 alike, not something central differences - // can approximate. Also skip the exact feature point itself (F1 = 0), covered separately. - if (f2 - f1 < 0.05f || f1 < 0.01f) + // can approximate. Also exclude a band around F2 = F3: edge (F2 - F1) kinks there too, + // since the second-nearest feature's identity swaps (Soul, cut 2a-i: an unfiltered sample + // hit a 0.78 edge error this way). Also skip the exact feature point itself (F1 = 0), + // covered separately. + if (f2 - f1 < 0.05f || f1 < 0.01f || f3 - f2 < 0.05f) continue; tested++; var result = math.cellular(p); + + // eps=0.0005f: unlike smin_grad, cellular's F1/F2 are Euclidean distance fields with real + // curvature, so central differences carry O(eps^2) truncation error and a larger eps is + // worse, not better (see CellularGradientTolerance's comment for the measured sweep). const float eps = 0.0005f; // Differentiate cellular's own returned .w values (F1): calling cellular again at p +/- @@ -221,8 +257,8 @@ public sealed class CellularAndSminGradTests float NearestValueAt(float3 q) => math.cellular(q).nearest.w; float EdgeValueAt(float3 q) => math.cellular(q).edge.w; - AssertGradientMatchesCentralDifference(NearestValueAt, p, eps, result.nearest); - AssertGradientMatchesCentralDifference(EdgeValueAt, p, eps, result.edge); + AssertGradientMatchesCentralDifference(NearestValueAt, p, eps, result.nearest, CellularGradientTolerance); + AssertGradientMatchesCentralDifference(EdgeValueAt, p, eps, result.edge, CellularGradientTolerance); } Assert.True(tested >= 100, $"only {tested} of {attempts} random points landed away from the F1=F2 band"); @@ -268,7 +304,7 @@ public sealed class CellularAndSminGradTests for (var i = 0; i < n; i++) { var cell = new float3(random.Next(-2000, 2000), random.Next(-2000, 2000), random.Next(-2000, 2000)); - var hash = math.pcg3d(cell); + var hash = math.pcg3d(math.int3(cell)); var jitter = new float3(Unit(hash.x), Unit(hash.y), Unit(hash.z)); ids[i] = math.cellular(cell + jitter).id; @@ -282,6 +318,60 @@ public sealed class CellularAndSminGradTests Assert.True(Math.Abs(Correlation(ids, jitterZ)) < 0.05, "id correlates with jitter.z"); } + [Fact] + public void JitterAndIdPassAChiSquareUniformityTest() + { + // F4/jitter uniformity: cellular hashes the winning cell's INTEGER coordinate (design.md, + // "Integer hashing"), never that coordinate's float32 bit pattern. Hashing the bit pattern of + // an integer-valued float shares more structure between adjacent cells (same exponent, a + // mantissa one increment apart) than the integers themselves do going into pcg3d's own + // mixing, and that structure survived into the jitter: Soul measured jitter.x's marginal + // chi-square statistic at 203.8 over 1,000,000 cells this way (15 degrees of freedom; + // critical value 37.7 at alpha = 0.001), soundly rejecting uniformity. 16 equal-width bins + // over [0, 1) give those same 15 degrees of freedom for jitter.x/y/z and id alike. + const int n = 200_000; + const int bins = 16; + const double criticalValue = 37.7; // chi-square, 15 df, alpha = 0.001 + + var jitterXCounts = new int[bins]; + var jitterYCounts = new int[bins]; + var jitterZCounts = new int[bins]; + var idCounts = new int[bins]; + var random = new System.Random(0xCE19); + + for (var i = 0; i < n; i++) + { + var cell = new float3(random.Next(-100000, 100000), random.Next(-100000, 100000), random.Next(-100000, 100000)); + var hash = math.pcg3d(math.int3(cell)); + jitterXCounts[Bin(Unit(hash.x), bins)]++; + jitterYCounts[Bin(Unit(hash.y), bins)]++; + jitterZCounts[Bin(Unit(hash.z), bins)]++; + + var idHash = math.pcg4d(new int4(math.int3(cell), 0)); + idCounts[Bin(Unit(idHash.w), bins)]++; + } + + var expected = (double)n / bins; + AssertUniform(jitterXCounts, expected, criticalValue, "jitter.x"); + AssertUniform(jitterYCounts, expected, criticalValue, "jitter.y"); + AssertUniform(jitterZCounts, expected, criticalValue, "jitter.z"); + AssertUniform(idCounts, expected, criticalValue, "id"); + } + + private static int Bin(float value, int binCount) => Math.Clamp((int)(value * binCount), 0, binCount - 1); + + private static void AssertUniform(int[] counts, double expected, double criticalValue, string name) + { + var chiSquare = 0.0; + foreach (var count in counts) + { + var delta = count - expected; + chiSquare += delta * delta / expected; + } + + Assert.True(chiSquare < criticalValue, $"{name} chi-square {chiSquare:G6} exceeds critical value {criticalValue} ({counts.Length - 1} df)"); + } + [Fact] public void NoNaNAtIntegerPointsNear2To24() { @@ -340,57 +430,80 @@ public sealed class CellularAndSminGradTests // F6, part 2 (Soul finding): a single-slice loop-bound mutant on one axis (e.g. dz < 2 // instead of dz <= 2, silently dropping the +2 layer on that axis alone) still passed the // random 20,000-point sweeps above. That is not those sweeps being weak; it is the - // exactness proof's own margin (F1, F2 always < sqrt(3), radius-2 cells' lower bound is - // exactly 2) making a genuine win from any single |offset component| = 2 cell astronomically - // rare over uniformly random points (zero hits in 3,000,000 samples, measured separately). + // exactness proof's own margin (F1, F2 always < sqrt(3), a radius-2 cell's lower bound is + // max(0, 2-1) = 1, not "exactly 2" as an earlier version of this comment claimed) making a + // genuine win from any single |offset component| = 2 cell rare over uniformly random points + // (Soul, cut 2a-i: measured about 1 hit in 28,000 samples, not the "zero in 3,000,000" an + // earlier version of this comment claimed). // Hunting for one by choosing candidate query POINTS is the wrong end of the search; this // instead scans candidate integer CELLS and, for each, builds the query point most likely to // realize that cell's neighbour as the winner: the closest point of the query cell's own unit // cube to the target neighbour's actual jittered feature. That construction finds a real // example for every one of the six axis-aligned radius-2 offsets within a couple thousand // candidate cells, deterministically (pcg3d has no randomness to get lucky or unlucky with). - var offsets = new[] + foreach (var offset in AxisAlignedRadiusTwoOffsets) { - new float3(2.0f, 0.0f, 0.0f), new float3(-2.0f, 0.0f, 0.0f), - new float3(0.0f, 2.0f, 0.0f), new float3(0.0f, -2.0f, 0.0f), - new float3(0.0f, 0.0f, 2.0f), new float3(0.0f, 0.0f, -2.0f), - }; - - foreach (var offset in offsets) - { - var found = false; - for (var qz = -80; qz <= 80 && !found; qz++) - for (var qy = -80; qy <= 80 && !found; qy++) - for (var qx = -80; qx <= 80 && !found; qx++) - { - var q = new float3(qx, qy, qz); - var targetFeature = FeaturePoint(q + offset); - var p = new float3( - Math.Clamp(targetFeature.x, q.x, q.x + 0.999f), - Math.Clamp(targetFeature.y, q.y, q.y + 0.999f), - Math.Clamp(targetFeature.z, q.z, q.z + 0.999f)); - - var (f1, f2) = FindF1F2(p); // independent 7x7x7 oracle - var d = math.distance(p, targetFeature); - if (MathF.Abs(d - f1) > 1.0e-4f && MathF.Abs(d - f2) > 1.0e-4f) - continue; // this cell didn't realize the target offset as F1 or F2; try another. - - found = true; - var result = math.cellular(p); - Assert.Equal(f1, result.nearest.w, precision: 3); - Assert.Equal(f2 - f1, result.edge.w, precision: 3); - } - - Assert.True(found, $"could not construct a query point where offset {offset.x},{offset.y},{offset.z} ever wins F1 or F2"); + var p = FindPointRealizingOffset(offset); + var (f1, f2) = FindF1F2(p); // independent 7x7x7 oracle + var result = math.cellular(p); + Assert.Equal(f1, result.nearest.w, precision: 3); + Assert.Equal(f2 - f1, result.edge.w, precision: 3); } } + // The six axis-aligned radius-2 offsets (F6, part 2): shared by + // SearchReachesEveryAxisAlignedRadiusTwoSlice and, via FindPointRealizingOffset, by + // HlslSourceCompatibilityTests' mirror-comparison inputs, so both pin the same construction. + internal static readonly float3[] AxisAlignedRadiusTwoOffsets = + { + new(2.0f, 0.0f, 0.0f), new(-2.0f, 0.0f, 0.0f), + new(0.0f, 2.0f, 0.0f), new(0.0f, -2.0f, 0.0f), + new(0.0f, 0.0f, 2.0f), new(0.0f, 0.0f, -2.0f), + }; + + // Hunting for a query point where a radius-2 offset wins F1 or F2 by choosing candidate query + // POINTS is the wrong end of the search (Soul, cut 2a-i: about 1 hit in 28,000 uniformly random + // points); this instead scans candidate integer CELLS and, for each, builds the query point most + // likely to realize that cell's neighbour as the winner: the closest point of the query cell's + // own unit cube to the target neighbour's actual jittered feature. That construction finds a + // real example for every one of the six axis-aligned radius-2 offsets within a couple thousand + // candidate cells, deterministically (pcg3d has no randomness to get lucky or unlucky with). + internal static float3 FindPointRealizingOffset(float3 offset) + { + for (var qz = -80; qz <= 80; qz++) + for (var qy = -80; qy <= 80; qy++) + for (var qx = -80; qx <= 80; qx++) + { + var q = new float3(qx, qy, qz); + var targetFeature = FeaturePoint(q + offset); + var p = new float3( + Math.Clamp(targetFeature.x, q.x, q.x + 0.999f), + Math.Clamp(targetFeature.y, q.y, q.y + 0.999f), + Math.Clamp(targetFeature.z, q.z, q.z + 0.999f)); + + var (f1, f2) = FindF1F2(p); // independent 7x7x7 oracle + var d = math.distance(p, targetFeature); + if (MathF.Abs(d - f1) > 1.0e-4f && MathF.Abs(d - f2) > 1.0e-4f) + continue; // this cell didn't realize the target offset as F1 or F2; try another. + + return p; + } + + throw new InvalidOperationException($"could not construct a query point where offset {offset.x},{offset.y},{offset.z} ever wins F1 or F2"); + } + [Fact] - public void PrunedSearchIsBitIdenticalToTheUnprunedFiveCubedSearch() + public void PrunedSearchIsBitIdenticalToTheUnprunedFiveCubedSearchBelowTwoToTheTwentyFifth() { // F6, part 2: the lower-bound prune inside math.cellular's loop must never change the // result. Compares against an independent, unpruned 5x5x5 reference that visits every cell - // in the same fixed order, over a large seeded sweep. + // in the same fixed order, over a large seeded sweep in [-50, 50]^3. + // + // This bit-identity has a domain: it holds for |p| < 2^25 (design.md, "Precision domain"). + // Above that, cell + offset starts to round in float32, so the pruned loop's recomputed + // neighbour and the unpruned loop's neighbour at the same nominal (dx, dy, dz) can land on + // different actual cells, and the two searches diverge (Soul, cut 2a-i: 180 of 1,100,000 + // points differ at |p| ~= 3.4e7). [-50, 50]^3 is comfortably inside the proven domain. var random = new System.Random(0xCE18); for (var i = 0; i < 20000; i++) { @@ -434,7 +547,7 @@ public sealed class CellularAndSminGradTests var grad1 = f1 > 0.0f ? (p - c1) / f1 : float3.zero; var grad2 = f2 > 0.0f ? (p - c2) / f2 : float3.zero; - var idHash = math.pcg4d(new float4(cellId, 0.0f)); + var idHash = math.pcg4d(new int4(math.int3(cellId), 0)); return new CultCellular(new float4(grad1, f1), new float4(grad2 - grad1, f2 - f1), Unit(idHash.w)); } @@ -468,9 +581,9 @@ public sealed class CellularAndSminGradTests // 3), strictly wider than production's own proven-sufficient 5x5x5 (F6), so it is a true // independent oracle rather than a second copy of the algorithm under test. ---- - private static float3 FeaturePoint(float3 cell) + internal static float3 FeaturePoint(float3 cell) { - var hash = math.pcg3d(cell); + var hash = math.pcg3d(math.int3(cell)); var jitter = new float3(Unit(hash.x), Unit(hash.y), Unit(hash.z)); return cell + jitter; } @@ -480,20 +593,34 @@ public sealed class CellularAndSminGradTests private static float Unit(int bits) => ((uint)bits >> 8) * (1.0f / 16777216.0f); private static (float f1, float f2) FindF1F2(float3 p) + { + var (f1, f2, _) = FindF1F2F3(p); + return (f1, f2); + } + + // Same independent 7x7x7 oracle as FindF1F2, plus the third-nearest distance. Central-difference + // tests need this: excluding only the F1=F2 gap is not enough, because a small F2=F3 gap is also + // a kink (the identity of the second-nearest feature swaps there), and hitting one produces a + // large, curvature-driven finite-difference error unrelated to eps or float noise (Soul, cut + // 2a-i: an unfiltered large sample hit an edge error of 0.78 this way, versus ~2e-4 once F2=F3 + // is also excluded). + private static (float f1, float f2, float f3) FindF1F2F3(float3 p) { var cell = new float3(MathF.Floor(p.x), MathF.Floor(p.y), MathF.Floor(p.z)); var f1 = float.PositiveInfinity; var f2 = float.PositiveInfinity; + var f3 = float.PositiveInfinity; for (var dz = -3; dz <= 3; dz++) for (var dy = -3; dy <= 3; dy++) for (var dx = -3; dx <= 3; dx++) { var d = math.distance(p, FeaturePoint(cell + new float3(dx, dy, dz))); - if (d < f1) { f2 = f1; f1 = d; } - else if (d < f2) { f2 = d; } + if (d < f1) { f3 = f2; f2 = f1; f1 = d; } + else if (d < f2) { f3 = f2; f2 = d; } + else if (d < f3) { f3 = d; } } - return (f1, f2); + return (f1, f2, f3); } private static (float3 offset, float f1, float f2) FindNearestOffset(float3 p) @@ -532,16 +659,16 @@ public sealed class CellularAndSminGradTests return bestCell; } - private static void AssertGradientMatchesCentralDifference(Func scalarField, float3 p, float eps, float4 valueAndGradient) + private static void AssertGradientMatchesCentralDifference(Func scalarField, float3 p, float eps, float4 valueAndGradient, float tolerance) { var numeric = new float3( (scalarField(p + new float3(eps, 0.0f, 0.0f)) - scalarField(p - new float3(eps, 0.0f, 0.0f))) / (2.0f * eps), (scalarField(p + new float3(0.0f, eps, 0.0f)) - scalarField(p - new float3(0.0f, eps, 0.0f))) / (2.0f * eps), (scalarField(p + new float3(0.0f, 0.0f, eps)) - scalarField(p - new float3(0.0f, 0.0f, eps))) / (2.0f * eps)); - AssertWithinTolerance(numeric.x, valueAndGradient.x, GradientTolerance); - AssertWithinTolerance(numeric.y, valueAndGradient.y, GradientTolerance); - AssertWithinTolerance(numeric.z, valueAndGradient.z, GradientTolerance); + AssertWithinTolerance(numeric.x, valueAndGradient.x, tolerance); + AssertWithinTolerance(numeric.y, valueAndGradient.y, tolerance); + AssertWithinTolerance(numeric.z, valueAndGradient.z, tolerance); } // A fixed decimal-place Assert.Equal snaps two numbers that agree to within a few times 1e-5 to