diff --git a/packages/cultmath/README.md b/packages/cultmath/README.md index 4f2018a3..61351da3 100644 --- a/packages/cultmath/README.md +++ b/packages/cultmath/README.md @@ -82,8 +82,9 @@ shader semantics. the analytic Hessian of `snoise(float3)` on the same lattice, hash and tie order as `snoise_grad` (a symmetric `float3x3`; GLSL `mat3`), and `integral_clamp_linear(a, b, lo, hi, t0, t1)`, the closed-form integral over `[t0, t1]` of `clamp(a + b t, - lo, hi)`; and, in C# for now (shaders follow), the octave stages `span_fbm` and - `span_ridged`, which return a `SpanNoise` (`field` as `float4(gradient, value)`, + lo, hi)`; and, in C# for now (shaders follow), the octave stages `span_fbm(path, stage, + octaves, lacunarity, gain)` and `span_ridged`, over a `SpanPath` and a `SpanStage` (`SpanStage.Identity` + is the plain fbm), which return a `SpanNoise` (`field` as `float4(gradient, value)`, `curvature`, `fold`) for a stretched, offset, shifted octave sum over a span. Three measured bounds describe `snoise`: `SNOISE_LIPSCHITZ` (gradient length, `NoiseBoundTests.MeasureLipschitz`), `SNOISE_HESSIAN` (Hessian spectral norm, `MeasureHessian`) and `SNOISE_THIRD` (third-derivative operator norm, diff --git a/packages/cultmath/docs/design.md b/packages/cultmath/docs/design.md index 0d7cd738..b7edc78b 100644 --- a/packages/cultmath/docs/design.md +++ b/packages/cultmath/docs/design.md @@ -809,34 +809,62 @@ configuration bound. A change to the kernels must re-run the measurements. ### Span stages `span_fbm` and `span_ridged` (`math.Span.cs`, C# only so far) are the octave stages -of a warp chain, each returning one `SpanNoise { float4 field; float curvature; float -fold; }`. The span is the path `q(t)`, `|t| <= halfLength`, through `q` with the given -`velocity` at its midpoint, `|q'| <= speed` and `|q''| <= acceleration` along it. -Octave `o` samples `snoise_grad` at `x_o(q) = f_o (M q + shift) + offset + o -octaveStep`, with `f_o = lacunarity^o`, amplitude `gain^o`, and `M = frequency (I + -(1/stretch - 1) a a^T)`, `a = normalize(axis)`; a zero axis or `stretch` 1 means no -stretch, and `stretch` must be positive and finite (unguarded, like `lacunarity`). -`field` is the octave sum and its gradient in `q` (`f_o M grad n_o`); at frequency 1, -stretch 1 and zero shift, offset and octaveStep it equals `fbm_grad` and -`ridged_grad` bit for bit. `span_ridged` sums `gain^o (1 - |n_o|)`, `ridged_grad`'s -convention; the Span Lab's `1 - 2|n|` is twice this minus the gain sum, which the -consumer forms. +of a warp chain, `span_fbm(path, stage, octaves, lacunarity, gain)`, each returning +one `SpanNoise { float4 field; float curvature; float fold; }`. The inputs are two +structs with named public fields and no positional constructor, so the three floats +of the span and the six fields of the transform cannot be swapped silently. +`SpanPath { midpoint, velocity, speed, acceleration, halfLength }` is the path +`q(t)`, `|t| <= halfLength`, through `midpoint` with `velocity` at `t = 0`, +`|q'| <= speed` and `|q''| <= acceleration` along it. `SpanStage { frequency, +stretch, axis, shift, offset, octaveStep }` makes octave `o` sample `snoise_grad` +at `x_o(q) = f_o (M q + shift) + offset + o octaveStep`, with `f_o = +lacunarity^o`, amplitude `gain^o`, and `M = frequency (I + (1/stretch - 1) a a^T)`, +`a = normalize(axis)`; a zero axis or `stretch` 1 means no stretch, and `stretch` +must be positive and finite (unguarded, like `lacunarity`). `default(SpanStage)` +has frequency 0; `SpanStage.Identity` is the stage whose field equals `fbm_grad` +and `ridged_grad` bit for bit. `octaves`, `lacunarity` and `gain` stay positional: +they are `fbm_grad`'s own trailing triple, and `octaves` is the one int among the +scalars. `field` is the octave sum and its gradient in `q` (`f_o M grad n_o`). +`span_ridged` sums `gain^o (1 - |n_o|)`, `ridged_grad`'s convention. The Span +Lab's stage is `1 - 2|n|` per octave: with `G` the sum of `gain^o` over the +octaves (signed, so the alternating sum for a negative gain), its value is `2 value +- G`, and its gradient, curvature and fold are twice `span_ridged`'s. Bounds come from the three pins only. With `|M| = |frequency| max(1, 1/stretch)` and `nu = f_o |M| speed` the noise-space speed, the Hessian of the noise on the span has -operator norm at most `Hl = min(SNOISE_HESSIAN, ||snoise_hessian(x_mid)||_F + +operator norm at most `Hl = min(SNOISE_HESSIAN, ||snoise_hessian(x_mid)||_F + e_h + SNOISE_THIRD nu halfLength)`, so `curvature = sum gain^o f_o^2 |M|^2 Hl` bounds the operator norm of the field's Hessian in `q` at every point of the span (the pin takes over once the third-derivative term exceeds it, and at `halfLength` 0 the local tier alone remains). For `span_ridged`, an octave leaves the branch it is on at the midpoint (`1 - s n`, `s` the midpoint's sign of `n`) only where `n` changes sign, and -there `|n| <= V - |n_mid|`, with `V = |grad n_mid . w| halfLength + (Hl nu |w| + Gl f_o -|M| acceleration) halfLength^2 / 2`, `w = f_o M velocity` and `Gl = min(SNOISE_LIPSCHITZ, -|grad n_mid| + Hl nu halfLength)`; `fold = sum 2 gain^o max(V - |n_mid|, 0)`, exactly 0 -for `span_fbm` and for an octave that cannot reach zero. Rounding is not widened: the -Frobenius norm exceeds the operator norm except on a rank-one Hessian, and the pins -carry a 10% margin. `SpanStageTests` encloses the Hessian norm and the branch -departure at 64 points of seeded spans, in the local and the global regime. +there `|n| <= V - |n_mid|`, with `V = (|grad n_mid . w| + e_g |w|) halfLength + (Hl nu |w| ++ Gl f_o |M| acceleration) halfLength^2 / 2 + 2 e_n`, `w = f_o M velocity` and `Gl = +min(SNOISE_LIPSCHITZ, |grad n_mid| + e_g + Hl nu halfLength)`; `fold = sum 2 gain^o +max(V - |n_mid|, 0)`, exactly 0 for `span_fbm` and for an octave that cannot reach zero. + +The float32 `snoise`, `snoise_grad` and `snoise_hessian` are not the exact functions, so +the quantities read at the midpoint carry allowances, by the derivation of `af_snoise`'s +comment. `S_o = |f_o| (|M q|_1 + |shift|_1) + |offset|_1 + |o octaveStep|_1 + sqrt(3) nu +halfLength` is the sum of the magnitudes the octave coordinate's roundings are relative +to, and `e_n = 2^-14 + L 2^-20 S_o` (value), `e_g = 2^-12 + H 2^-20 S_o` (gradient +length) and `e_h = 2^-10 + T 2^-20 S_o` (Frobenius norm of the Hessian), with `L`, `H`, +`T` the three pins; they add no pin. Measured against a double kernel with the float +hash (6000 points per band, `|x|` from 2^-2 to 2^23), the error is at most 0.057 of +`e_n`, 0.046 of `e_g` and 0.075 of `e_h` (the last at the smallest band, where the +absolute Hessian error is 2.3e-5 and the absolute term carries it). Errors grow +linearly in `|x|`: the value's is 2e-4 at 2^8, 1.5e-2 at 2^14 and 0.94 at 2^20. +Rounding in the Frobenius-to-operator-norm step needs no allowance: the Frobenius norm +exceeds the operator norm. + +Domain. No coordinate range is checked, clamped or refused. The allowances grow with +`S_o`, so `curvature` saturates at `SNOISE_HESSIAN` (true everywhere) by `S_o` of about +1.1e5, and `fold` turns positive for every octave whose `V + 2 e_n` reaches `|n_mid|`: +the bounds stay sound at every magnitude and stop being informative where float32 noise +is (see "Precision domain" above). `SpanStageTests` encloses the exact (double) Hessian +norm and both the float and the exact branch departure at 64 points of seeded spans in +bands of `S_o` from under 1e2 to over 1e6, and in a fixed span where the local tier +with its third-derivative term halved fails. 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/SpanPath.cs b/packages/cultmath/src/CultMath/SpanPath.cs new file mode 100644 index 00000000..cc6524d9 --- /dev/null +++ b/packages/cultmath/src/CultMath/SpanPath.cs @@ -0,0 +1,26 @@ +namespace CultMath; + +/// +/// The span and bound over: the path q(t), +/// |t| <= , through at t = 0 with there, +/// |q'| <= and |q''| <= along it. speed is at least +/// |velocity|; speed, acceleration and halfLength are non-negative. Name every field: the three floats +/// have no order a caller could mistake for a convention. +/// +public struct SpanPath +{ + /// q(0), the point the field and its gradient are evaluated at. + public float3 midpoint; + + /// q'(0). + public float3 velocity; + + /// An upper bound on |q'(t)| over the span; at least |velocity|. + public float speed; + + /// An upper bound on |q''(t)| over the span. + public float acceleration; + + /// The span covers |t| <= halfLength; 0 is the midpoint alone. + public float halfLength; +} diff --git a/packages/cultmath/src/CultMath/SpanPath.cs.meta b/packages/cultmath/src/CultMath/SpanPath.cs.meta new file mode 100644 index 00000000..3c58b0e9 --- /dev/null +++ b/packages/cultmath/src/CultMath/SpanPath.cs.meta @@ -0,0 +1,2 @@ +fileFormatVersion: 2 +guid: 3a8c1d52e7b94f06a1d5c2b9e04f7a31 diff --git a/packages/cultmath/src/CultMath/SpanStage.cs b/packages/cultmath/src/CultMath/SpanStage.cs new file mode 100644 index 00000000..faeccd3e --- /dev/null +++ b/packages/cultmath/src/CultMath/SpanStage.cs @@ -0,0 +1,33 @@ +namespace CultMath; + +/// +/// The transform of one octave stage of a warp chain, as and +/// read it. Octave o samples snoise at x_o = f_o (M q + shift) + offset + +/// o octaveStep, with f_o = lacunarity^o and M = (I + (1/ - 1) a a^T), +/// a = normalize(); a zero axis, or a stretch of 1, means no stretch. stretch must be +/// positive and finite, and nothing guards it. default(SpanStage) has frequency 0, which samples +/// one point: start from . +/// +public struct SpanStage +{ + /// The stage whose span_fbm/span_ridged field equals fbm_grad/ridged_grad: frequency 1, stretch 1, every vector zero. + public static readonly SpanStage Identity = new SpanStage { frequency = 1.0f, stretch = 1.0f }; + + /// The scale of M; its sign is kept. + public float frequency; + + /// The factor the axis direction is stretched by (M scales it by 1/stretch); 1 is none. + public float stretch; + + /// The direction stretch applies along; any length, zero for none. + public float3 axis; + + /// Added to M q before the per-octave scaling. + public float3 shift; + + /// Added after the per-octave scaling. + public float3 offset; + + /// Added to the octave coordinate once per octave index: o octaveStep. + public float3 octaveStep; +} diff --git a/packages/cultmath/src/CultMath/SpanStage.cs.meta b/packages/cultmath/src/CultMath/SpanStage.cs.meta new file mode 100644 index 00000000..7edaee93 --- /dev/null +++ b/packages/cultmath/src/CultMath/SpanStage.cs.meta @@ -0,0 +1,2 @@ +fileFormatVersion: 2 +guid: 9d46b0f27c1e4a85b3e8f1a60c2d5b97 diff --git a/packages/cultmath/src/CultMath/math.Span.cs b/packages/cultmath/src/CultMath/math.Span.cs index 4df9cd21..73bf199e 100644 --- a/packages/cultmath/src/CultMath/math.Span.cs +++ b/packages/cultmath/src/CultMath/math.Span.cs @@ -135,55 +135,78 @@ public static partial class math // Plain fbm of snoise over a span: the octave stage of a warp chain, its value and gradient at the // span's midpoint and the bounds that hold over the whole span. See span_stage for the span, the - // transform and the bound rules. At frequency 1, stretch 1 and zero shift, offset and octaveStep the - // field equals fbm_grad(q, octaves, lacunarity, gain) bit for bit. - public static SpanNoise span_fbm(float3 q, float3 velocity, float speed, float acceleration, float halfLength, - float frequency, float stretch, float3 axis, float3 shift, float3 offset, float3 octaveStep, int octaves, float lacunarity, float gain) => - span_stage(false, q, velocity, speed, acceleration, halfLength, frequency, stretch, axis, shift, offset, octaveStep, octaves, lacunarity, gain); + // transform and the bound rules. At SpanStage.Identity the field equals fbm_grad(path.midpoint, + // octaves, lacunarity, gain) bit for bit. + public static SpanNoise span_fbm(SpanPath path, SpanStage stage, int octaves, float lacunarity, float gain) => + span_stage(false, path, stage, octaves, lacunarity, gain); - // Ridged fbm over a span: octave o contributes gain^o (1 - |n_o|), ridged_grad's convention (the Span - // Lab's 1 - 2|n| is twice this minus the gain sum, the consumer's to form). fold bounds how far the - // sum can leave the branch it is on at the midpoint. At the identity transform the field equals - // ridged_grad bit for bit. - public static SpanNoise span_ridged(float3 q, float3 velocity, float speed, float acceleration, float halfLength, - float frequency, float stretch, float3 axis, float3 shift, float3 offset, float3 octaveStep, int octaves, float lacunarity, float gain) => - span_stage(true, q, velocity, speed, acceleration, halfLength, frequency, stretch, axis, shift, offset, octaveStep, octaves, lacunarity, gain); + // Ridged fbm over a span: octave o contributes gain^o (1 - |n_o|), ridged_grad's convention. fold + // bounds how far the sum can leave the branch it is on at the midpoint. At SpanStage.Identity the field + // equals ridged_grad bit for bit. + // + // The Span Lab's stage is 1 - 2|n| per octave. With G the sum of gain^o over the octaves (signed: a + // negative gain gives its alternating sum), the lab's value is 2 value - G, and its gradient, + // curvature and fold are twice this one's. A consumer in the lab form applies that mapping and does + // not derive its own bounds. + public static SpanNoise span_ridged(SpanPath path, SpanStage stage, int octaves, float lacunarity, float gain) => + span_stage(true, path, stage, octaves, lacunarity, gain); - // The span is the path q(t), |t| <= halfLength, through q at t = 0 with the given velocity there, - // |q'| <= speed and |q''| <= acceleration along it (speed at least |velocity|; all three and halfLength - // non-negative). Octave o samples snoise_grad at x_o(q) = f_o (M q + shift) + offset + o octaveStep with - // f_o = lacunarity^o, amplitude gain^o and M = frequency (I + (1/stretch - 1) a a^T), a = normalize(axis); - // a zero axis, or stretch 1, means no stretch. stretch must be positive and finite, and nothing guards - // it, as nothing guards lacunarity. octaves clamps to [0, 16] as fbm_grad's does. The field is the - // octave sum in q: gradient f_o M grad n_o (M is symmetric), the same octave loop as fbm_grad and - // ridged_grad. + private static float span_l1(float3 v) => abs(v.x) + abs(v.y) + abs(v.z); + + // The span is SpanPath: q(t), |t| <= halfLength, through midpoint with velocity at t = 0, |q'| <= speed + // and |q''| <= acceleration along it (speed at least |velocity|; all three and halfLength non-negative). + // The transform is SpanStage: octave o samples snoise_grad at x_o(q) = f_o (M q + shift) + offset + + // o octaveStep with f_o = lacunarity^o, amplitude gain^o and M = frequency (I + (1/stretch - 1) a a^T), + // a = normalize(axis); a zero axis, or stretch 1, means no stretch. stretch must be positive and + // finite, and nothing guards it, as nothing guards lacunarity. octaves clamps to [0, 16] as fbm_grad's + // does. The field is the octave sum in q: gradient f_o M grad n_o (M is symmetric), the same octave + // loop as fbm_grad and ridged_grad. // // Bounds, per octave, from SNOISE_LIPSCHITZ (L), SNOISE_HESSIAN (H) and SNOISE_THIRD (T) only. With // |M| = |frequency| max(1, 1/stretch) (1 without stretch) and nu = f_o |M| speed, the noise-space speed: // * |x(t) - x_mid| <= nu halfLength, so the Hessian of n at x(t) has operator norm at most - // Hl = min(H, ||hessian(x_mid)||_F + T nu halfLength) (the Hessian is T-Lipschitz in operator norm, - // and the Frobenius norm bounds the operator norm), and |grad n(x(t))| at most - // Gl = min(L, |grad n_mid| + Hl nu halfLength). + // Hl = min(H, ||hessian(x_mid)||_F + e_h + T nu halfLength) (the Hessian is T-Lipschitz in operator + // norm, and the Frobenius norm bounds the operator norm), and |grad n(x(t))| at most + // Gl = min(L, |grad n_mid| + e_g + Hl nu halfLength). // * curvature = sum gain^o f_o^2 |M|^2 Hl: the Hessian in q of n(x_o(q)) is f_o^2 M H M. // * fold, ridged only. Over the span the octave's value departs from its midpoint branch // gain^o (1 - s n) (s the sign of n at the midpoint) only where n has changed sign, and there // |n(t)| <= V - |n_mid|, V bounding |n(t) - n_mid|. Writing n(t) - n_mid as the integral of // grad n . x', x' = w + rho with w = f_o M velocity and |rho| <= f_o |M| acceleration |t|, and // grad n = grad n_mid + delta with |delta| <= Hl nu |t|: - // V = |grad n_mid . w| halfLength + (Hl nu |w| + Gl f_o |M| acceleration) halfLength^2 / 2, + // V = (|grad n_mid . w| + e_g |w|) halfLength + (Hl nu |w| + Gl f_o |M| acceleration) halfLength^2 / 2 + // + 2 e_n, // and the octave's fold is 2 gain^o max(V - |n_mid|, 0), zero when it cannot reach zero. - // Rounding is not widened: the Frobenius norm exceeds the operator norm except on a rank-one - // Hessian, and the pins carry a 10% margin. - private static SpanNoise span_stage(bool ridged, float3 q, float3 velocity, float speed, float acceleration, float halfLength, - float frequency, float stretch, float3 axis, float3 shift, float3 offset, float3 octaveStep, int octaves, float lacunarity, float gain) + // + // The float32 snoise, snoise_grad and snoise_hessian are not the exact functions, and the field a + // consumer reads is the float one, so what is read at the midpoint carries an allowance, by the + // derivation af_snoise's comment gives (math.Affine.cs): the skew and the corner offsets round the + // coordinate by about 2u of the magnitudes it is formed from, which moves a result by its derivative + // times that, and the caller's scaling of the point and a GPU's evaluation take 2u to 16u = 2^-20. + // S_o is the sum of those magnitudes, |f_o| (|M q|_1 + |shift|_1) + |offset|_1 + |o octaveStep|_1, plus + // the span's reach sqrt(3) nu halfLength (the coordinate of a path point). The allowances are the + // value's e_n = 2^-14 + L 2^-20 S_o, the gradient length's e_g = 2^-12 + H 2^-20 S_o (H is the + // gradient's own derivative's pin) and the Hessian's Frobenius norm's e_h = 2^-10 + T 2^-20 S_o (T is + // the Hessian's derivative's pin): af_snoise's absolute terms for the first two, and 2^-10 for the + // Hessian's, which the kernel's rounding at small coordinates (2.3e-5 absolute, measured against a + // double kernel) sits well under. The same probe (6000 points per band, |x| from 2^-2 to 2^23) puts + // the error at most 0.057, 0.046 and 0.075 of e_n, e_g and e_h. No new pin enters. + // + // Domain: no coordinate range is checked, clamped or refused. The allowances grow with S_o, so + // curvature saturates at H (true everywhere) by S_o of about 1.1e5, and fold turns positive for every + // octave whose V + 2 e_n reaches |n_mid|: the bounds stay sound at every magnitude and stop being + // informative where float32 noise is (design.md, Precision domain). + private static SpanNoise span_stage(bool ridged, SpanPath path, SpanStage stage, int octaves, float lacunarity, float gain) { + var q = path.midpoint; + var halfLength = path.halfLength; octaves = clamp(octaves, 0, 16); - var stretched = dot(axis, axis) > 0.0f && stretch != 1.0f; - var direction = stretched ? normalize(axis) : new float3(0.0f, 0.0f, 0.0f); - var squeeze = stretched ? 1.0f / stretch - 1.0f : 0.0f; - var norm = abs(frequency) * (stretched ? max(1.0f, 1.0f / stretch) : 1.0f); - var mq = stretched ? frequency * (q + squeeze * dot(direction, q) * direction) : frequency * q; - var mv = stretched ? frequency * (velocity + squeeze * dot(direction, velocity) * direction) : frequency * velocity; + var stretched = dot(stage.axis, stage.axis) > 0.0f && stage.stretch != 1.0f; + var direction = stretched ? normalize(stage.axis) : new float3(0.0f, 0.0f, 0.0f); + var squeeze = stretched ? 1.0f / stage.stretch - 1.0f : 0.0f; + var norm = abs(stage.frequency) * (stretched ? max(1.0f, 1.0f / stage.stretch) : 1.0f); + var mq = stretched ? stage.frequency * (q + squeeze * dot(direction, q) * direction) : stage.frequency * q; + var mv = stretched ? stage.frequency * (path.velocity + squeeze * dot(direction, path.velocity) * direction) : stage.frequency * path.velocity; var amplitude = 1.0f; var scale = 1.0f; var value = 0.0f; @@ -192,10 +215,11 @@ public static partial class math var fold = 0.0f; for (var i = 0; i < octaves; i++) { - var x = scale * (mq + shift) + (offset + (float)i * octaveStep); + var step = (float)i * stage.octaveStep; + var x = scale * (mq + stage.shift) + (stage.offset + step); var n = snoise_grad(x); var g = new float3(n.x, n.y, n.z); - var mg = stretched ? frequency * (g + squeeze * dot(direction, g) * direction) : frequency * g; + var mg = stretched ? stage.frequency * (g + squeeze * dot(direction, g) * direction) : stage.frequency * g; if (ridged) { var s = (float)sign(n.w); @@ -210,15 +234,19 @@ public static partial class math var weight = abs(amplitude); var factor = abs(scale) * norm; - var nu = factor * speed; - var local = min(SNOISE_HESSIAN, span_frobenius(snoise_hessian(x)) + SNOISE_THIRD * nu * halfLength); + var nu = factor * path.speed; + var rounding = 9.5367431640625e-7f * (abs(scale) * (span_l1(mq) + span_l1(stage.shift)) + span_l1(stage.offset) + span_l1(step) + 1.7320508f * nu * halfLength); + var valueSlack = 6.103515625e-5f + SNOISE_LIPSCHITZ * rounding; + var gradientSlack = 2.44140625e-4f + SNOISE_HESSIAN * rounding; + var hessianSlack = 9.765625e-4f + SNOISE_THIRD * rounding; + var local = min(SNOISE_HESSIAN, span_frobenius(snoise_hessian(x)) + hessianSlack + SNOISE_THIRD * nu * halfLength); curvature += weight * factor * factor * local; if (ridged) { var w = scale * mv; var reach = local * nu * halfLength; - var gl = min(SNOISE_LIPSCHITZ, length(g) + reach); - var travel = abs(dot(g, w)) * halfLength + (local * nu * length(w) + gl * factor * acceleration) * (0.5f * halfLength * halfLength); + var gl = min(SNOISE_LIPSCHITZ, length(g) + gradientSlack + reach); + var travel = (abs(dot(g, w)) + gradientSlack * length(w)) * halfLength + (local * nu * length(w) + gl * factor * path.acceleration) * (0.5f * halfLength * halfLength) + 2.0f * valueSlack; fold += 2.0f * weight * max(travel - abs(n.w), 0.0f); } diff --git a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs index 8d61150a..0962a357 100644 --- a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs @@ -633,7 +633,7 @@ public sealed class NoiseBoundTests // gradient direction p_c, in double, with snoise_hessian's lattice, step() tie order and (float) hash, // so one corner table serves the Hessian and its derivative. Key names the lattice cell and traversal // order: two points with the same key are differentiated through the same corners. - private static (double[][] X, double[][] P, double[] Key) Corners(double vx, double vy, double vz) + internal static (double[][] X, double[][] P, double[] Key) Corners(double vx, double vy, double vz) { const double cx = 1.0 / 6.0; const double cy = 1.0 / 3.0; @@ -686,7 +686,7 @@ public sealed class NoiseBoundTests // The Hessian of snoise at v in double, entry i * 3 + j: snoise_hessian's per-corner derivation, summed // over Corners. A corner with m = 0.5 - |x|^2 <= 0 contributes nothing. - private static double[] HessianDouble(double vx, double vy, double vz) + internal static double[] HessianDouble(double vx, double vy, double vz) { var (x, p, _) = Corners(vx, vy, vz); var hessian = new double[9]; @@ -705,6 +705,23 @@ public sealed class NoiseBoundTests return hessian; } + // The value of snoise at v in double: 105 sum m^4 s over Corners, the kernel's exact value (the tied + // kernel, as HessianDouble is its exact Hessian). + internal static double ValueDouble(double vx, double vy, double vz) + { + var (x, p, _) = Corners(vx, vy, vz); + var total = 0.0; + 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; + total += 105.0 * m * m * m * m * (p[c][0] * x[c][0] + p[c][1] * x[c][1] + p[c][2] * x[c][2]); + } + + return total; + } + 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 bc6c6eeb..8bd68193 100644 --- a/packages/cultmath/tests/CultMath.Tests/SpanStageTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/SpanStageTests.cs @@ -32,48 +32,76 @@ public sealed class SpanStageTests return normalize(v); } - // One draw of the stage parameters: every one of them off its identity value. - private readonly record struct Stage(float Frequency, float Stretch, float3 Axis, float3 Shift, float3 Offset, float3 Step, int Octaves, float Lacunarity, float Gain) + private static float LogUniform(System.Random random, float lo, float hi) => MathF.Exp(Uniform(random, MathF.Log(lo), MathF.Log(hi))); + + private static float Signed(System.Random random, float magnitude) => random.Next(2) == 0 ? magnitude : -magnitude; + + private static readonly float3 None = new(0.0f, 0.0f, 0.0f); + + private static double L1(double[] v) => Math.Abs(v[0]) + Math.Abs(v[1]) + Math.Abs(v[2]); + + // One draw of the stage parameters: every one of them off its identity value, over the ranges the stage + // takes in use (stage-parameter-range): |frequency| 1e-3..10, stretch 1e-2..1e2, |lacunarity| 0.5..3, + // |gain| 0.1..1.2, the last three and frequency with random signs where a sign is meaningful, + // octaveStep nonzero in every component. + private readonly record struct Stage(SpanStage S, int Octaves, float Lacunarity, float Gain) { - public static Stage Draw(System.Random random, int octaves) => new( - Uniform(random, 0.3f, 1.6f), - random.Next(2) == 0 ? Uniform(random, 0.2f, 0.8f) : Uniform(random, 1.25f, 5.0f), - UnitVector(random) * Uniform(random, 0.5f, 3.0f), - Vector(random, -8.0f, 8.0f), - Vector(random, -8.0f, 8.0f), - Vector(random, -3.0f, 3.0f), - octaves, - Uniform(random, 1.5f, 2.3f), - Uniform(random, 0.35f, 0.7f)); + public static Stage Draw(System.Random random, int maxOctaves = 8, int minOctaves = 1) => new( + new SpanStage + { + frequency = Signed(random, LogUniform(random, 1.0e-3f, 10.0f)), + stretch = LogUniform(random, 1.0e-2f, 1.0e2f), + axis = UnitVector(random) * Uniform(random, 0.5f, 3.0f), + shift = Vector(random, -8.0f, 8.0f), + offset = Vector(random, -8.0f, 8.0f), + octaveStep = new float3(Signed(random, Uniform(random, 0.1f, 3.0f)), Signed(random, Uniform(random, 0.1f, 3.0f)), Signed(random, Uniform(random, 0.1f, 3.0f))), + }, + random.Next(minOctaves, maxOctaves + 1), + Signed(random, Uniform(random, 0.5f, 3.0f)), + Signed(random, Uniform(random, 0.1f, 1.2f))); - public SpanNoise Fbm(float3 q, float3 velocity, float speed, float acceleration, float halfLength) => - span_fbm(q, velocity, speed, acceleration, halfLength, Frequency, Stretch, Axis, Shift, Offset, Step, Octaves, Lacunarity, Gain); + // A draw whose octave coordinates stay under `limit` for every q with |q|_1 <= qL1: the fastest + // octave's |f_o| |M| is capped by rejection. For the derivative checks, where float rounding of the + // coordinate must not swamp a difference quotient. + public static Stage DrawSmooth(System.Random random, float limit, float qL1, int maxOctaves = 8, int minOctaves = 1) + { + while (true) + { + var stage = Draw(random, maxOctaves, minOctaves); + var worst = 0.0; + for (var o = 0; o < stage.Octaves; o++) + worst = Math.Max(worst, Math.Abs(stage.Scale(o)) * (stage.Norm() * qL1 + L1(Components(stage.S.shift))) + L1(Components(stage.S.offset)) + o * L1(Components(stage.S.octaveStep))); + if (worst < limit) + return stage; + } + } - public SpanNoise Ridged(float3 q, float3 velocity, float speed, float acceleration, float halfLength) => - span_ridged(q, velocity, speed, acceleration, halfLength, Frequency, Stretch, Axis, Shift, Offset, Step, Octaves, Lacunarity, Gain); + public SpanNoise Fbm(SpanPath path) => span_fbm(path, S, Octaves, Lacunarity, Gain); + + public SpanNoise Ridged(SpanPath path) => span_ridged(path, S, Octaves, Lacunarity, Gain); // M = frequency (I + (1/stretch - 1) a a^T) in double, a the normalised axis (none for a zero axis). public double[,] Matrix() { var m = new double[3, 3]; - var a = new double[] { Axis.x, Axis.y, Axis.z }; + var a = Components(S.axis); var length = Math.Sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2]); for (var i = 0; i < 3; i++) for (var j = 0; j < 3; j++) - m[i, j] = Frequency * ((i == j ? 1.0 : 0.0) + (length > 0.0 ? (1.0 / Stretch - 1.0) * a[i] * a[j] / (length * length) : 0.0)); + m[i, j] = S.frequency * ((i == j ? 1.0 : 0.0) + (length > 0.0 && S.stretch != 1.0f ? (1.0 / S.stretch - 1.0) * a[i] * a[j] / (length * length) : 0.0)); return m; } - // Octave o's noise-space point for the span point q. - public float3 Point(float3 q, int octave) + // Octave o's noise-space point for the span point q, in double from the same inputs. + public double[] PointD(float3 q, int octave) { var m = Matrix(); var scale = Math.Pow(Lacunarity, octave); var x = new double[3]; - var qv = new double[] { q.x, q.y, q.z }; - var shift = new double[] { Shift.x, Shift.y, Shift.z }; - var offset = new double[] { Offset.x, Offset.y, Offset.z }; - var step = new double[] { Step.x, Step.y, Step.z }; + var qv = Components(q); + var shift = Components(S.shift); + var offset = Components(S.offset); + var step = Components(S.octaveStep); for (var i = 0; i < 3; i++) { var mq = 0.0; @@ -82,28 +110,59 @@ public sealed class SpanStageTests x[i] = scale * (mq + shift[i]) + offset[i] + octave * step[i]; } + return x; + } + + public float3 Point(float3 q, int octave) + { + var x = PointD(q, octave); return new float3((float)x[0], (float)x[1], (float)x[2]); } + // S_o: the magnitudes octave o's coordinate roundings are relative to, plus the span's reach. + public double Magnitude(SpanPath path, int octave) + { + var m = Matrix(); + var qv = Components(path.midpoint); + var mq = new double[3]; + for (var i = 0; i < 3; i++) + for (var j = 0; j < 3; j++) + mq[i] += m[i, j] * qv[j]; + var scale = Math.Abs(Scale(octave)); + return scale * (L1(mq) + L1(Components(S.shift))) + L1(Components(S.offset)) + octave * L1(Components(S.octaveStep)) + + Math.Sqrt(3.0) * scale * Norm() * path.speed * path.halfLength; + } + + public double MaxMagnitude(SpanPath path) => Enumerable.Range(0, Octaves).Max(o => Magnitude(path, o)); + + // Octave o as a stage of its own: one octave at the scaled frequency, shift and offset. + public Stage OneOctave(int octave) => new( + new SpanStage { frequency = S.frequency * (float)Scale(octave), stretch = S.stretch, axis = S.axis, shift = S.shift * (float)Scale(octave), offset = S.offset + (float)octave * S.octaveStep, octaveStep = None }, + 1, 1.0f, 1.0f); + public double Amplitude(int octave) => Math.Pow(Gain, octave); public double Scale(int octave) => Math.Pow(Lacunarity, octave); // The spectral norm of M: frequency times the larger of 1 and 1/stretch. - public double Norm() => Math.Abs(Frequency) * Math.Max(1.0, 1.0 / Stretch); + public double Norm() => Math.Abs(S.frequency) * Math.Max(1.0, 1.0 / S.stretch); } + private static double[] Components(float3 v) => new double[] { v.x, v.y, v.z }; + + private static SpanPath PathAt(float3 q) => new() { midpoint = q, velocity = None, speed = 0.0f, acceleration = 0.0f, halfLength = 0.0f }; + // A span through q: velocity v, constant acceleration c of norm `acceleration`, so q(t) = q + v t + c t^2 / 2; // the speed bound is |v| plus what c can add over halfLength. - private readonly record struct Span(float3 Q, float3 Velocity, float3 Acceleration, float Speed, float AccelerationBound, float HalfLength) + private readonly record struct Span(SpanPath Path, float3 Acceleration) { public static Span Draw(System.Random random, float3 q, float halfLength, float acceleration) { var velocity = UnitVector(random) * Uniform(random, 0.2f, 3.0f); - return new Span(q, velocity, UnitVector(random) * acceleration, length(velocity) + acceleration * halfLength, acceleration, halfLength); + return new Span(new SpanPath { midpoint = q, velocity = velocity, speed = length(velocity) + acceleration * halfLength, acceleration = acceleration, halfLength = halfLength }, UnitVector(random) * acceleration); } - public float3 At(float t) => Q + Velocity * t + Acceleration * (0.5f * t * t); + public float3 At(float t) => Path.midpoint + Path.velocity * t + Acceleration * (0.5f * t * t); } // The operator norm of a symmetric 3x3 matrix (entries in row-major order): the larger |eigenvalue|. @@ -121,28 +180,25 @@ public sealed class SpanStageTests return Math.Max(Math.Abs(q + 2.0 * s * Math.Cos(phi)), Math.Abs(q + 2.0 * s * Math.Cos(phi + 2.0 * Math.PI / 3.0))); } - // The Hessian in q of the stage field at span point q (signed per octave for a ridged sum), by the chain - // rule from snoise_hessian: sum over octaves of weight_o f_o^2 M H_o M. + // The exact Hessian in q of the stage field at span point q (signed per octave for a ridged sum), by the + // chain rule from the double kernel's Hessian: sum over octaves of weight_o f_o^2 M H_o M, with every + // octave point computed in double from the same SpanStage. private static double[] FieldHessian(Stage stage, float3 q, bool ridged) { var m = stage.Matrix(); var total = new double[9]; for (var o = 0; o < stage.Octaves; o++) { - var x = stage.Point(q, o); - var h = snoise_hessian(x); - var weight = stage.Amplitude(o) * stage.Scale(o) * stage.Scale(o) * (ridged ? -Math.Sign(snoise(x)) : 1.0); - var hd = new double[3, 3]; - for (var i = 0; i < 3; i++) - for (var j = 0; j < 3; j++) - hd[i, j] = h[i][j]; + var x = stage.PointD(q, o); + var h = NoiseBoundTests.HessianDouble(x[0], x[1], x[2]); + var weight = stage.Amplitude(o) * stage.Scale(o) * stage.Scale(o) * (ridged ? -Math.Sign(NoiseBoundTests.ValueDouble(x[0], x[1], x[2])) : 1.0); for (var i = 0; i < 3; i++) for (var j = 0; j < 3; j++) { var sum = 0.0; for (var k = 0; k < 3; k++) for (var l = 0; l < 3; l++) - sum += m[i, k] * hd[k, l] * m[l, j]; + sum += m[i, k] * h[k * 3 + l] * m[l, j]; total[i * 3 + j] += weight * sum; } } @@ -168,24 +224,41 @@ public sealed class SpanStageTests } /// - /// At frequency 1, stretch 1 (with and without an axis) and zero shift, offset and octaveStep, the field of + /// The exact reference the enclosure tests lean on is the kernel: the double ValueDouble is snoise to + /// within float rounding at small coordinates, so a sign taken from it is the kernel's sign. + /// + [Fact] + public void ExactReferenceIsTheKernel() + { + var random = new System.Random(0xE4AC7); + for (var k = 0; k < 300; k++) + { + var p = Vector(random, -20.0f, 20.0f); + Assert.InRange(NoiseBoundTests.ValueDouble(p.x, p.y, p.z) - snoise(p), -2.0e-4, 2.0e-4); + } + } + + /// + /// At SpanStage.Identity (also with an axis, stretch 1) and zero shift, offset and octaveStep, the field of /// span_fbm and span_ridged is fbm_grad's and ridged_grad's, bit for bit, for 0 to 20 octaves (the clamp /// included), every lacunarity and gain drawn, over random and simplex-tie points: any change to the octave - /// loop's order, scaling or accumulation fails. + /// loop's order, scaling or accumulation, and any bit the rounding allowances put in the field, fails. /// [Fact] public void SpanStagesAtIdentityAreFbmAndRidged() { var random = new System.Random(0x57A6E1); - var none = new float3(0.0f, 0.0f, 0.0f); foreach (var p in Points(random)) { var octaves = random.Next(0, 21); var lacunarity = Uniform(random, 1.4f, 2.6f); var gain = Uniform(random, 0.3f, 0.8f); - var axis = random.Next(2) == 0 ? none : UnitVector(random) * 2.0f; - var fbm = span_fbm(p, UnitVector(random), 1.0f, 0.5f, 0.25f, 1.0f, 1.0f, axis, none, none, none, octaves, lacunarity, gain); - var ridged = span_ridged(p, UnitVector(random), 1.0f, 0.5f, 0.25f, 1.0f, 1.0f, axis, none, none, none, octaves, lacunarity, gain); + var stage = SpanStage.Identity; + if (random.Next(2) == 1) + stage.axis = UnitVector(random) * 2.0f; + var path = new SpanPath { midpoint = p, velocity = UnitVector(random), speed = 1.0f, acceleration = 0.5f, halfLength = 0.25f }; + var fbm = span_fbm(path, stage, octaves, lacunarity, gain); + var ridged = span_ridged(path, stage, octaves, lacunarity, gain); Assert.True(SameBits(fbm.field, fbm_grad(p, octaves, lacunarity, gain)), $"span_fbm at {p}, {octaves} octaves: {fbm.field} is not fbm_grad {fbm_grad(p, octaves, lacunarity, gain)}"); Assert.True(SameBits(ridged.field, ridged_grad(p, octaves, lacunarity, gain)), $"span_ridged at {p}, {octaves} octaves: {ridged.field} is not ridged_grad {ridged_grad(p, octaves, lacunarity, gain)}"); Assert.Equal(0.0f, fbm.fold); @@ -194,10 +267,12 @@ public sealed class SpanStageTests /// /// The gradient is the derivative of the value through the whole transform: with every parameter off its - /// identity (frequency, stretch both ways, axis, shift, offset, octaveStep, lacunarity, gain), field.xyz - /// equals the fourth-order central difference of field.w along each axis, within a tolerance built from the - /// third-derivative pin (truncation) and float32 rounding of the coordinates. A dropped M or lacunarity - /// factor, or M applied to the wrong side, moves it by the gradient itself. + /// identity (frequency of either sign, stretch both ways, axis, shift, offset, octaveStep, lacunarity and + /// gain of either sign, 4 to 8 octaves) field.xyz equals the fourth-order central difference of field.w + /// along each axis, within a tolerance built from the third-derivative pin (truncation) and float32 + /// rounding of the coordinates, and both agree with an independent octave sum. A dropped M or lacunarity + /// factor, M applied to the wrong side, a clamp on frequency, stretch or the octaveStep index, or a dropped + /// sign (abs of the scale or the amplitude) moves the field by the gradient itself. /// [Theory] [InlineData(false)] @@ -207,15 +282,23 @@ public sealed class SpanStageTests var random = new System.Random(ridged ? 0x6A1D : 0x6A1C); var worst = 0.0; var checkedCases = 0; - for (var k = 0; k < 400; k++) + var spans = 0; + for (var k = 0; k < 6000 && spans < 80; k++) { - var stage = Stage.Draw(random, ridged ? random.Next(1, 3) : random.Next(1, 4)); + var stage = Stage.DrawSmooth(random, 300.0f, 9.0f, 8, 4); var q = Vector(random, -3.0f, 3.0f); // The fold's crease at n = 0 is not differentiable: keep every octave's stencil clear of it. - if (ridged && Enumerable.Range(0, stage.Octaves).Any(o => abs(snoise(stage.Point(q, o))) < 0.3f)) + if (ridged && Enumerable.Range(0, stage.Octaves).Any(o => abs(snoise(stage.Point(q, o))) < 0.15f)) continue; + spans++; var span = Span.Draw(random, q, 0.1f, 0.1f); - SpanNoise At(float3 p) => ridged ? stage.Ridged(p, span.Velocity, span.Speed, span.AccelerationBound, span.HalfLength) : stage.Fbm(p, span.Velocity, span.Speed, span.AccelerationBound, span.HalfLength); + SpanNoise At(float3 p) + { + var path = span.Path; + path.midpoint = p; + return ridged ? stage.Ridged(path) : stage.Fbm(path); + } + var centre = At(q); // The value and gradient from an independent assembly of the octaves: pins shift, offset, octaveStep, // M, lacunarity and gain themselves, which a value/gradient self-consistency check cannot see. @@ -233,18 +316,16 @@ public sealed class SpanStageTests referenceGradient[i] += sign * stage.Amplitude(o) * stage.Scale(o) * matrix[i, j] * g[j]; } - var weight = 0.0; - for (var o = 0; o < stage.Octaves; o++) - weight += stage.Amplitude(o) * stage.Scale(o) * stage.Norm(); - Assert.True(Math.Abs(referenceValue - centre.field.w) <= 5.0e-4 * (1.0 + weight), $"stage {stage}: the value is {centre.field.w:R}, an independent octave sum gives {referenceValue:R}"); - for (var i = 0; i < 3; i++) - Assert.True(Math.Abs(referenceGradient[i] - centre.field[i]) <= 5.0e-3 * (1.0 + weight), $"stage {stage}: gradient {"xyz"[i]} is {centre.field[i]:R}, an independent octave sum gives {referenceGradient[i]:R}"); - // A step of 0.02 in the noise space of the fastest octave. - var fastest = stage.Norm() * Math.Pow(stage.Lacunarity, stage.Octaves - 1); - var h = (float)(0.02 / fastest); var scale = 0.0; for (var o = 0; o < stage.Octaves; o++) - scale += stage.Amplitude(o) * stage.Scale(o) * stage.Norm(); + scale += Math.Abs(stage.Amplitude(o) * stage.Scale(o)) * stage.Norm(); + var amplitudes = Enumerable.Range(0, stage.Octaves).Sum(o => Math.Abs(stage.Amplitude(o))); + Assert.True(Math.Abs(referenceValue - centre.field.w) <= 2.0e-3 * (1.0 + amplitudes), $"stage {stage}: the value is {centre.field.w:R}, an independent octave sum gives {referenceValue:R}"); + for (var i = 0; i < 3; i++) + Assert.True(Math.Abs(referenceGradient[i] - centre.field[i]) <= 1.0e-2 * (1.0 + scale), $"stage {stage}: gradient {"xyz"[i]} is {centre.field[i]:R}, an independent octave sum gives {referenceGradient[i]:R}"); + // A step of 0.02 in the noise space of the fastest octave. + var fastest = stage.Norm() * Enumerable.Range(0, stage.Octaves).Max(o => Math.Abs(stage.Scale(o))); + var h = (float)(0.02 / fastest); for (var axis = 0; axis < 3; axis++) { var e = new float3(axis == 0 ? h : 0.0f, axis == 1 ? h : 0.0f, axis == 2 ? h : 0.0f); @@ -276,33 +357,46 @@ public sealed class SpanStageTests var random = new System.Random(0x57E7C4); for (var k = 0; k < 60; k++) { - var stage = Stage.Draw(random, 3); + var stage = Stage.DrawSmooth(random, 300.0f, 15.0f, 3, 3); var q = Vector(random, -5.0f, 5.0f); var span = Span.Draw(random, q, 0.2f, 0.3f); - var plain = stage with { Stretch = 1.0f }; - var noAxis = stage with { Axis = new float3(0.0f, 0.0f, 0.0f), Stretch = 3.0f }; - var reference = plain.Ridged(q, span.Velocity, span.Speed, span.AccelerationBound, span.HalfLength); - var isotropic = (plain with { Axis = new float3(0.0f, 0.0f, 0.0f) }).Ridged(q, span.Velocity, span.Speed, span.AccelerationBound, span.HalfLength); - var viaNoAxis = noAxis.Ridged(q, span.Velocity, span.Speed, span.AccelerationBound, span.HalfLength); + var plain = stage with { S = stage.S with { stretch = 1.0f } }; + var noAxis = stage with { S = stage.S with { axis = None, stretch = 3.0f } }; + var reference = plain.Ridged(span.Path); + var isotropic = (plain with { S = plain.S with { axis = None } }).Ridged(span.Path); + var viaNoAxis = noAxis.Ridged(span.Path); Assert.True(SameBits(reference.field, isotropic.field) && reference.curvature == isotropic.curvature && reference.fold == isotropic.fold, "stretch 1 with an axis differs from no axis"); Assert.True(SameBits(isotropic.field, viaNoAxis.field) && isotropic.curvature == viaNoAxis.curvature && isotropic.fold == viaNoAxis.fold, "a zero axis with stretch 3 differs from the isotropic stage"); - var unit = stage with { Axis = new float3(0.0f, 1.0f, 0.0f) }; - var scaled = stage with { Axis = new float3(0.0f, 2.0f, 0.0f) }; - var a = unit.Fbm(q, span.Velocity, span.Speed, span.AccelerationBound, span.HalfLength); - var b = scaled.Fbm(q, span.Velocity, span.Speed, span.AccelerationBound, span.HalfLength); + var unit = stage with { S = stage.S with { axis = new float3(0.0f, 1.0f, 0.0f) } }; + var scaled = stage with { S = stage.S with { axis = new float3(0.0f, 2.0f, 0.0f) } }; + var a = unit.Fbm(span.Path); + var b = scaled.Fbm(span.Path); Assert.True(SameBits(a.field, b.field) && a.curvature == b.curvature, "a non-unit axis differs from its normalisation"); - var stretched = unit.Fbm(q, span.Velocity, span.Speed, span.AccelerationBound, span.HalfLength); - Assert.False(SameBits(stretched.field, (unit with { Stretch = 1.0f }).Fbm(q, span.Velocity, span.Speed, span.AccelerationBound, span.HalfLength).field), "a stretch of " + unit.Stretch + " did not change the field"); + var stretched = unit.Fbm(span.Path); + Assert.False(SameBits(stretched.field, (unit with { S = unit.S with { stretch = 1.0f } }).Fbm(span.Path).field), "a stretch of " + unit.S.stretch + " did not change the field"); } } - // The stage parameters and per-octave bookkeeping the curvature checks share. + // The coordinate bands the enclosure tests draw: a large offset (0 keeps the drawn one) lifts S_o through + // the decades; Bucket names the decade a span's largest S_o is in. + private static readonly float[] OffsetScales = { 0.0f, 40.0f, 400.0f, 4.0e3f, 4.0e4f, 4.0e5f, 4.0e6f }; + + private static int Bucket(double magnitude) => magnitude < 1.0e2 ? 0 : magnitude < 1.0e3 ? 1 : magnitude < 1.0e4 ? 2 : magnitude < 1.0e5 ? 3 : magnitude < 1.0e6 ? 4 : 5; + + private static Stage Banded(System.Random random, float offsetScale, int octaves) + { + var stage = Stage.Draw(random, octaves); + if (offsetScale > 0.0f) + stage = stage with { S = stage.S with { offset = Vector(random, -offsetScale, offsetScale) } }; + return stage; + } + private static void CheckCurvature(Stage stage, bool ridged, Span span, SpanNoise result, int samples, ref double tightest, ref int checkedPoints) { for (var s = 0; s < samples; s++) { - var t = span.HalfLength * (2.0f * s / (samples - 1) - 1.0f); + var t = span.Path.halfLength * (2.0f * s / (samples - 1) - 1.0f); var norm = SymmetricNorm(FieldHessian(stage, span.At(t), ridged)); Assert.True(norm <= result.curvature * (1.0 + 1.0e-3), $"stage {stage}, span {span}: the Hessian norm at t = {t:R} is {norm:R}, above curvature {result.curvature:R}"); tightest = Math.Min(tightest, result.curvature > 0.0 ? 1.0 - norm / result.curvature : 1.0); @@ -311,13 +405,15 @@ public sealed class SpanStageTests } /// - /// span-bounds-enclose, curvature: over seeded spans the operator norm of the field's q-Hessian (chain rule - /// from snoise_hessian, M and the octave scales) at 64 points per span never exceeds curvature, in the - /// local regime (short spans, where the third-derivative term is the bound) and the global one (long spans, - /// where SNOISE_HESSIAN is); at halfLength 0 curvature is the local tier, the Frobenius norm of the - /// midpoint Hessians scaled by amplitude, f_o^2 and |M|^2, to within float rounding; and a span through - /// MeasureHessian's witness stays enclosed in 200 directions, which the Frobenius norm at the midpoint alone - /// does not guarantee. + /// span-bounds-enclose, curvature: over seeded spans the operator norm of the exact field q-Hessian (chain + /// rule from the double kernel's Hessian, every octave point in double, M and the octave scales) at 64 + /// points per span never exceeds curvature, in the local regime (short spans, where the third-derivative + /// term is the bound) and the global one (long spans, where SNOISE_HESSIAN is), and in every coordinate + /// band from S_o under 1e2 to over 1e6, where the float32 kernel is no longer the exact one; at + /// halfLength 0 curvature is the widened local tier, the Frobenius norm of the midpoint Hessians plus + /// e_h, capped at the pin, scaled by amplitude, f_o^2 and |M|^2; and a span through MeasureHessian's + /// witness stays enclosed in 200 directions, which the Frobenius norm at the midpoint alone does not + /// guarantee. /// [Theory] [InlineData(false)] @@ -327,28 +423,34 @@ public sealed class SpanStageTests var random = new System.Random(ridged ? 0xC0DE2 : 0xC0DE1); var tightest = 1.0; var checkedPoints = 0; - foreach (var halfLength in new[] { 0.002f, 0.02f, 0.2f, 1.0f, 8.0f }) - for (var k = 0; k < 24; k++) - { - var stage = Stage.Draw(random, random.Next(1, 4)); - var span = Span.Draw(random, Vector(random, -30.0f, 30.0f), halfLength, Uniform(random, 0.0f, 2.0f)); - var result = ridged ? stage.Ridged(span.Q, span.Velocity, span.Speed, span.AccelerationBound, span.HalfLength) : stage.Fbm(span.Q, span.Velocity, span.Speed, span.AccelerationBound, span.HalfLength); - CheckCurvature(stage, ridged, span, result, 64, ref tightest, ref checkedPoints); - } + var buckets = new int[6]; + foreach (var offsetScale in OffsetScales) + foreach (var halfLength in new[] { 0.02f, 0.5f, 8.0f }) + for (var k = 0; k < 8; k++) + { + var stage = Banded(random, offsetScale, 3); + var span = Span.Draw(random, Vector(random, -30.0f, 30.0f), halfLength, Uniform(random, 0.0f, 2.0f)); + var result = ridged ? stage.Ridged(span.Path) : stage.Fbm(span.Path); + buckets[Bucket(stage.MaxMagnitude(span.Path))]++; + CheckCurvature(stage, ridged, span, result, 64, ref tightest, ref checkedPoints); + } - // halfLength 0: the local tier exactly. + Assert.True(buckets.All(count => count >= 3), $"the coordinate bands drawn are {string.Join(", ", buckets)} spans: a band is not reached"); + + // halfLength 0: the widened local tier. for (var k = 0; k < 40; k++) { - var stage = Stage.Draw(random, random.Next(1, 4)); + var stage = Stage.DrawSmooth(random, 100.0f, 90.0f, 3); var q = Vector(random, -30.0f, 30.0f); var span = Span.Draw(random, q, 0.0f, 1.0f); - var result = ridged ? stage.Ridged(q, span.Velocity, span.Speed, span.AccelerationBound, 0.0f) : stage.Fbm(q, span.Velocity, span.Speed, span.AccelerationBound, 0.0f); + var result = ridged ? stage.Ridged(span.Path) : stage.Fbm(span.Path); var expected = 0.0; for (var o = 0; o < stage.Octaves; o++) { var h = snoise_hessian(stage.Point(q, o)); var frobenius = Math.Sqrt(Enumerable.Range(0, 9).Sum(i => (double)h[i / 3][i % 3] * h[i / 3][i % 3])); - expected += stage.Amplitude(o) * stage.Scale(o) * stage.Scale(o) * stage.Norm() * stage.Norm() * Math.Min(SNOISE_HESSIAN, frobenius); + var slack = Math.Pow(2.0, -10) + SNOISE_THIRD * Math.Pow(2.0, -20) * stage.Magnitude(span.Path, o); + expected += Math.Abs(stage.Amplitude(o)) * stage.Scale(o) * stage.Scale(o) * stage.Norm() * stage.Norm() * Math.Min(SNOISE_HESSIAN, frobenius + slack); } Assert.InRange(result.curvature, expected * (1.0 - 1.0e-3), expected * (1.0 + 1.0e-3)); @@ -357,27 +459,50 @@ public sealed class SpanStageTests // Through the Hessian witness: one octave at the identity transform, the span's midpoint 0.1 to 0.5 away, // moving toward the witness at speed 1, so the witness lies on the span. var witness = new float3(-13.7920685f, -9.672134f, 61.563934f); - var identity = new Stage(1.0f, 1.0f, new float3(0.0f, 0.0f, 0.0f), new float3(0.0f, 0.0f, 0.0f), new float3(0.0f, 0.0f, 0.0f), new float3(0.0f, 0.0f, 0.0f), 1, 2.0f, 0.5f); + var identity = new Stage(SpanStage.Identity, 1, 2.0f, 0.5f); var through = 0.0; for (var k = 0; k < 200; k++) { var u = UnitVector(random); var distance = Uniform(random, 0.1f, 0.5f); - var span = new Span(witness + u * distance, -u, new float3(0.0f, 0.0f, 0.0f), 1.0f, 0.0f, 0.6f); - var result = ridged ? identity.Ridged(span.Q, span.Velocity, 1.0f, 0.0f, 0.6f) : identity.Fbm(span.Q, span.Velocity, 1.0f, 0.0f, 0.6f); + var span = new Span(new SpanPath { midpoint = witness + u * distance, velocity = -u, speed = 1.0f, acceleration = 0.0f, halfLength = 0.6f }, None); + var result = ridged ? identity.Ridged(span.Path) : identity.Fbm(span.Path); CheckCurvature(identity, ridged, span, result, 64, ref tightest, ref checkedPoints); through = Math.Max(through, SymmetricNorm(FieldHessian(identity, witness, ridged))); } Assert.True(through > 40.0, $"the witness Hessian norm is {through:R}: the through-the-witness spans do not reach the curvature"); - output.WriteLine($"IV-REPORT span curvature ({(ridged ? "ridged" : "fbm")}): {checkedPoints} points enclosed, least relative slack {tightest:R}, witness norm {through:R}"); + output.WriteLine($"IV-REPORT span curvature ({(ridged ? "ridged" : "fbm")}): {checkedPoints} points enclosed, bands {string.Join("/", buckets)}, least relative slack {tightest:R}, witness norm {through:R}"); + } + + /// + /// local-tier-coefficients-unpinned: a fixed span (found by an audit, kept) at SpanStage.Identity, one + /// octave, where the Hessian's operator norm at the end of the span (11.768) exceeds the local tier built + /// with half of SNOISE_THIRD's term (about 10.1) and the exact norm along the 64 sampled points stays + /// under curvature (12.526 before the allowance). + /// + [Theory] + [InlineData(false)] + [InlineData(true)] + public void SpanCurvatureEnclosesThePinnedRankOneSpan(bool ridged) + { + var stage = new Stage(SpanStage.Identity, 1, 2.0f, 0.5f); + var span = new Span(new SpanPath { midpoint = new float3(144.93297f, -286.08800f, -53.59360f), velocity = new float3(-0.02139f, 0.17122f, -0.98500f), speed = 1.0f, acceleration = 0.0f, halfLength = 0.01f }, None); + var result = ridged ? stage.Ridged(span.Path) : stage.Fbm(span.Path); + var tightest = 1.0; + var points = 0; + CheckCurvature(stage, ridged, span, result, 64, ref tightest, ref points); + var end = SymmetricNorm(FieldHessian(stage, span.At(span.Path.halfLength), ridged)); + Assert.True(end > 10.5 && result.curvature >= end, $"the end-of-span norm is {end:R} and curvature {result.curvature:R}: the span no longer pins the local tier"); } /// /// span-bounds-enclose, fold: for span_ridged, the value's departure from the midpoint branch (each /// octave's gain^o (s n - |n|), s the sign of n at the midpoint, summed) never exceeds fold at 64 points of - /// every span; fold is positive on spans that do cross zero (and the departure there is not trivially - /// zero), exactly 0 when no octave can reach zero and for span_fbm, and 0 at halfLength 0. + /// every span, measured twice: for the float field span_ridged itself evaluates at the path point, and for + /// the exact octave sum (double kernel); in every coordinate band from S_o under 1e2 to over 1e6. fold is + /// positive on spans that do cross zero (and the departure there is not trivially zero), exactly 0 for + /// span_fbm and when no octave can reach zero. /// [Fact] public void SpanFoldEnclosesTheBranchDeviation() @@ -385,62 +510,77 @@ public sealed class SpanStageTests var random = new System.Random(0xF01D); var crossing = 0; var largest = 0.0; - var tightest = 1.0; - for (var k = 0; k < 400; k++) - { - var stage = Stage.Draw(random, random.Next(1, 4)); - var halfLength = new[] { 0.02f, 0.1f, 0.5f, 2.0f }[random.Next(4)]; - // Draw midpoints with small noise values so spans cross zero, as well as arbitrary ones. - var q = Vector(random, -30.0f, 30.0f); - if (k % 2 == 0) + var tightestFloat = 1.0; + var tightestExact = 1.0; + var buckets = new int[6]; + foreach (var offsetScale in OffsetScales) + for (var k = 0; k < 60; k++) { - var best = q; - var smallest = float.MaxValue; - for (var trial = 0; trial < 40; trial++) + var stage = Banded(random, offsetScale, 3); + var halfLength = new[] { 0.02f, 0.1f, 0.5f, 2.0f }[random.Next(4)]; + // Draw midpoints with small noise values so spans cross zero, as well as arbitrary ones. + var q = Vector(random, -30.0f, 30.0f); + if (k % 2 == 0) { - var candidate = Vector(random, -30.0f, 30.0f); - var value = abs(snoise(stage.Point(candidate, 0))); - if (value < smallest) - (smallest, best) = (value, candidate); + var best = q; + var smallest = float.MaxValue; + for (var trial = 0; trial < 40; trial++) + { + var candidate = Vector(random, -30.0f, 30.0f); + var value = abs(snoise(stage.Point(candidate, 0))); + if (value < smallest) + (smallest, best) = (value, candidate); + } + + q = best; } - q = best; - } - - var span = Span.Draw(random, q, halfLength, Uniform(random, 0.0f, 2.0f)); - var fold = stage.Ridged(q, span.Velocity, span.Speed, span.AccelerationBound, halfLength).fold; - Assert.Equal(0.0f, stage.Fbm(q, span.Velocity, span.Speed, span.AccelerationBound, halfLength).fold); - Assert.Equal(0.0f, stage.Ridged(q, span.Velocity, span.Speed, span.AccelerationBound, 0.0f).fold); - var signs = Enumerable.Range(0, stage.Octaves).Select(o => Math.Sign(snoise(stage.Point(q, o)))).ToArray(); - var peak = 0.0; - for (var s = 0; s < 64; s++) - { - var t = halfLength * (2.0f * s / 63 - 1.0f); - var point = span.At(t); - var departure = 0.0; - for (var o = 0; o < stage.Octaves; o++) + var span = Span.Draw(random, q, halfLength, Uniform(random, 0.0f, 2.0f)); + buckets[Bucket(stage.MaxMagnitude(span.Path))]++; + var fold = stage.Ridged(span.Path).fold; + Assert.Equal(0.0f, stage.Fbm(span.Path).fold); + var signs = Enumerable.Range(0, stage.Octaves).Select(o => Math.Sign(stage.OneOctave(o).Fbm(PathAt(q)).field.w)).ToArray(); + var peakFloat = 0.0; + var peakExact = 0.0; + for (var s = 0; s < 64; s++) { - var n = (double)snoise(stage.Point(point, o)); - departure += stage.Amplitude(o) * (signs[o] * n - Math.Abs(n)); + var t = halfLength * (2.0f * s / 63 - 1.0f); + var point = span.At(t); + var departureFloat = 0.0; + var departureExact = 0.0; + for (var o = 0; o < stage.Octaves; o++) + { + var octave = stage.OneOctave(o); + var nFloat = (double)octave.Fbm(PathAt(point)).field.w; + var absFloat = 1.0 - octave.Ridged(PathAt(point)).field.w; + departureFloat += stage.Amplitude(o) * (signs[o] * nFloat - absFloat); + var x = stage.PointD(point, o); + var nExact = NoiseBoundTests.ValueDouble(x[0], x[1], x[2]); + departureExact += stage.Amplitude(o) * (signs[o] * nExact - Math.Abs(nExact)); + } + + peakFloat = Math.Max(peakFloat, Math.Abs(departureFloat)); + peakExact = Math.Max(peakExact, Math.Abs(departureExact)); + Assert.True(Math.Abs(departureFloat) <= fold * (1.0 + 1.0e-3) + 1.0e-6, $"stage {stage}, span {span}: the float branch departure at t = {t:R} is {departureFloat:R}, above fold {fold:R}"); + Assert.True(Math.Abs(departureExact) <= fold * (1.0 + 1.0e-3) + 1.0e-6, $"stage {stage}, span {span}: the exact branch departure at t = {t:R} is {departureExact:R}, above fold {fold:R}"); } - peak = Math.Max(peak, Math.Abs(departure)); - Assert.True(Math.Abs(departure) <= fold * (1.0 + 1.0e-3), $"stage {stage}, span {span}: the branch departure at t = {t:R} is {departure:R}, above fold {fold:R}"); + if (fold > 0.0) + crossing++; + largest = Math.Max(largest, Math.Max(peakFloat, peakExact)); + if (fold > 0.0 && peakFloat > 0.0) + tightestFloat = Math.Min(tightestFloat, peakFloat / fold); + if (fold > 0.0 && peakExact > 0.0) + tightestExact = Math.Min(tightestExact, peakExact / fold); } - if (fold > 0.0) - crossing++; - largest = Math.Max(largest, peak); - if (fold > 0.0 && peak > 0.0) - tightest = Math.Min(tightest, peak / fold); - } - - Assert.True(crossing >= 100, $"only {crossing} of 400 spans had a positive fold"); + Assert.True(buckets.All(count => count >= 3), $"the coordinate bands drawn are {string.Join(", ", buckets)} spans: a band is not reached"); + Assert.True(crossing >= 100, $"only {crossing} of {OffsetScales.Length * 60} spans had a positive fold"); Assert.True(largest > 0.1, $"the largest branch departure drawn is {largest:R}: the check has no teeth"); - output.WriteLine($"IV-REPORT span fold: {crossing} spans with fold > 0, largest departure {largest:R}, tightest departure / fold {tightest:R}"); + output.WriteLine($"IV-REPORT span fold: {crossing} spans with fold > 0, bands {string.Join("/", buckets)}, largest departure {largest:R}, tightest departure / fold {tightestFloat:R} (float) {tightestExact:R} (exact)"); // An octave that cannot reach zero: a midpoint where |n| is large and a short span give fold exactly 0. - var calm = new Stage(1.0f, 1.0f, new float3(0.0f, 0.0f, 0.0f), new float3(0.0f, 0.0f, 0.0f), new float3(0.0f, 0.0f, 0.0f), new float3(0.0f, 0.0f, 0.0f), 1, 2.0f, 0.5f); + var calm = new Stage(SpanStage.Identity, 1, 2.0f, 0.5f); var found = 0; for (var k = 0; k < 400 && found < 20; k++) { @@ -448,7 +588,7 @@ public sealed class SpanStageTests if (abs(snoise(q)) < 0.6f) continue; found++; - Assert.Equal(0.0f, calm.Ridged(q, new float3(1.0f, 0.0f, 0.0f), 1.0f, 0.0f, 0.01f).fold); + Assert.Equal(0.0f, calm.Ridged(new SpanPath { midpoint = q, velocity = new float3(1.0f, 0.0f, 0.0f), speed = 1.0f, acceleration = 0.0f, halfLength = 0.01f }).fold); } Assert.True(found >= 5, $"only {found} calm midpoints found"); @@ -464,8 +604,7 @@ public sealed class SpanStageTests public void SpanFoldCoversAccelerationAlone() { var random = new System.Random(0xACCE1); - var identity = new Stage(1.0f, 1.0f, new float3(0.0f, 0.0f, 0.0f), new float3(0.0f, 0.0f, 0.0f), new float3(0.0f, 0.0f, 0.0f), new float3(0.0f, 0.0f, 0.0f), 1, 2.0f, 0.5f); - var still = new float3(0.0f, 0.0f, 0.0f); + var identity = new Stage(SpanStage.Identity, 1, 2.0f, 0.5f); var cases = 0; var crossings = 0; var tightest = double.MaxValue; @@ -480,8 +619,8 @@ public sealed class SpanStageTests if (length(gradient) < 1.5f || abs(n.w) > 0.4f * length(gradient) * acceleration * halfLength * halfLength * 0.5f) continue; var c = -(float)Math.Sign(n.w) * acceleration * normalize(gradient); - var span = new Span(q, still, c, acceleration * halfLength, acceleration, halfLength); - var fold = identity.Ridged(q, still, span.Speed, acceleration, halfLength).fold; + var span = new Span(new SpanPath { midpoint = q, velocity = None, speed = acceleration * halfLength, acceleration = acceleration, halfLength = halfLength }, c); + var fold = identity.Ridged(span.Path).fold; var departure = 0.0; for (var s = 0; s < 64; s++) {