yt-project / unyt

Handle, manipulate, and convert data with units in Python
https://unyt.readthedocs.io
BSD 3-Clause "New" or "Revised" License
364 stars 48 forks source link

An error occurs when the exponent is an array #522

Open xshaokun opened 6 days ago

xshaokun commented 6 days ago

Description

In the power operation of unyt arrays, when the exponent is also an array, it leads to an error.

In [1]: from unyt import unyt_array

In [2]: a = unyt_array([1,2,3,4])

In [3]: b = unyt_array([1,2,3,4])

In [4]: a ** b
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[4], line 1
----> 1 a ** b

File ~/Apps/anaconda3/lib/python3.9/site-packages/unyt/array.py:1745, in unyt_array.__pow__(self, p, mod)
   1737 def __pow__(self, p, mod=None, /):
   1738     """
   1739     Power function, over-rides the ufunc as
   1740     ``numpy`` passes the ``_ones_like`` ufunc for
   1741     powers of zero (in some cases), which leads to an
   1742     incorrect unit calculation.
   1743     """
-> 1745     if p == 0.0:
   1746         ret = self.ua
   1747         ret.units = Unit("dimensionless")

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

While for numpy arrays, it behaviors like

In [5]: import numpy as np

In [6]: a = np.array([1,2,3,4])

In [7]: b = np.array([1,2,3,4])

In [8]: a ** b
Out[8]: array([  1,   4,  27, 256])

I now use to_ndarray() to make it work but lose the units.

neutrinoceros commented 6 days ago

Thanks for reporting this, it's clearly a defect ! It might sufficient to replace p == 0.0 with np.all(p == 0.0) and add a couple test cases. We should also make sure that a TypeError is raised if p is another (non dimensionless) unyt array. Are you interested in giving it a try ? Otherwise I should have time to do it myself later this week.

xshaokun commented 5 days ago

Thank you for the response. From the perspective of dimensions, p as an exponent here should not be a quantity with dimensions, so I don't think we need to worry about p being non-dimensionless. Since I will be quite busy for a while and may not be able to fix this bug, I appreciate you would make it :)

neutrinoceros commented 5 days ago

Hum, maybe I was a bit quick in my previous answer: the only cases where it makes sense to support raising a unyt_array a to an array power p is when p is a 1-element array or has a broadcastable shape with a and a unique value. The reason is that the resulting array can only have a uniform unit, which is only possible if all elements are raised to the same power... so, even though the specific example you illustrated the issue with doesn't make sense to me (but the error should be clearer), I'll add support for the few cases that do make sense.

xshaokun commented 4 days ago

Oh, in some cases, when the base is a dimensionless number, the example I gave should make sense. For instance, I’m currently trying to calculate synchrotron radiation, where the spectral indices vary across different regions. The following formula shows that the base in the final part is a dimensionless quantity, where the p is the power-law index of non-thermal electrons.

image

(from Rybicki and Lightman 1986)

neutrinoceros commented 4 days ago

Yes, if the power itself is dimensionless, it should work. See #524 and please let me know if this helps !