fix(simd): take the ssimulacra2 edge-diff difference in double, and gate ISA invariance

CHANGELOG.md

@@ -25852,6 +25852,31 @@ mirroring the `sse_line_16_c` reference in `core/src/feature/integer_psnr.c`.
   PR).
 
 
+- **`ssimulacra2` scored differently on a host with SIMD than on one
+  without.** The four `edge_diff_map` kernels (AVX2, AVX-512, NEON, SVE2)
+  computed the per-pixel difference `|img - blur(img)|` with a *float*
+  subtract and promoted to double afterwards, while the scalar reference —
+  and each kernel's own scalar tail — promote first and subtract in double,
+  which is exact for two floats. Because the ssimulacra2 pipeline is
+  ill-conditioned downstream (that difference is a catastrophic cancellation
+  and pooling takes a 4-norm), the rounding surfaced as a ~1.6e-09 score
+  difference between a SIMD host and a scalar one. The per-feature SIMD test
+  could not see it: it compares each kernel against a scalar reference defined
+  inside the test file, at a 33x21 fixture on which the float subtraction
+  happens to be exact. Fixed by folding the subtraction into the per-lane loop
+  in double, matching the reference exactly. The fork-added
+  `ssimulacra2_test.py` snapshot is unchanged to six decimals, so no snapshot
+  was regenerated.
+
+- **New gate: a score may not depend on the host instruction set.**
+  `core/test/test_feature_isa_invariance.c` drives the public API twice over
+  one fixture — once with the host ISA, once with `cpumask` disabling every
+  SIMD flag — and asserts bit-identical scores across ten features. It uses no
+  internal symbols, so unlike the per-feature SIMD tests it cannot drift away
+  from the shipped code. It found the `edge_diff_map` defect above on its
+  first run; the other nine features already passed.
+
+
 Remove two stale `//TODO: dedupe` comments from `predict.c` and `libvmaf.c`
 that pointed at each other. The deduplication was completed in PR #1067
 (ADR-0480, `bootstrap_names.h` shared header); replace each marker with a

changelog.d/fixed/ssimulacra2-edge-diff-double-subtract.md

@@ -0,0 +1,23 @@
+- **`ssimulacra2` scored differently on a host with SIMD than on one
+  without.** The four `edge_diff_map` kernels (AVX2, AVX-512, NEON, SVE2)
+  computed the per-pixel difference `|img - blur(img)|` with a *float*
+  subtract and promoted to double afterwards, while the scalar reference —
+  and each kernel's own scalar tail — promote first and subtract in double,
+  which is exact for two floats. Because the ssimulacra2 pipeline is
+  ill-conditioned downstream (that difference is a catastrophic cancellation
+  and pooling takes a 4-norm), the rounding surfaced as a ~1.6e-09 score
+  difference between a SIMD host and a scalar one. The per-feature SIMD test
+  could not see it: it compares each kernel against a scalar reference defined
+  inside the test file, at a 33x21 fixture on which the float subtraction
+  happens to be exact. Fixed by folding the subtraction into the per-lane loop
+  in double, matching the reference exactly. The fork-added
+  `ssimulacra2_test.py` snapshot is unchanged to six decimals, so no snapshot
+  was regenerated.
+
+- **New gate: a score may not depend on the host instruction set.**
+  `core/test/test_feature_isa_invariance.c` drives the public API twice over
+  one fixture — once with the host ISA, once with `cpumask` disabling every
+  SIMD flag — and asserts bit-identical scores across ten features. It uses no
+  internal symbols, so unlike the per-feature SIMD tests it cannot drift away
+  from the shipped code. It found the `edge_diff_map` defect above on its
+  first run; the other nine features already passed.

core/src/feature/AGENTS.md

@@ -1576,3 +1576,24 @@ gain limit, so it moves `adm` scores directly. The four historical spellings
 of this predicate disagreed on about 4e-5 of near-parallel scale-0 band
 quadruples; see
 [`docs/research/2030-adm-angle-flag-fp64-free.md`](../../../docs/research/2030-adm-angle-flag-fp64-free.md).
