Repository: yt-project/unyt

Issue:
hdf5 reading for unyt.quantity.
Hi I believe that there is a bug whenever I try to read a `unyt.unyt_quantity ` that was stored in a hdf5 file using `write_hdf5`.  Basically I think it treats the quantity as a non-zero dimensional array, and uses ellipses to get the data.  I have tried to use `unyt.unyt_quantity.from_hdf5()` and `unyt.unyt_array.from_hdf5()`. 


```
----> 1 unyt.unyt_quantity.from_hdf5("...")
File (.../python3.9/site-packages/unyt/array.py:1557), in unyt_array.from_hdf5(cls, filename, dataset_name, group_name)
   1555     g = f
  1555     g = f
   1556 dataset = g[dataset_name]
-> 1557 data = dataset[:]
   1558 units = dataset.attrs.get("units", "")
   1559 unit_lut = default_unit_symbol_lut.copy()

File (objects.pyx:54), in h5py._objects.with_phil.wrapper()

File (objects.pyx:55), in h5py._objects.with_phil.wrapper()

File(..../h5py/_hl/dataset.py:817), in Dataset.__getitem__(self, args, new_dtype)
    815 if self.shape == ():
    816     fspace = self.id.get_space()
--> 817     selection = sel2.select_read(fspace, args)
    818     if selection.mshape is None:
    819         arr = numpy.zeros((), dtype=new_dtype)

File (../h5py/_hl/selections2.py:101), in select_read(fspace, args)
     96 """ Top-level dispatch function for reading.
     97 
     98 At the moment, only supports reading from scalar datasets.
     99 """
    100 if fspace.shape == ():
--> 101     return ScalarReadSelection(fspace, args)
    103 raise NotImplementedError()

File [.../h5py/_hl/selections2.py:86), in ScalarReadSelection.__init__(self, fspace, args)
     84     self.mshape = ()
     85 else:
---> 86     raise ValueError("Illegal slicing argument for scalar dataspace")
     88 self.mspace = h5s.create(h5s.SCALAR)
     89 self.fspace = fspace

ValueError: Illegal slicing argument for scalar dataspace
```

Hints from issue comments:
@neutrinoceros: Can you provide a self-contained reproducer script ?@aaramburo: Hi yes, here is a simple script that for me reproduces the issue

```python
import unyt

variable_quantity = 5*unyt.cm
variable_quantity.write_hdf5("example.hdf5",dataset_name="variable", group_name="example_quantity")
in_quantity = unyt.unyt_quantity.from_hdf5("example.hdf5",dataset_name="variable", group_name="example_quantity")
```@neutrinoceros: Thank you for reporting this. I opened a fix at #560

Solution patch:
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 30a4e2fd..5d6ecf0a 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -22,6 +22,7 @@ jobs:
           - '3.10'
           - '3.11'
           - '3.12'
+          - '3.13'
         # Test all on ubuntu, test ends on macos and windows
         include:
           - os: macos-13
@@ -30,9 +31,9 @@ jobs:
           - os: windows-latest
             python-version: '3.9'
           - os: macos-latest
-            python-version: '3.12'
+            python-version: '3.13'
           - os: windows-latest
-            python-version: '3.12'
+            python-version: '3.13'
 
 
     steps:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 13800da8..243e0ae1 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -29,7 +29,7 @@ repos:
             additional_dependencies: [black==24.10.0]
 
     - repo: https://github.com/astral-sh/ruff-pre-commit
-      rev: v0.9.4
+      rev: v0.9.9
       hooks:
           - id: ruff-format
           - id: ruff
diff --git a/tox.ini b/tox.ini
index 4d2010d0..65895cda 100644
--- a/tox.ini
+++ b/tox.ini
@@ -1,5 +1,5 @@
 [tox]
-envlist = py39-docs,begin,py39-dependencies,py39-versions,py{39,310,311,312},py39-unyt-module-test-function,end
+envlist = py39-docs,begin,py39-dependencies,py39-versions,py{39,310,311,312,313},py39-unyt-module-test-function,end
 isolated_build = True
 
 [gh-actions]
@@ -8,6 +8,7 @@ python =
     3.10: py310
     3.11: py311
     3.12: py312
+    3.13: py313
 
 [testenv]
 package = wheel
diff --git a/unyt/_array_functions.py b/unyt/_array_functions.py
index 56283a84..5a85fd6c 100644
--- a/unyt/_array_functions.py
+++ b/unyt/_array_functions.py
@@ -488,72 +488,72 @@ def block(arrays):
 
 
 @implements(np.fft.fft)
