diff --git a/packages/cultmath/docs/design.md b/packages/cultmath/docs/design.md index 15cb34f8..38a14215 100644 --- a/packages/cultmath/docs/design.md +++ b/packages/cultmath/docs/design.md @@ -745,8 +745,14 @@ refined 50.9549, witness (-245.19885, -53.01323, -91.47363). `SNOISE2_LIPSCHITZ = 8.117002` is the same for the gradient of `snoise(float2)` (`MeasureLipschitz2`: sampled 7.2575, refined 7.3791, witness (-67.12454, -129.53296), `|snoise2|` at most 0.99964), and it backs `iv_snoise_ball(float2, -float)`. Both are empirical with a margin, not proofs; the enclosure tests are -the defence, and a change to the kernels must re-run the measurements. +float)`. Both are empirical with a margin, not proofs. A random ball cannot see a +curvature constant 20% low, so `SNOISE_HESSIAN` is defended by +`HessianConstantPinsSampledCurvature` and `SnoiseFormEnclosesTheCurvatureWitness` +(`af_snoise` at the point where the curvature is realised, along 4,000 directions: a +constant of 40 misses there). The margin over the function itself: a scan of the whole +289^3 hash period found a realised maximum of 50.58, and bounding each corner's gradient +independently over every hash configuration gives 52.22, so 56.05 sits 7% above the +configuration bound. A change to the kernels must re-run the measurements. The slice. `af_frustum_ball(centreSlope, z0, z1, footprintPerDepth, warpVariation)` returns the centre `(m_c z_m, z_m)` and a radius `z1 diff --git a/packages/cultmath/src/CultMath/math.Affine.cs b/packages/cultmath/src/CultMath/math.Affine.cs index 6b392683..85228e01 100644 --- a/packages/cultmath/src/CultMath/math.Affine.cs +++ b/packages/cultmath/src/CultMath/math.Affine.cs @@ -25,8 +25,13 @@ public static partial class math // differences of snoise_grad, largest absolute eigenvalue) over 1e6 seeded points, refines the 1e4 // largest by ascent, and multiplies the refined maximum by 1.10; // NoiseBoundTests.HessianConstantPinsSampledCurvature pins this value against a fast re-measure and - // the ascent from MeasureHessian's witness. Empirical with a margin, not a proof; the enclosure - // tests are the defence. A change to the snoise kernel must re-run MeasureHessian and re-pin it. + // the ascent from MeasureHessian's witness. Empirical with a margin, not a proof. Random balls + // cannot see a constant 20% low, so the defence is those two curvature tests plus + // NoiseBoundTests.SnoiseFormEnclosesTheCurvatureWitness, where af_snoise misses at the point the + // curvature is realised if the constant is under it. The margin over the function: a scan of the + // whole 289^3 hash period found a realised maximum of 50.58, and bounding each corner's gradient + // independently over every hash configuration gives 52.22 (Soul pass s1). A change to the snoise + // kernel must re-run MeasureHessian and re-pin it. public const float SNOISE_HESSIAN = 56.050385f; /// The point form (x, 0, 0): the value x at every eps. Exact. diff --git a/packages/cultmath/tests/CultMath.Tests/AffineTests.cs b/packages/cultmath/tests/CultMath.Tests/AffineTests.cs index 482b5ecf..612de682 100644 --- a/packages/cultmath/tests/CultMath.Tests/AffineTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/AffineTests.cs @@ -612,4 +612,62 @@ public sealed class AffineTests Assert.True(misses.Count == 0, $"{misses.Count} misses:" + Environment.NewLine + string.Join(Environment.NewLine, misses.Take(8))); } + + /// + /// af_snoise's width, pinned both ways. Over 6,000 seeded (centre, axis, radius) with reach R from 1e-3 + /// to 3 (centres in [-20, 20]^3, a third out to [-1000, 1000]^3, a third of the axes and a third of the + /// radii 0), the form is the documented rule evaluated in double from snoise_grad: x0 = n; the centred + /// form when |g . axis| + |g| radius + M R^2 / 2 is under L R (checked away from that boundary), else + /// the Lipschitz form; e = that half-width + the float32 allowance (2^-14 + 2^-12 R + L 2^-20 (|centre|_1 + R)) + /// plus the form's own rounding absorption (|n| + |a| + e) 2^-20; each of x0, a and e within a part in + /// 1e6 of it. The enclosure tests pass on an allowance or an absorption of any size on the CPU, where + /// float32 is exact; this is the test that sees them. + /// + [Fact] + public void SnoiseFormCarriesItsDocumentedWidth() + { + var u20 = Math.ScaleB(1.0, -20); + double lipschitz = SNOISE_LIPSCHITZ, hessian = SNOISE_HESSIAN; + var random = new System.Random(0xAF13); + var (centred, wide) = (0, 0); + for (var t = 0; t < 6_000; t++) + { + var spread = t % 3 == 2 ? 1000.0f : 20.0f; + var centre = new float3(Uniform(random, -spread, spread), Uniform(random, -spread, spread), Uniform(random, -spread, spread)); + var axis = random.Next(3) == 0 ? new float3(0.0f, 0.0f, 0.0f) : UnitVector(random) * LogUniform(random, 1.0e-3f, 3.0f); + var radius = random.Next(3) == 0 ? 0.0f : LogUniform(random, 1.0e-3f, 3.0f); + var g = snoise_grad(centre); + var reach = (double)length(axis) + radius; + var lip = lipschitz * reach; + var a = (double)g.x * axis.x + (double)g.y * axis.y + (double)g.z * axis.z; + var remainder = Math.Sqrt((double)g.x * g.x + (double)g.y * g.y + (double)g.z * g.z) * radius + 0.5 * hessian * reach * reach; + var allowance = Math.ScaleB(1.0, -14) + Math.ScaleB(reach, -12) + + lipschitz * u20 * (Math.Abs((double)centre.x) + Math.Abs((double)centre.y) + Math.Abs((double)centre.z) + reach); + + double Width(bool useCentred, out double slope) + { + slope = useCentred ? a : 0.0; + var e = (useCentred ? remainder : lip) + allowance; + return e + (Math.Abs((double)g.w) + Math.Abs(slope) + e) * u20; + } + + var zCentred = Width(true, out var slopeCentred); + var zLipschitz = Width(false, out _); + var form = af_snoise(centre, axis, radius); + var observedCentred = Math.Abs(form.z - zCentred) <= Math.Abs(form.z - zLipschitz); + var expectedCentred = Math.Abs(a) + remainder < lip; + if (Math.Abs(Math.Abs(a) + remainder - lip) > 1.0e-3 * lip) + Assert.True(observedCentred == expectedCentred || zCentred == zLipschitz, $"af_snoise({centre}, {axis}, {radius:R}) = {form} took the {(observedCentred ? "centred" : "Lipschitz")} form; the rule takes the {(expectedCentred ? "centred" : "Lipschitz")} one (|a| + remainder = {Math.Abs(a) + remainder:R}, L R = {lip:R})"); + var z = observedCentred ? zCentred : zLipschitz; + var slope = observedCentred ? slopeCentred : 0.0; + Assert.True(form.x == g.w, $"af_snoise({centre}, {axis}, {radius:R}).x0 = {form.x:R}, snoise = {g.w:R}"); + Assert.True(Math.Abs(form.y - slope) <= 1.0e-5 * Math.Abs((double)length(new float3(g.x, g.y, g.z)) * length(axis)) + 1.0e-12, $"af_snoise({centre}, {axis}, {radius:R}).a = {form.y:R}, g . axis = {slope:R}"); + Assert.True(Math.Abs(form.z - z) <= 1.0e-6 * z, $"af_snoise({centre}, {axis}, {radius:R}) = {form}: e is {form.z:R}, the documented rule gives {z:R}"); + centred += observedCentred ? 1 : 0; + wide += !observedCentred ? 1 : 0; + } + + output.WriteLine($"AF-REPORT af_snoise width: {centred} centred cases, {wide} Lipschitz cases"); + Assert.True(centred > 1000 && wide > 1000, $"the cases do not exercise both forms: {centred} centred, {wide} Lipschitz"); + } } diff --git a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs index cc7caebb..5fb3672a 100644 --- a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs @@ -535,7 +535,9 @@ public sealed class NoiseBoundTests // Hill-climbs f from p as AscendGradientNorm does: the ascent direction is f's central difference, // the step grows on success and halves on failure. Returns the largest value reached. - private static double Ascend(Func f, float3 p) + private static double Ascend(Func f, float3 p) => AscendPoint(f, p).Best; + + private static (double Best, float3 At) AscendPoint(Func f, float3 p) { const float h = 1.0e-3f; var step = 1.0e-2f; @@ -562,7 +564,7 @@ public sealed class NoiseBoundTests } } - return best; + return (best, p); } private static double Ascend2(Func f, float2 p) => Ascend(q => f(new float2(q.x, q.y)), new float3(p.x, p.y, 0.0f)); @@ -633,6 +635,39 @@ public sealed class NoiseBoundTests Assert.True(Math.Abs(ratio - 1.0) <= 0.005, $"SNOISE_HESSIAN is {ratio:R} times 1.10 x the re-measured {refined:R}; re-run MeasureHessian and re-pin"); } + /// + /// SNOISE_HESSIAN where the curvature is realised. At the point the ascent from MeasureHessian's witness + /// reaches, the absolute second difference of snoise along 4,000 seeded directions u at step 0.05 has a + /// largest value within 10% of SNOISE_HESSIAN / 1.10 (the realised curvature is there, so the check + /// has teeth), and af_snoise(witness, r u, 0) encloses snoise(witness + eps r u) at eps in + /// {-1, -1/2, 0, 1/2, 1} for r in {0.1, 0.3} and every one of those directions, no tolerance: a constant + /// under the realised curvature (40 against about 50.6) misses on the directions near the top + /// eigenvector. The enclosure tests on random balls cannot see a constant 20% low; this and + /// HessianConstantPinsSampledCurvature are its defence. + /// + [Fact] + public void SnoiseFormEnclosesTheCurvatureWitness() + { + var (_, witness) = AscendPoint(HessianNorm, HessianWitness); + var random = new System.Random(0xCA7E); + var directions = Enumerable.Range(0, 4000).Select(_ => UnitVector(random)).ToArray(); + const float step = 0.05f; + var centreValue = snoise(witness); + var realised = directions.Max(u => Math.Abs(snoise(witness + u * step) + (double)snoise(witness - u * step) - 2.0 * centreValue) / ((double)step * step)); + output.WriteLine($"IV-REPORT curvature witness ({witness.x:R}, {witness.y:R}, {witness.z:R}): largest second difference {realised:R} against SNOISE_HESSIAN {SNOISE_HESSIAN:R}"); + Assert.True(realised >= 0.9 * SNOISE_HESSIAN / 1.10, $"the largest second difference at the witness is {realised:R}, under 0.9 SNOISE_HESSIAN / 1.10 = {0.9 * SNOISE_HESSIAN / 1.10:R}: the check has no teeth, or the constant moved"); + foreach (var u in directions) + foreach (var r in new[] { 0.1f, 0.3f }) + { + var form = af_snoise(witness, u * r, 0.0f); + foreach (var eps in new[] { -1.0f, -0.5f, 0.0f, 0.5f, 1.0f }) + { + var value = snoise(witness + u * (r * eps)); + Assert.True(Math.Abs(value - (form.x + (double)form.y * eps)) <= form.z, $"af_snoise(witness, {u * r}, 0) = {form} misses snoise = {value:R} at eps {eps:R}, direction {u}"); + } + } + } + /// /// The provenance of SNOISE2_LIPSCHITZ: the largest gradient length of snoise(float2) (Gradient2Norm) /// over 1e6 seeded points of [-256, 256]^2, the 1e4 largest refined by ascent, times 1.10; and the