+
+## A score must not depend on the host ISA (ADR-1207, ADR-1208)
+
+Every SIMD kernel in this tree is required to be bit-exact with its scalar
+reference, so the same input must produce the same bits on an AVX-512 host, an
+AVX2 host and a host with no SIMD. Two things follow for anyone touching a
+kernel here:
+
+- The per-feature `test_<feature>_simd.c` files compare a SIMD kernel against a
+  scalar reference **defined inside the test TU**, because the shipped scalar
+  functions are `static`. That reference can drift away from the shipped one —
+  it did, twice (ADR-1205, ADR-1208). Passing `test_<feature>_simd` is
+  therefore necessary but not sufficient.
+- `core/test/test_feature_isa_invariance.c` is the gate that actually compares
+  the shipped SIMD path against the shipped scalar path, end to end, via
+  `VmafConfiguration.cpumask`. Run it after any kernel change.
+
+Concretely, when a kernel promotes `float` inputs to `double`, do the promotion
+**before** the arithmetic, not after. `(double)a - (double)b` is exact for two
+floats; `(double)(a - b)` is not, and mixing the two between a vector body and
+its scalar tail makes the result depend on the vector width.

core/src/feature/arm64/ssimulacra2_neon.c

@@ -325,15 +325,25 @@ void ssimulacra2_edge_diff_map_neon(const float *img1, const float *mu1, const f
             const float32x4_t a2 = vld1q_f32(r2 + i);
             const float32x4_t am1 = vld1q_f32(rm1 + i);
             const float32x4_t am2 = vld1q_f32(rm2 + i);
-            const float32x4_t d1 = vabsq_f32(vsubq_f32(a1, am1));
-            const float32x4_t d2 = vabsq_f32(vsubq_f32(a2, am2));
-            alignas(16) float d1f[4];
-            alignas(16) float d2f[4];
-            vst1q_f32(d1f, d1);
-            vst1q_f32(d2f, d2);
+            /* ADR-1208: the reference difference is taken in DOUBLE.
+             * `a` and `am` are floats, so `(double)a - (double)am` is exact,
+             * whereas subtracting in float rounds first. The scalar
+             * `edge_diff_map`, this function's own scalar tail, and the test's
+             * reference all promote before subtracting; vectorising the
+             * subtract in float made the ssimulacra2 score depend on whether
+             * the host had SIMD. The per-lane loop below is scalar anyway, so
+             * nothing is lost by folding the subtraction into it. */
+            alignas(16) float a1f[4];
+            alignas(16) float am1f[4];
+            alignas(16) float a2f[4];
+            alignas(16) float am2f[4];
+            vst1q_f32(a1f, a1);
+            vst1q_f32(am1f, am1);
+            vst1q_f32(a2f, a2);
+            vst1q_f32(am2f, am2);
             for (int k = 0; k < 4; k++) {
-                double ed1 = (double)d1f[k];
-                double ed2 = (double)d2f[k];
+                double ed1 = fabs((double)a1f[k] - (double)am1f[k]);
+                double ed2 = fabs((double)a2f[k] - (double)am2f[k]);
                 double d = (1.0 + ed2) / (1.0 + ed1) - 1.0;
                 double art = d > 0.0 ? d : 0.0;
                 double det = d < 0.0 ? -d : 0.0;

core/src/feature/arm64/ssimulacra2_sve2.c

@@ -361,15 +361,25 @@ void ssimulacra2_edge_diff_map_sve2(const float *img1, const float *mu1, const f
             const svfloat32_t a2 = svld1_f32(pg, r2 + i);
             const svfloat32_t am1 = svld1_f32(pg, rm1 + i);
             const svfloat32_t am2 = svld1_f32(pg, rm2 + i);
-            const svfloat32_t d1 = svabs_f32_x(pg, svsub_f32_x(pg, a1, am1));
-            const svfloat32_t d2 = svabs_f32_x(pg, svsub_f32_x(pg, a2, am2));
-            alignas(16) float d1f[4] = {0.f, 0.f, 0.f, 0.f};
-            alignas(16) float d2f[4] = {0.f, 0.f, 0.f, 0.f};
-            svst1_f32(pg, d1f, d1);
-            svst1_f32(pg, d2f, d2);
+            /* ADR-1208: the reference difference is taken in DOUBLE.
+             * `a` and `am` are floats, so `(double)a - (double)am` is exact,
+             * whereas subtracting in float rounds first. The scalar
+             * `edge_diff_map`, this function's own scalar tail, and the test's
+             * reference all promote before subtracting; vectorising the
+             * subtract in float made the ssimulacra2 score depend on whether
+             * the host had SIMD. The per-lane loop below is scalar anyway, so
+             * nothing is lost by folding the subtraction into it. */
+            alignas(16) float a1f[4] = {0.f, 0.f, 0.f, 0.f};
+            alignas(16) float am1f[4] = {0.f, 0.f, 0.f, 0.f};
+            alignas(16) float a2f[4] = {0.f, 0.f, 0.f, 0.f};
+            alignas(16) float am2f[4] = {0.f, 0.f, 0.f, 0.f};
+            svst1_f32(pg, a1f, a1);
+            svst1_f32(pg, am1f, am1);
+            svst1_f32(pg, a2f, a2);
+            svst1_f32(pg, am2f, am2);
             for (int k = 0; k < 4; k++) {
-                double ed1 = (double)d1f[k];
-                double ed2 = (double)d2f[k];
+                double ed1 = fabs((double)a1f[k] - (double)am1f[k]);
+                double ed2 = fabs((double)a2f[k] - (double)am2f[k]);
                 double d = (1.0 + ed2) / (1.0 + ed1) - 1.0;
                 double art = d > 0.0 ? d : 0.0;
                 double det = d < 0.0 ? -d : 0.0;

core/src/feature/x86/ssimulacra2_avx2.c

@@ -374,7 +374,6 @@ void ssimulacra2_edge_diff_map_avx2(const float *img1, const float *mu1, const f
 {
     const size_t plane = (size_t)w * (size_t)h;
     const double one_per_pixels = 1.0 / (double)plane;
-    const __m256 vsignmask = _mm256_castsi256_ps(_mm256_set1_epi32(0x7FFFFFFF));
 
     for (int c = 0; c < 3; c++) {
         double s0 = 0.0;
@@ -392,16 +391,25 @@ void ssimulacra2_edge_diff_map_avx2(const float *img1, const float *mu1, const f
             const __m256 a2 = _mm256_loadu_ps(r2 + i);
             const __m256 am1 = _mm256_loadu_ps(rm1 + i);
             const __m256 am2 = _mm256_loadu_ps(rm2 + i);
-            const __m256 d1 = _mm256_and_ps(vsignmask, _mm256_sub_ps(a1, am1));
-            const __m256 d2 = _mm256_and_ps(vsignmask, _mm256_sub_ps(a2, am2));
-            /* Promote to double per-lane for the divide, matching scalar. */
-            alignas(32) float d1f[8];
-            alignas(32) float d2f[8];
-            _mm256_store_ps(d1f, d1);
-            _mm256_store_ps(d2f, d2);
+            /* ADR-1208: the reference difference is taken in DOUBLE.
+             * `a` and `am` are floats, so `(double)a - (double)am` is exact,
+             * whereas subtracting in float rounds first. The scalar
+             * `edge_diff_map`, this function's own scalar tail, and the test's
+             * reference all promote before subtracting; vectorising the
+             * subtract in float made the ssimulacra2 score depend on whether
+             * the host had SIMD. The per-lane loop below is scalar anyway, so
+             * nothing is lost by folding the subtraction into it. */
+            alignas(32) float a1f[8];
+            alignas(32) float am1f[8];
+            alignas(32) float a2f[8];
+            alignas(32) float am2f[8];
+            _mm256_store_ps(a1f, a1);
+            _mm256_store_ps(am1f, am1);
+            _mm256_store_ps(a2f, a2);
+            _mm256_store_ps(am2f, am2);
             for (int k = 0; k < 8; k++) {
-                double ed1 = (double)d1f[k];
-                double ed2 = (double)d2f[k];
+                double ed1 = fabs((double)a1f[k] - (double)am1f[k]);
+                double ed2 = fabs((double)a2f[k] - (double)am2f[k]);
                 double d = (1.0 + ed2) / (1.0 + ed1) - 1.0;
                 double art = d > 0.0 ? d : 0.0;
                 double det = d < 0.0 ? -d : 0.0;

core/src/feature/x86/ssimulacra2_avx512.c

@@ -318,7 +318,6 @@ void ssimulacra2_edge_diff_map_avx512(const float *img1, const float *mu1, const
 {
     const size_t plane = (size_t)w * (size_t)h;
     const double one_per_pixels = 1.0 / (double)plane;
-    const __m512 vsignmask = _mm512_castsi512_ps(_mm512_set1_epi32(0x7FFFFFFF));
 
     for (int c = 0; c < 3; c++) {
         double s0 = 0.0;
@@ -336,15 +335,25 @@ void ssimulacra2_edge_diff_map_avx512(const float *img1, const float *mu1, const
             const __m512 a2 = _mm512_loadu_ps(r2 + i);
             const __m512 am1 = _mm512_loadu_ps(rm1 + i);
             const __m512 am2 = _mm512_loadu_ps(rm2 + i);
-            const __m512 d1 = _mm512_and_ps(vsignmask, _mm512_sub_ps(a1, am1));
-            const __m512 d2 = _mm512_and_ps(vsignmask, _mm512_sub_ps(a2, am2));
-            alignas(64) float d1f[16];
-            alignas(64) float d2f[16];
-            _mm512_store_ps(d1f, d1);
-            _mm512_store_ps(d2f, d2);
+            /* ADR-1208: the reference difference is taken in DOUBLE.
+             * `a` and `am` are floats, so `(double)a - (double)am` is exact,
+             * whereas subtracting in float rounds first. The scalar
+             * `edge_diff_map`, this function's own scalar tail, and the test's
+             * reference all promote before subtracting; vectorising the
+             * subtract in float made the ssimulacra2 score depend on whether
+             * the host had SIMD. The per-lane loop below is scalar anyway, so
+             * nothing is lost by folding the subtraction into it. */
+            alignas(64) float a1f[16];
+            alignas(64) float am1f[16];
+            alignas(64) float a2f[16];
+            alignas(64) float am2f[16];
+            _mm512_store_ps(a1f, a1);
+            _mm512_store_ps(am1f, am1);
+            _mm512_store_ps(a2f, a2);
+            _mm512_store_ps(am2f, am2);
             for (int k = 0; k < 16; k++) {
-                double ed1 = (double)d1f[k];
-                double ed2 = (double)d2f[k];
+                double ed1 = fabs((double)a1f[k] - (double)am1f[k]);
+                double ed2 = fabs((double)a2f[k] - (double)am2f[k]);
                 double d = (1.0 + ed2) / (1.0 + ed1) - 1.0;
                 double art = d > 0.0 ? d : 0.0;
                 double det = d < 0.0 ? -d : 0.0;

core/test/meson.build

@@ -1325,6 +1325,21 @@ _ssimulacra2_test_x86_fma_args = []
 if host_machine.cpu_family() in ['x86_64', 'x86']
     _ssimulacra2_test_x86_fma_args = ['-mavx2', '-mfma']
 endif
+# ADR-1207 — a score must not depend on the host instruction set.
+# Drives the public API twice over one fixture, once with the host ISA and
+# once with cpumask disabling every SIMD flag, and asserts bit-identical
+# scores. Complements the per-feature *_simd tests, which compare a SIMD
+# kernel against a scalar reference defined inside the test TU rather than
+# against the shipped one.
+test_feature_isa_invariance = executable('test_feature_isa_invariance',
+    ['test.c', 'test_feature_isa_invariance.c'],
+    include_directories : [libvmaf_inc, test_inc],
+    link_with : get_option('default_library') == 'both' ? libvmaf.get_static_lib() : libvmaf,
+    dependencies : [pthread_dependency, math_lib],
+)
+test('test_feature_isa_invariance', test_feature_isa_invariance,
+    suite : ['fast'])
+
 test_ssimulacra2_simd = executable('test_ssimulacra2_simd',
     ['test.c', 'test_ssimulacra2_simd.c',
      '../src/mem.cpp', '../src/picture.c', '../src/ref.cpp',