-def ftt_fft(a, *args, **kwargs):
+def fft_fft(a, *args, **kwargs):
     return np.fft.fft._implementation(np.asarray(a), *args, **kwargs) * a.units
 
 
 @implements(np.fft.fft2)
-def ftt_fft2(a, *args, **kwargs):
+def fft_fft2(a, *args, **kwargs):
     return np.fft.fft2._implementation(np.asarray(a), *args, **kwargs) * a.units
 
 
 @implements(np.fft.fftn)
-def ftt_fftn(a, *args, **kwargs):
+def fft_fftn(a, *args, **kwargs):
     return np.fft.fftn._implementation(np.asarray(a), *args, **kwargs) * a.units
 
 
 @implements(np.fft.hfft)
-def ftt_hfft(a, *args, **kwargs):
+def fft_hfft(a, *args, **kwargs):
     return np.fft.hfft._implementation(np.asarray(a), *args, **kwargs) * a.units
 
 
 @implements(np.fft.rfft)
-def ftt_rfft(a, *args, **kwargs):
+def fft_rfft(a, *args, **kwargs):
     return np.fft.rfft._implementation(np.asarray(a), *args, **kwargs) * a.units
 
 
 @implements(np.fft.rfft2)
-def ftt_rfft2(a, *args, **kwargs):
+def fft_rfft2(a, *args, **kwargs):
     return np.fft.rfft2._implementation(np.asarray(a), *args, **kwargs) * a.units
 
 
 @implements(np.fft.rfftn)
-def ftt_rfftn(a, *args, **kwargs):
+def fft_rfftn(a, *args, **kwargs):
     return np.fft.rfftn._implementation(np.asarray(a), *args, **kwargs) * a.units
 
 
 @implements(np.fft.ifft)
-def ftt_ifft(a, *args, **kwargs):
+def fft_ifft(a, *args, **kwargs):
     return np.fft.ifft._implementation(np.asarray(a), *args, **kwargs) * a.units
 
 
 @implements(np.fft.ifft2)
-def ftt_ifft2(a, *args, **kwargs):
+def fft_ifft2(a, *args, **kwargs):
     return np.fft.ifft2._implementation(np.asarray(a), *args, **kwargs) * a.units
 
 
 @implements(np.fft.ifftn)
-def ftt_ifftn(a, *args, **kwargs):
+def fft_ifftn(a, *args, **kwargs):
     return np.fft.ifftn._implementation(np.asarray(a), *args, **kwargs) * a.units
 
 
 @implements(np.fft.ihfft)
-def ftt_ihfft(a, *args, **kwargs):
+def fft_ihfft(a, *args, **kwargs):
     return np.fft.ihfft._implementation(np.asarray(a), *args, **kwargs) * a.units
 
 
 @implements(np.fft.irfft)
-def ftt_irfft(a, *args, **kwargs):
+def fft_irfft(a, *args, **kwargs):
     return np.fft.irfft._implementation(np.asarray(a), *args, **kwargs) * a.units
 
 
 @implements(np.fft.irfft2)
-def ftt_irfft2(a, *args, **kwargs):
+def fft_irfft2(a, *args, **kwargs):
     return np.fft.irfft2._implementation(np.asarray(a), *args, **kwargs) * a.units
 
 
 @implements(np.fft.irfftn)
-def ftt_irfftn(a, *args, **kwargs):
+def fft_irfftn(a, *args, **kwargs):
     return np.fft.irfftn._implementation(np.asarray(a), *args, **kwargs) * a.units
 
 
diff --git a/unyt/array.py b/unyt/array.py
index 8970312e..87341db4 100644
--- a/unyt/array.py
+++ b/unyt/array.py
@@ -1562,7 +1562,7 @@ def from_hdf5(cls, filename, dataset_name=None, group_name=None):
         else:
             g = f
         dataset = g[dataset_name]
-        data = dataset[:]
+        data = dataset[...]
         units = dataset.attrs.get("units", "")
         unit_lut = default_unit_symbol_lut.copy()
         unit_lut_load = pickle.loads(dataset.attrs["unit_registry"].tobytes())
