Segmentation fault with u.Msun

[kyleaoman]

• unyt version: https://github.com/kyleaoman/unyt/tree/hist_bugfix commit 13b6da0
• Python version: 3.12
• Operating System: Linux (Fedora 40)

Description

Trying to histogram a unyt array initialized with unyt_array(arr, units) in density mode when units are u.Msun crashes with a segmentation fault... but not when units are u.solMass, and not when the array is initialized as np.array(...) * u.Msun.

What I Did

This is rearing its head in the updated handling of the histogramdd argument density but is pretty clearly more fundamental than anything I've directly introduced. Here's as much as I've managed to isolate it so far:

import numpy as np
import unyt as u

sample = [
    u.unyt_array(np.arange(3), u.Msun),
]
np.histogramdd(
    sample,
    density=True,
)

That exits with a segmentation fault.

This however exits cleanly:

import numpy as np
import unyt as u

sample = [
    u.unyt_array(np.arange(3), u.solMass),
]
np.histogramdd(
    sample,
    density=True,
)

So does this:

import numpy as np
import unyt as u

sample = [
    np.arange(3) * u.Msun,
]
np.histogramdd(
    sample,
    density=True,
)

[neutrinoceros]

Can you share the exact version of Python you're using (down to the patch level), your provider (Python.org, conda, something else ?) and the version of numpy you are running with ? unyt itself is pure Python so I highly doubt that the memory access error is on our side.

[neutrinoceros]

FWIW, running your first example on my system produces a RecursionError, not a segfault. This is very likely showing a problem in unyt itself.

[neutrinoceros]

Note that the key difference betwee unyt.msun and unyt.solMass is that the former is a quantity (whose value is stretching the limits of floating point arithmetic) while the latter is a unit.

[neutrinoceros]

I'm able to avoid this RecursionError with the following patch

diff --git a/unyt/_array_functions.py b/unyt/_array_functions.py
index 62aca69..a752fd8 100644
--- a/unyt/_array_functions.py
+++ b/unyt/_array_functions.py
@@ -288,7 +288,8 @@ def _histogramdd(
     if density:
         for s in sample:
             counts /= getattr(s, "units", 1)
-    counts *= getattr(weights, "units", 1)
+    if weights is not None and hasattr(weights, "units"):
+        counts *= weights.units
     return counts, tuple(_bin * getattr(s, "units", 1) for _bin, s in zip(bins, sample))

(it also seems to work when I pass weights, regardless if it has units or not). I'm going to propose patching this case and similar ones from your previous PR.

[kyleaoman]

That seems like it should be equivalent, wonder what's triggering the recursion... also wonder if it's a deeper issue that might crop up elsewhere, or if this fix is enough.

(Also just to mention - when I first encountered this in a pytest run it came with a partial stack trace before the segfault, that indeed looked like something recursing a lot.)

[neutrinoceros]

That seems like it should be equivalent, wonder what's triggering the recursion...

multipliying an array by a scalar 1 calls in some ufunc/array func logic, but at this point we're already deep in unyt_array.__array_function__. I think that's where the issue is coming from, though I must admit I didn't look much closer after I found a working solution. PR incoming.