diff --git a/packages/cultmath/docs/design.md b/packages/cultmath/docs/design.md index 3511d755..15a28cba 100644 --- a/packages/cultmath/docs/design.md +++ b/packages/cultmath/docs/design.md @@ -243,6 +243,13 @@ the main families' reach is at least about 1e-3, where the second-order remainde allowance still passes. `GoldenFixtureMatchesCSharp` requires each of the two to hold cases that fail on C#'s own results once the allowance is removed. +Four more appended families, the coordinate-band `snoise_grad` and `snoise_hessian` reads (256 cases each, the magnitude of the +point log-uniform from 1e-4 to 1e2 and from 1e2 to 4e6), carry an allowance `check`: with `S = |x|_1` and the allowances +`math.span_slack(math.span_rounding(1, x, zero stage, 0, 0, 0))` (the owners of the span stage's allowances in C#, HLSL and GLSL), +the result must lie within the allowance of the kernel's exact value, gradient or Hessian in double, no tolerance. The two small +families must also fail it in at least 16 cases once the absolute term is removed. `GoldenFixtureMatchesCSharp` applies the check +to C#'s own results, and a device capture is judged by the same `Verify`. + `GoldenFixtureMatchesCSharp` evaluates C# on the committed fixture's own arguments. Arguments are drawn only when the fixture is regenerated (`CULTMATH_WRITE_GLSL=1`), because the draw itself calls `MathF.Pow`, whose diff --git a/packages/cultmath/tests/CultMath.Tests/GlslMirrorTests.cs b/packages/cultmath/tests/CultMath.Tests/GlslMirrorTests.cs index f643d2c3..982e49c3 100644 --- a/packages/cultmath/tests/CultMath.Tests/GlslMirrorTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/GlslMirrorTests.cs @@ -408,6 +408,74 @@ static int4 clamp(int4 x, int4 a, int4 b) { if (a.x > b.x || a.y > b.y || a.z > return $"least enclosure slack {enclosure:E3} of the radius, least slack above R0 + d {lower:E3} of R0, below its widening {upper:F3} of 2^-19 (R0 + |c0|_1 + spread)"; } + private const string SpanAllowanceRead = + "Read at the case's point x, with S = |x|_1 and (e_n, e_g, e_h) = math.span_slack(math.span_rounding(1, x, zero stage, 0, 0, 0)), " + + "against the kernel's exact values in double: |n - n_exact| <= e_n and |g - g_exact|_2 <= e_g for snoise_grad, " + + "||H - H_exact||_F <= e_h for snoise_hessian, no tolerance. " + + "Report the largest fraction of each allowance per decade of S and the largest absolute error at S < 1e-2."; + + // SpanAllowanceRead applied to a family's results (C#'s own, or a device's): every read is within the + // allowance math.span_slack gives at its point, the owner's own numbers and no restated literal. A small + // family (bands 1e-4 to 1e2) must also fail the check once the absolute term is removed in at least 16 + // cases, and hold at least 48 cases with S <= 1e-2: the allowance's absolute term is what decides there. + private static string AssertSpanAllowance(JsonArray cases, bool hessian, bool small) + { + var absolute = math.span_slack(0.0f); + var fractions = new SortedDictionary(); + double smallestError = 0.0; + int absoluteDecides = 0, nearOrigin = 0; + foreach (var entry in cases) + { + var sides = ((string)entry!).Split(" -> "); + double[] x = Floats(sides[0]), result = Floats(sides[1]); + var s = Math.Abs(x[0]) + Math.Abs(x[1]) + Math.Abs(x[2]); + var slack = math.span_slack(math.span_rounding(1.0f, new float3((float)x[0], (float)x[1], (float)x[2]), default, new float3(0.0f, 0.0f, 0.0f), 0.0f, 0.0f)); + double valueError = 0.0, gradientError, valueFraction = 0.0, gradientFraction; + bool needsAbsolute; + if (hessian) + { + var exact = NoiseBoundTests.HessianDouble(x[0], x[1], x[2]); + var sum = 0.0; + for (var i = 0; i < 9; i++) + sum += (result[i] - exact[i]) * (result[i] - exact[i]); + gradientError = Math.Sqrt(sum); + Assert.True(gradientError <= slack.z, $"snoise_hessian case {entry}: ||H - H_exact||_F is {gradientError:R}, past the allowance {slack.z:R}"); + gradientFraction = gradientError / slack.z; + needsAbsolute = gradientError > slack.z - absolute.z; + } + else + { + var exact = NoiseBoundTests.GradDouble(x[0], x[1], x[2]); + valueError = Math.Abs(result[3] - NoiseBoundTests.ValueDouble(x[0], x[1], x[2])); + gradientError = Math.Sqrt((result[0] - exact[0]) * (result[0] - exact[0]) + (result[1] - exact[1]) * (result[1] - exact[1]) + (result[2] - exact[2]) * (result[2] - exact[2])); + Assert.True(valueError <= slack.x, $"snoise_grad case {entry}: |n - n_exact| is {valueError:R}, past the allowance {slack.x:R}"); + Assert.True(gradientError <= slack.y, $"snoise_grad case {entry}: |g - g_exact|_2 is {gradientError:R}, past the allowance {slack.y:R}"); + valueFraction = valueError / slack.x; + gradientFraction = gradientError / slack.y; + needsAbsolute = valueError > slack.x - absolute.x || gradientError > slack.y - absolute.y; + } + + var decade = (int)Math.Floor(Math.Log10(s)); + fractions.TryGetValue(decade, out var largest); + fractions[decade] = (Math.Max(largest.Value, valueFraction), Math.Max(largest.Gradient, gradientFraction)); + if (needsAbsolute) + absoluteDecides++; + if (s < 1.0e-2) + smallestError = Math.Max(smallestError, Math.Max(valueError, gradientError)); + if (s <= 1.0e-2) + nearOrigin++; + } + + if (small) + { + Assert.True(absoluteDecides >= 16, $"only {absoluteDecides} of {cases.Count} cases fail once the absolute term is removed: the family cannot tell a device that needs it from one that does not"); + Assert.True(nearOrigin >= 48, $"only {nearOrigin} of {cases.Count} cases have |x|_1 <= 1e-2"); + } + + var report = string.Join(", ", fractions.Select(f => hessian ? $"1e{f.Key}: {f.Value.Gradient:F3}" : $"1e{f.Key}: {f.Value.Value:F3}/{f.Value.Gradient:F3}")); + return $"largest fraction of {(hessian ? "e_h" : "e_n/e_g")} per decade of |x|_1: {report}; largest absolute error at |x|_1 < 1e-2 {smallestError:E3}; {absoluteDecides} of {cases.Count} cases need the absolute term"; + } + // A committed case's arguments, decoded into the C# method's parameter types in Hex's field order. private static object[] Arguments(MethodInfo method, string hex) { @@ -575,7 +643,8 @@ static int4 clamp(int4 x, int4 a, int4 b) { if (a.x > b.x || a.y > b.y || a.z > // Every family the site's WebGL2 readback evaluates, in a fixed order. A new family is appended, so the // seeded stream and every earlier family's cases stay as they were. pcg3d and pcg4d are "exact"; the // float families are ulp-bounded, with the bound measured on the device by the consumer, and - // iv_frustum_ball also carries an enclosure check. phacelle is in CultMath.Phacelle.glsl, so its + // iv_frustum_ball also carries an enclosure check; the four coordinate-band snoise_grad and snoise_hessian + // families carry an allowance check, applied to C#'s own results here and to a device's by the device judge. phacelle is in CultMath.Phacelle.glsl, so its // evaluator concatenates that file after CultMath.glsl; the site, which does not vendor it, skips it. private static readonly Family[] Families = new Family[] { @@ -720,8 +789,24 @@ static int4 clamp(int4 x, int4 a, int4 b) { if (a.x > b.x || a.y > b.y || a.z > // midpoint a Point3, one span in four a point. new Family("SpanNoise cultmath_span_fbm(SpanPath, SpanStage, int, float, float)", "ulp-bounded", SpanArguments), new Family("SpanNoise cultmath_span_ridged(SpanPath, SpanStage, int, float, float)", "ulp-bounded", SpanArguments), + // The kernel's reads over coordinate bands, the magnitude log-uniform in the band: the allowance of + // math.span_slack is checked against the exact kernel at every case (SpanAllowanceRead). + new Family("vec4 cultmath_snoise_grad(vec3)", "ulp-bounded, coordinate bands 1e-4 to 1e2", r => new object[] { BandPoint3(r, 1.0e-4, 1.0e2) }, SpanAllowanceRead, Verify: cases => AssertSpanAllowance(cases, hessian: false, small: true)), + new Family("vec4 cultmath_snoise_grad(vec3)", "ulp-bounded, coordinate bands 1e2 to 4e6", r => new object[] { BandPoint3(r, 1.0e2, 4.0e6) }, SpanAllowanceRead, Verify: cases => AssertSpanAllowance(cases, hessian: false, small: false)), + new Family("mat3 cultmath_snoise_hessian(vec3)", "ulp-bounded, coordinate bands 1e-4 to 1e2", r => new object[] { BandPoint3(r, 1.0e-4, 1.0e2) }, SpanAllowanceRead, Verify: cases => AssertSpanAllowance(cases, hessian: true, small: true)), + new Family("mat3 cultmath_snoise_hessian(vec3)", "ulp-bounded, coordinate bands 1e2 to 4e6", r => new object[] { BandPoint3(r, 1.0e2, 4.0e6) }, SpanAllowanceRead, Verify: cases => AssertSpanAllowance(cases, hessian: true, small: false)), }).ToArray(); + // A point of magnitude m, m log-uniform in [lo, hi]: one component at +-m, the other two uniform in [-m, m], + // the axis of the large one random. + private static float3 BandPoint3(System.Random r, double lo, double hi) + { + var m = lo * Math.Pow(hi / lo, r.NextDouble()); + var v = new double[] { (r.NextDouble() * 2.0 - 1.0) * m, (r.NextDouble() * 2.0 - 1.0) * m, (r.NextDouble() * 2.0 - 1.0) * m }; + v[r.Next(3)] = r.Next(2) == 0 ? m : -m; + return new float3((float)v[0], (float)v[1], (float)v[2]); + } + private static object[] SpanArguments(System.Random r) { var stage = SpanStageTests.Stage.Draw(r); diff --git a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs index f97d5bfe..643c994a 100644 --- a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs @@ -722,6 +722,44 @@ public sealed class NoiseBoundTests return total; } + // The gradient of snoise at v in double: 105 sum (m^4 p - 8 m^3 s x) over Corners, the derivative of + // ValueDouble (every corner's x moves one for one with v, dm = -2 x, ds = p). + internal static double[] GradDouble(double vx, double vy, double vz) + { + var (x, p, _) = Corners(vx, vy, vz); + var gradient = new double[3]; + for (var c = 0; c < 4; c++) + { + var m = 0.5 - (x[c][0] * x[c][0] + x[c][1] * x[c][1] + x[c][2] * x[c][2]); + if (m <= 0.0) + continue; + var s = p[c][0] * x[c][0] + p[c][1] * x[c][1] + p[c][2] * x[c][2]; + for (var i = 0; i < 3; i++) + gradient[i] += 105.0 * (m * m * m * m * p[c][i] - 8.0 * m * m * m * s * x[c][i]); + } + + return gradient; + } + + /// + /// GradDouble is the kernel's gradient: within float rounding of snoise_grad (L1 distance at most 2e-4) + /// at 300 seeded points of [-20, 20]^3, the way ExactReferenceIsTheKernel pins ValueDouble. + /// + [Fact] + public void GradDoubleIsTheKernelsGradient() + { + var random = new System.Random(0x6BAD); + double Coordinate() => random.NextDouble() * 40.0 - 20.0; + for (var k = 0; k < 300; k++) + { + var p = new float3((float)Coordinate(), (float)Coordinate(), (float)Coordinate()); + var exact = GradDouble(p.x, p.y, p.z); + var kernel = snoise_grad(p); + var distance = Math.Abs(exact[0] - kernel.x) + Math.Abs(exact[1] - kernel.y) + Math.Abs(exact[2] - kernel.z); + Assert.True(distance <= 2.0e-4, $"GradDouble({p}) is {exact[0]:R}, {exact[1]:R}, {exact[2]:R}, snoise_grad {kernel.x:R}, {kernel.y:R}, {kernel.z:R}"); + } + } + private static double Delta(int a, int b) => a == b ? 1.0 : 0.0; // The symmetric third-derivative tensor of snoise at p, entry (i * 3 + j) * 3 + k, in closed form: diff --git a/packages/cultmath/tests/CultMath.Tests/SpanStageTests.cs b/packages/cultmath/tests/CultMath.Tests/SpanStageTests.cs index 8324b5a1..dc1784f9 100644 --- a/packages/cultmath/tests/CultMath.Tests/SpanStageTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/SpanStageTests.cs @@ -752,6 +752,54 @@ public sealed class SpanStageTests Assert.True(crossing >= 100, $"only {crossing} of 600 midpoints had an exact sign opposite to the float one: the check has no teeth"); } + /// + /// The allowance a span computes covers the float kernel's read at every point of its path: the identity + /// stage, the midpoint within |q|_inf 1, a velocity of random direction and speed 1 to 1e3 (log-uniform), + /// a halfLength 1 to 1e3, so path points lie out to about 1e6. At 16 points of the straight path + /// q + v t the value, the gradient's length and the Hessian's Frobenius norm of the float kernel are + /// within math.span_slack(math.span_rounding(...)) of the exact kernel's, the reach term of S_o included. + /// Teeth: the largest value error exceeds the allowance computed with nu = 0 by more than 10 times. + /// + [Fact] + public void SpanAllowanceCoversReadsAlongTheSpan() + { + var random = new System.Random(0x5EA4); + double worstValue = 0.0, worstGradient = 0.0, worstHessian = 0.0, teeth = 0.0; + for (var k = 0; k < 200; k++) + { + var q = Vector(random, -1.0f, 1.0f); + var speed = LogUniform(random, 1.0f, 1.0e3f); + var velocity = UnitVector(random) * speed; + var halfLength = LogUniform(random, 1.0f, 1.0e3f); + var slack = span_slack(span_rounding(1.0f, q, SpanStage.Identity, None, speed, halfLength)); + var midpointSlack = span_slack(span_rounding(1.0f, q, SpanStage.Identity, None, 0.0f, halfLength)); + for (var j = 0; j < 16; j++) + { + var x = q + velocity * (halfLength * (j / 7.5f - 1.0f)); + var value = NoiseBoundTests.ValueDouble(x.x, x.y, x.z); + var gradient = NoiseBoundTests.GradDouble(x.x, x.y, x.z); + var hessian = NoiseBoundTests.HessianDouble(x.x, x.y, x.z); + var read = snoise_grad(x); + var h = snoise_hessian(x); + var valueError = Math.Abs(read.w - value); + var gradientError = Math.Sqrt(Math.Pow(read.x - gradient[0], 2) + Math.Pow(read.y - gradient[1], 2) + Math.Pow(read.z - gradient[2], 2)); + var hessianError = Math.Sqrt(Math.Pow(h._m00 - hessian[0], 2) + Math.Pow(h._m01 - hessian[1], 2) + Math.Pow(h._m02 - hessian[2], 2) + + Math.Pow(h._m10 - hessian[3], 2) + Math.Pow(h._m11 - hessian[4], 2) + Math.Pow(h._m12 - hessian[5], 2) + + Math.Pow(h._m20 - hessian[6], 2) + Math.Pow(h._m21 - hessian[7], 2) + Math.Pow(h._m22 - hessian[8], 2)); + Assert.True(valueError <= slack.x, $"midpoint {q}, speed {speed:R}, halfLength {halfLength:R}, point {x}: the value error {valueError:R} is past the allowance {slack.x:R}"); + Assert.True(gradientError <= slack.y, $"midpoint {q}, speed {speed:R}, halfLength {halfLength:R}, point {x}: the gradient error {gradientError:R} is past the allowance {slack.y:R}"); + Assert.True(hessianError <= slack.z, $"midpoint {q}, speed {speed:R}, halfLength {halfLength:R}, point {x}: the Hessian error {hessianError:R} is past the allowance {slack.z:R}"); + worstValue = Math.Max(worstValue, valueError / slack.x); + worstGradient = Math.Max(worstGradient, gradientError / slack.y); + worstHessian = Math.Max(worstHessian, hessianError / slack.z); + teeth = Math.Max(teeth, valueError / midpointSlack.x); + } + } + + output.WriteLine($"IV-REPORT span allowance: largest fraction of e_n {worstValue:F3}, e_g {worstGradient:F3}, e_h {worstHessian:F3}; largest value error over the allowance without the reach {teeth:F1}"); + Assert.True(teeth > 10.0, $"the largest value error is only {teeth:F1} times the allowance without the reach: the check has no teeth"); + } + /// /// The fold's acceleration term, alone, at the identity stage and under a transform (frequency 2.5, stretch /// 0.5, shift and offset): the midpoint velocity is zero and the acceleration points along the gradient