@@ -1772,8 +1772,8 @@ def __pow__(self, p, mod=None, /):
         Power function
         """
         # see https://github.com/yt-project/unyt/issues/203
-        if p == 0.0:
-            ret = self.ua
+        if np.isscalar(p) and p == 0.0:
+            ret = self.unit_array
             ret.units = Unit("dimensionless")
             return ret
         else:
@@ -1862,17 +1862,32 @@ def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
                 u1 = Unit(registry=getattr(u0, "registry", None))
             elif ufunc is power:
                 u1 = inp1
-                if inp0.shape != () and inp1.shape != ():
-                    raise UnitOperationError(ufunc, u0, u1)
-                if isinstance(u1, unyt_array):
-                    if u1.units.is_dimensionless:
-                        pass
-                    else:
+                if inp0.shape == () or inp1.shape == ():
+                    if isinstance(u1, unyt_array) and not u1.units.is_dimensionless:
                         raise UnitOperationError(ufunc, u0, u1.units)
-                if u1.shape == ():
-                    u1 = float(u1)
+                    if u1.shape == ():
+                        u1 = float(u1)
+                    else:
+                        u1 = 1.0
+                elif inp0.shape == inp1.shape:
+                    if isinstance(u1, unyt_array) and not u1.units.is_dimensionless:
+                        raise UnitOperationError(ufunc, u0, getattr(u1, "units", None))
+
+                    if (
+                        (isinstance(u0, Unit) and not u0.is_dimensionless)
+                        or isinstance(u0, unyt_array)
+                        and not u0.units.is_dimensionless
+                    ):
+                        # u0 has units
+                        if np.ptp(u1) != 0:
+                            raise UnitOperationError(
+                                ufunc, u0, getattr(u1, "units", None)
+                            )
+
+                    first_element_slice = (0,) * u1.ndim
+                    u1 = float(u1[first_element_slice])
                 else:
-                    u1 = 1.0
+                    raise UnitOperationError(ufunc, u0, u1)
             unit_operator = self._ufunc_registry[ufunc]
 
             if (


Test patch:
diff --git a/unyt/tests/test_unyt_array.py b/unyt/tests/test_unyt_array.py
index de0d1a3c..07dd097c 100644
--- a/unyt/tests/test_unyt_array.py
+++ b/unyt/tests/test_unyt_array.py
@@ -528,6 +528,44 @@ def test_power():
     assert_array_equal_units(cm_arr**0, unyt_quantity(1.0, "dimensionless"))
 
 
+@pytest.mark.parametrize(
+    "p",
+    [
+        pytest.param(unyt_array(0), id="scalar power 0"),
+        pytest.param(unyt_array(1), id="scalar power 1"),
+        pytest.param(unyt_array([0, 0, 0, 0]), id="uniform power 0"),
+        pytest.param(unyt_array([1, 1, 1, 1]), id="uniform power 1"),
+        pytest.param(unyt_array([0, 1, 2, 3]), id="non-uniform power"),
+    ],
+)
+def test_dimensionless_array_power(p):
+    # see https://github.com/yt-project/unyt/issues/522
+    a = unyt_array([1, 2, 3, 4])
+    expected = a.value**p.value
+    assert_array_equal(expected, (a**p).value)
+
+
+@pytest.mark.parametrize(
+    "a, p",
+    [
+        pytest.param(
+            unyt_array([0, 1, 2, 3], "cm"),
+            unyt_array([0, 1, 2, 3]),
+            id="non-uniform power with arr non dimensionless",
+        ),
+        pytest.param(
+            unyt_array([0, 1, 2, 3]),
+            unyt_array([0, 0, 0, 0], "cm"),
+            id="non dimensionless power",
+        ),
+    ],
+)
+def test_invalid_array_power(a, p):
+    # see https://github.com/yt-project/unyt/issues/522
+    with pytest.raises(UnitOperationError):
+        a**p
+
+
 def test_comparisons():
     """
     Test numpy ufunc comparison operators for unit consistency.
@@ -1587,6 +1625,25 @@ def test_h5_io():
     shutil.rmtree(tmpdir)
 
 
+def test_h5_quantity_roundtrip(tmp_path):
+    # regression test for https://github.com/yt-project/unyt/issues/559
+    pytest.importorskip("h5py")
+
+    q1 = 5 * cm
+    savefile = tmp_path / "savefile.hdf5"
+    q1.write_hdf5(
+        savefile,
+        dataset_name="test_ds",
+        group_name="test_group",
+    )
+    q2 = unyt_quantity.from_hdf5(
+        savefile,
+        dataset_name="test_ds",
+        group_name="test_group",
+    )
+    assert_array_equal_units(q2, q1)
+
+
 def test_equivalencies():
     import unyt as u