[BUG] Off-axis projections for octree code lead to out-of-memory errors

[cphyc]

Bug report

Bug summary

Currently, the off-axis projection routine for octree dataset precomputes the look of a cell from the camera angle. It does so at different resolutions, ranging from cells that span ~2pixels to the largest one.

When doing a deep zoom, however, the largest cell may be thousands of pixels across after projection, which requires a correspondingly large array to be allocated, thus leading to an out-of-memory error.

Code for reproduction

import yt

ds = yt.load("output_00080")

sp = ds.sphere("max", (200, "kpc"))
sp_smol = ds.sphere(sp.center, (5, "kpc"))
sp.set_field_parameter("bulk_velocity", sp_smol.quantities.bulk_velocity(use_gas=False, use_particles=True))
n = np.array([0.99, 0.01, 0.01])
n /= np.linalg.norm(n)

p = yt.ProjectionPlot(ds, n, ("gas", "density"), center=sp.center, width=(20, "kpc"), weight_field=("gas", "dx"))

Actual outcome

yt : [INFO     ] 2024-12-06 11:51:34,873 Parameters: current_time              = 11.925285011256845 Gyr
yt : [INFO     ] 2024-12-06 11:51:34,873 Parameters: domain_dimensions         = [64 64 64]
yt : [INFO     ] 2024-12-06 11:51:34,873 Parameters: domain_left_edge          = [0. 0. 0.]
yt : [INFO     ] 2024-12-06 11:51:34,873 Parameters: domain_right_edge         = [1. 1. 1.]
yt : [INFO     ] 2024-12-06 11:51:34,874 Parameters: cosmological_simulation   = True
yt : [INFO     ] 2024-12-06 11:51:34,874 Parameters: current_redshift          = 0.14255728632206321
yt : [INFO     ] 2024-12-06 11:51:34,874 Parameters: omega_lambda              = 0.723999977111816
yt : [INFO     ] 2024-12-06 11:51:34,874 Parameters: omega_matter              = 0.276000022888184
yt : [INFO     ] 2024-12-06 11:51:34,874 Parameters: omega_radiation           = 0.0
yt : [INFO     ] 2024-12-06 11:51:34,874 Parameters: hubble_constant           = 0.703000030517578
yt : [WARNING  ] 2024-12-06 11:51:34,878 Detected 2 extra particle fields assuming kind `double`. Consider using the `extra_particle_fields` keyword argument if you have unexpected behavior.
yt : [WARNING  ] 2024-12-06 11:51:34,884 This cooling file format is no longer supported. Cooling field loading skipped.
yt : [INFO     ] 2024-12-06 11:51:35,603 Identified    16/   16 intersecting domains (   16 through hilbert key indexing)
yt : [INFO     ] 2024-12-06 11:51:36,285 max value is 7.98220e-24 at 0.7520217895507812 0.7684097290039062 0.7971725463867188
yt : [INFO     ] 2024-12-06 11:51:36,302 Identified     2/   16 intersecting domains (    2 through hilbert key indexing)
yt : [INFO     ] 2024-12-06 11:51:36,337 xlim = -0.000161 0.000161
yt : [INFO     ] 2024-12-06 11:51:36,337 ylim = -0.000161 0.000161
yt : [INFO     ] 2024-12-06 11:51:36,337 zlim = -0.883346 0.883346
yt : [INFO     ] 2024-12-06 11:51:36,342 Making a fixed resolution buffer of (('gas', 'density')) 800 by 800
yt : [INFO     ] 2024-12-06 11:51:36,419 Identified    16/   16 intersecting domains (   16 through hilbert key indexing)
Traceback (most recent call last):
  File "/home/cphyc/Documents/prog/yt/debug/debug_offaxis2.py", line 12, in <module>
    p = yt.ProjectionPlot(ds, n, ("gas", "density"), center=sp.center, width=(20, "kpc"), weight_field=("gas", "dx"))
  File "/home/cphyc/Documents/prog/yt/yt/visualization/plot_window.py", line 2545, in __init__
    PWViewerMPL.__init__(
    ~~~~~~~~~~~~~~~~~~~~^
        self,
        ^^^^^
    ...<7 lines>...
        buff_size=buff_size,
        ^^^^^^^^^^^^^^^^^^^^
    )
    ^
  File "/home/cphyc/Documents/prog/yt/yt/visualization/plot_window.py", line 880, in __init__
    PlotWindow.__init__(self, *args, **kwargs)
    ~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/cphyc/Documents/prog/yt/yt/visualization/plot_window.py", line 268, in __init__
    self._setup_plots()
    ~~~~~~~~~~~~~~~~~^^
  File "/home/cphyc/Documents/prog/yt/yt/visualization/plot_window.py", line 1081, in _setup_plots
    image = self.frb.get_image(f)
  File "/home/cphyc/Documents/prog/yt/yt/visualization/fixed_resolution.py", line 213, in get_image
    self._generate_image_and_mask(key)
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^
  File "/home/cphyc/Documents/prog/yt/yt/visualization/fixed_resolution.py", line 643, in _generate_image_and_mask
    buff = off_axis_projection(
        dd.dd,
    ...<11 lines>...
        method=dd.method,
    )
  File "/home/cphyc/Documents/prog/yt/yt/visualization/volume_rendering/off_axis_projection.py", line 475, in off_axis_projection
    add_cells_to_image_offaxis(
    ~~~~~~~~~~~~~~~~~~~~~~~~~~^
        Xp=xyz,
        ^^^^^^^
    ...<7 lines>...
        Ny=resolution[1],
        ^^^^^^^^^^^^^^^^^
    )
    ^
  File "yt/utilities/lib/image_utilities.pyx", line 202, in yt.utilities.lib.image_utilities.add_cells_to_image_offaxis
    cdef np.ndarray[np.float64_t, ndim=3] stamp_arr = np.zeros((max_depth, Nsx, Nsy), dtype=float)
numpy._core._exceptions._ArrayMemoryError: Unable to allocate 2.37 TiB for an array with shape (18, 134527, 134527) and data type float64

Expected outcome

No error.

Version Information

• Operating System:
• Python Version:
• yt version: 4.4+
• Other Libraries (if applicable):

[matthewturk]

Any ideas how we might address it?

[cphyc]

Yes. At present, what the code does is:

1. Compute how many pixels $N$ the largest cell covers,
2. Compute once what a cell of this size seen from the current angle looks like, and store the result in a buffer of size $N^2$.
3. Do the same for cells that are half this big until the cell apparent size if <1 pixel.
4. Iterate over all cells, and use the precomputed buffer as a kernel for each.

Overall, the buffer is a 3D array of size $\mathcal{O}(N^2\log_2 N)$. The issue is that, if a cell is large compared to the field of view, the buffer may take a huge amount of space.

Possible solutions include:

• put an upper limit of how large the buffer can be, e.g. $N\lessapprox 100\mathrm{pixel}$ and interpolate for cells that would be larger than this,
• directly compute ray-cube intersection for large cells, and resort to the aforementioned approach for smaller ones,
• since the buffer is actually a piece-wise 2D linear function[1], we could also represent it analytically rather than using a 3D array.

[1]: The value only depends on the projected-plane position of each of the 8 nodes making a cell.