Merge pull request #2909 from Xarthisius/minor_fix_for_2731

Fix indexing issues in sph pixelization

tests/tests.yaml

@@ -12,7 +12,7 @@ answer_tests:
     - yt/frontends/amrvac/tests/test_outputs.py:test_riemann_cartesian_175D
     - yt/frontends/amrvac/tests/test_outputs.py:test_rmi_cartesian_dust_2D
 
-  local_arepo_006:
+  local_arepo_007:  # PR 2909
     - yt/frontends/arepo/tests/test_outputs.py:test_arepo_bullet
     - yt/frontends/arepo/tests/test_outputs.py:test_arepo_tng59
     - yt/frontends/arepo/tests/test_outputs.py:test_index_override
@@ -59,7 +59,7 @@ answer_tests:
     - yt/frontends/flash/tests/test_outputs.py:test_wind_tunnel
     - yt/frontends/flash/tests/test_outputs.py:test_fid_1to3_b1
 
-  local_gadget_004:
+  local_gadget_005:  # PR 2909
     - yt/frontends/gadget/tests/test_outputs.py:test_iso_collapse
     - yt/frontends/gadget/tests/test_outputs.py:test_pid_uniqueness
     - yt/frontends/gadget/tests/test_outputs.py:test_bigendian_field_access
@@ -73,7 +73,7 @@ answer_tests:
   local_gdf_001:
     - yt/frontends/gdf/tests/test_outputs.py:test_sedov_tunnel
 
-  local_gizmo_005:
+  local_gizmo_006:  # PR 2909
     - yt/frontends/gizmo/tests/test_outputs.py:test_gizmo_64
 
   local_halos_009:
@@ -82,7 +82,7 @@ answer_tests:
     - yt/frontends/gadget_fof/tests/test_outputs.py:test_fields_g5
     - yt/frontends/gadget_fof/tests/test_outputs.py:test_fields_g42
 
-  local_owls_005:
+  local_owls_006:  # PR 2909
     - yt/frontends/owls/tests/test_outputs.py:test_snapshot_033
     - yt/frontends/owls/tests/test_outputs.py:test_OWLS_particlefilter
 
@@ -94,7 +94,7 @@ answer_tests:
     - yt/visualization/tests/test_particle_plot.py:test_particle_phase_answers
     - yt/visualization/tests/test_raw_field_slices.py:test_raw_field_slices
 
-  local_tipsy_006:
+  local_tipsy_007:  # PR 2909
     - yt/frontends/tipsy/tests/test_outputs.py:test_pkdgrav
     - yt/frontends/tipsy/tests/test_outputs.py:test_gasoline_dmonly
     - yt/frontends/tipsy/tests/test_outputs.py:test_tipsy_galaxy

yt/utilities/answer_testing/framework.py

@@ -36,6 +36,7 @@
 from yt.utilities.exceptions import YTCloudError, YTNoAnswerNameSpecified, YTNoOldAnswer
 from yt.utilities.logger import disable_stream_logging
 from yt.visualization import (
+    image_writer as image_writer,
     particle_plots as particle_plots,
     plot_window as pw,
     profile_plotter as profile_plotter,
@@ -669,7 +670,11 @@ def compare(self, new_result, old_result):
             # weight_field does not have units, so we do not directly compare them
             if k == "weight_field_sum":
                 continue
-            assert_allclose_units(new_result[k], old_result[k], 1e-10)
+            try:
+                assert_allclose_units(new_result[k], old_result[k], 1e-10)
+            except AssertionError:
+                dump_images(new_result[k], old_result[k])
+                raise
 
 
 class PixelizedParticleProjectionValuesTest(PixelizedProjectionValuesTest):
@@ -782,6 +787,23 @@ def compare(self, new_result, old_result):
             assert newc == oldc
 
 
+def dump_images(new_result, old_result, decimals=10):
+    tmpfd, old_image = tempfile.mkstemp(suffix=".png")
+    os.close(tmpfd)
+    tmpfd, new_image = tempfile.mkstemp(suffix=".png")
+    os.close(tmpfd)
+    image_writer.write_projection(new_result, new_image)
+    image_writer.write_projection(old_result, old_image)
+    results = compare_images(old_image, new_image, 10 ** (-decimals))
+    if results is not None:
+        tempfiles = [
+            line.strip() for line in results.split("\n") if line.endswith(".png")
+        ]
+        for fn in tempfiles:
+            sys.stderr.write(f"\n[[ATTACHMENT|{fn}]]")
+        sys.stderr.write("\n")
+
+
 def compare_image_lists(new_result, old_result, decimals):
     fns = []
     for _ in range(2):

yt/utilities/lib/pixelization_routines.pyx

@@ -997,7 +997,7 @@ def pixelize_sph_kernel_projection(
     cdef np.float64_t q_ij2, posx_diff, posy_diff, ih_j2
     cdef np.float64_t x, y, dx, dy, idx, idy, h_j2, px, py
     cdef np.float64_t period_x, period_y
-    cdef int index, i, j
+    cdef int index, i, j, ii, jj
     cdef np.float64_t[:] _weight_field
     cdef int xiter[2]
     cdef int yiter[2]
@@ -1040,30 +1040,27 @@ def pixelize_sph_kernel_projection(
 
             xiter[1] = yiter[1] = 999
 
-            px = posx[j]
-            py = posy[j]
-
             if check_period == 1:
-                if px - hsml[j] < x_min:
+                if posx[j] - hsml[j] < x_min:
                     xiter[1] = +1
                     xiterv[1] = period_x
-                elif px + hsml[j] > x_max:
+                elif posx[j] + hsml[j] > x_max:
                     xiter[1] = -1
                     xiterv[1] = -period_x
-                if py - hsml[j] < y_min:
+                if posy[j] - hsml[j] < y_min:
                     yiter[1] = +1
                     yiterv[1] = period_y
-                elif py + hsml[j] > y_max:
+                elif posy[j] + hsml[j] > y_max:
                     yiter[1] = -1
                     yiterv[1] = -period_y
 
-            for xi in range(2):
-                if xiter[xi] == 999: continue
-                px += xiterv[xi]
+            for ii in range(2):
+                if xiter[ii] == 999: continue
+                px = posx[j] + xiterv[ii]
                 if (px + hsml[j] < x_min) or (px - hsml[j] > x_max): continue
-                for yi in range(2):
-                    if yiter[yi] == 999: continue
-                    py += yiterv[yi]
+                for jj in range(2):
+                    if yiter[jj] == 999: continue
+                    py = posy[j] + yiterv[jj]
                     if (py + hsml[j] < y_min) or (py - hsml[j] > y_max): continue
 
                     # here we find the pixels which this particle contributes to
@@ -1296,7 +1293,7 @@ def pixelize_sph_kernel_slice(
     cdef np.int64_t xi, yi, x0, x1, y0, y1
     cdef np.float64_t q_ij, posx_diff, posy_diff, ih_j
     cdef np.float64_t x, y, dx, dy, idx, idy, h_j2, h_j, px, py
-    cdef int index, i, j
+    cdef int index, i, j, ii, jj
     cdef np.float64_t period_x, period_y
     cdef int xiter[2]
     cdef int yiter[2]
@@ -1331,30 +1328,27 @@ def pixelize_sph_kernel_slice(
 
             xiter[1] = yiter[1] = 999
 
-            px = posx[j]
-            py = posy[j]
-
             if check_period == 1:
-                if px - hsml[j] < x_min:
+                if posx[j] - hsml[j] < x_min:
                     xiter[1] = +1
                     xiterv[1] = period_x
-                elif px + hsml[j] > x_max:
+                elif posx[j] + hsml[j] > x_max:
                     xiter[1] = -1
                     xiterv[1] = -period_x
-                if py - hsml[j] < y_min:
+                if posy[j] - hsml[j] < y_min:
                     yiter[1] = +1
                     yiterv[1] = period_y
-                elif py + hsml[j] > y_max:
+                elif posy[j] + hsml[j] > y_max:
                     yiter[1] = -1
                     yiterv[1] = -period_y
 
-            for xi in range(2):
-                if xiter[xi] == 999: continue
-                px += xiterv[xi]
+            for ii in range(2):
+                if xiter[ii] == 999: continue
+                px = posx[j] + xiterv[ii]
                 if (px + hsml[j] < x_min) or (px - hsml[j] > x_max): continue
-                for yi in range(2):
-                    if yiter[yi] == 999: continue
-                    py += yiterv[yi]
+                for jj in range(2):
+                    if yiter[jj] == 999: continue
+                    py = posy[j] + yiterv[jj]
                     if (py + hsml[j] < y_min) or (py - hsml[j] > y_max): continue
 
                     x0 = <np.int64_t> ( (px - hsml[j] - x_min) * idx)
@@ -1417,7 +1411,7 @@ def pixelize_sph_kernel_arbitrary_grid(np.float64_t[:, :, :] buff,
     cdef np.int64_t xi, yi, zi, x0, x1, y0, y1, z0, z1
     cdef np.float64_t q_ij, posx_diff, posy_diff, posz_diff, px, py, pz
     cdef np.float64_t x, y, z, dx, dy, dz, idx, idy, idz, h_j3, h_j2, h_j, ih_j
-    cdef int index, i, j, k
+    cdef int index, i, j, k, ii, jj, kk
     cdef np.float64_t period_x, period_y, period_z
 
     cdef int xiter[2]
@@ -1462,41 +1456,37 @@ def pixelize_sph_kernel_arbitrary_grid(np.float64_t[:, :, :] buff,
 
             xiter[1] = yiter[1] = ziter[1] = 999
 
-            px = posx[j]
-            py = posy[j]
-            pz = posz[j]
-
             if check_period == 1:
-                if px - hsml[j] < x_min:
+                if posx[j] - hsml[j] < x_min:
                     xiter[1] = +1
                     xiterv[1] = period_x
-                elif px + hsml[j] > x_max:
+                elif posx[j] + hsml[j] > x_max:
                     xiter[1] = -1
                     xiterv[1] = -period_x
-                if py - hsml[j] < y_min:
+                if posy[j] - hsml[j] < y_min:
                     yiter[1] = +1
                     yiterv[1] = period_y
-                elif py + hsml[j] > y_max:
+                elif posy[j] + hsml[j] > y_max:
                     yiter[1] = -1
                     yiterv[1] = -period_y
-                if pz - hsml[j] < z_min:
+                if posz[j] - hsml[j] < z_min:
                     ziter[1] = +1
                     ziterv[1] = period_z
-                elif pz + hsml[j] > z_max:
+                elif posz[j] + hsml[j] > z_max:
                     ziter[1] = -1
                     ziterv[1] = -period_z
 
-            for xi in range(2):
-                if xiter[xi] == 999: continue
-                px += xiterv[xi]
+            for ii in range(2):
+                if xiter[ii] == 999: continue
+                px = posx[j] + xiterv[ii]
                 if (px + hsml[j] < x_min) or (px - hsml[j] > x_max): continue
-                for yi in range(2):
-                    if yiter[yi] == 999: continue
-                    py += yiterv[yi]
+                for jj in range(2):
+                    if yiter[jj] == 999: continue
+                    py = posy[j] + yiterv[jj]
                     if (py + hsml[j] < y_min) or (py - hsml[j] > y_max): continue
-                    for zi in range(2):
-                        if ziter[zi] == 999: continue
-                        pz += ziterv[zi]
+                    for kk in range(2):
+                        if ziter[kk] == 999: continue
+                        pz = posz[j] + ziterv[kk]
                         if (pz + hsml[j] < z_min) or (pz - hsml[j] > z_max): continue
 
                         x0 = <np.int64_t> ( (px - hsml[j] - x_min) * idx)