Repository: yt-project/unyt

Issue:
Subtracting temperatures gives degC rather than delta_degC 
* unyt version:
'v2.8.0'
* Python version:
3.8.12
* Operating System:
Mac

### Description

Working with change in temperature. I would expect a difference between two temperatures to return a delta_degC rather than a degC. 

### What I Did

Actual
```
>>> import unyt as u
>>> delta = 2 * u.degC - 1 * u.degC 
>>> delta 
unyt_quantity(1, 'degC')
>>> delta.to('K')
unyt_quantity(274.15, 'K')
```

Expected

```
>>> import unyt as u
>>> delta = 2 * u.degC - 1 * u.degC 
>>> delta 
unyt_quantity(1, 'delta_degC')
>>> delta.to('K')
unyt_quantity(1, 'K')
```

Discussion:
ngoldbaum: Let me know if you'd like a hand with this.
neutrinoceros: Sure ! I'm actually having a hard time understanding how my changes affect tests `unyt/tests/test_unyt_array.py`
ngoldbaum: > I'm actually having a hard time understanding how my changes affect tests `unyt/tests/test_unyt_array.py`

See [d708386](https://github.com/yt-project/unyt/pull/413/commits/d708386c3bd643ba2bff1429a16f9fa926aa9ccc)
neutrinoceros: awesome, well played ! The one remaining failure (in docs) touches a case that I'm not sure how to treat. Should `1*delta_degC - 1*degC` remain a valid operation ?
ngoldbaum: It looks like Pint returns data in units of `degC`, so I think we should probably follow that. I guess one way to justify it is:

```
1*delta_degC - 1*degC = -1 * (1*degC - 1*delta_degC)
```

So it makes sense for both operations to return results in the same units.

It also means we're not raising a new error without a deprecation, which is nice.

Patch:
diff --git a/tox.ini b/tox.ini
index 8d42841..f220e14 100644
--- a/tox.ini
+++ b/tox.ini
@@ -21,7 +21,7 @@ deps =
     pint
     astropy
     coverage[toml]>=5.0
-    packaging>=20.9,
+    packaging>=20.9
     pytest-cov
     pytest-doctestplus
     matplotlib!=3.5.0
diff --git a/unyt/array.py b/unyt/array.py
index 517dae0..ef06fb4 100644
--- a/unyt/array.py
+++ b/unyt/array.py
@@ -127,6 +127,7 @@ from unyt.unit_registry import (
     _sanitize_unit_system,
     default_unit_registry,
 )
+from unyt.unit_symbols import delta_degC, delta_degF
 
 from ._deprecation import warn_deprecated
 
@@ -194,6 +195,42 @@ def _preserve_units(unit1, unit2=None):
     return 1, unit1
 
 
+@lru_cache(maxsize=128, typed=False)
+def _difference_units(unit1, unit2=None):
+    if unit1.dimensions is not temperature:
+        return _preserve_units(unit1, unit2)
+
+    s1 = repr(unit1)
+    if unit2 is not None and unit2 != unit1:
+        s2 = repr(unit2)
+        if s1 in s2 and s2.startswith("delta_"):
+            return 1, unit1
+        elif s2 in s1 and s1.startswith("delta_"):
+            return 1, unit2
+        else:
+            raise InvalidUnitOperation(
+                "Quantities with units of Fahrenheit or Celsius "
+                "cannot be multiplied, divided, subtracted or "
+                "added with data that has different units."
+            )
+
+    if unit1.base_offset == 0.0:
+        return 1, unit1
+
+    if s1 == "degF":
+        return 1, delta_degF
+    elif s1 == "degC":
+        return 1, delta_degC
+    else:
+        # This is supposed to be unreachable
+        raise RuntimeError(
+            "Could not determine difference temperature units "
+            f"in operation ({unit1} - {unit2}).\n"
+            "If you see this error please file an issue at "
+            "https://github.com/yt-project/unyt/issues/new"
+        )
+
+
 @lru_cache(maxsize=128, typed=False)
 def _power_unit(unit, power):
     return 1, unit**power
@@ -444,7 +481,7 @@ class unyt_array(np.ndarray):
 
     _ufunc_registry = {
         add: _preserve_units,
-        subtract: _preserve_units,
+        subtract: _difference_units,
         multiply: _multiply_units,
         divide: _divide_units,
         logaddexp: _return_without_unit,
@@ -1834,7 +1871,12 @@ class unyt_array(np.ndarray):
             ):
                 raise UnitOperationError(ufunc, u0, u1)
 
-            if unit_operator in (_preserve_units, _comparison_unit, _arctan2_unit):
+            if unit_operator in (
+                _preserve_units,
+                _comparison_unit,
+                _arctan2_unit,
+                _difference_units,
+            ):
                 # check "is" equality first for speed
                 if u0 is not u1 and u0 != u1:
                     # we allow adding, multiplying, comparisons with


Test patch:
diff --git a/unyt/tests/test_array_functions.py b/unyt/tests/test_array_functions.py
index 95ca32f..82ae07f 100644
--- a/unyt/tests/test_array_functions.py
+++ b/unyt/tests/test_array_functions.py
@@ -4,14 +4,13 @@ import re
 import numpy as np
 import pytest
 
-from unyt import A, K, cm, degC, degF, delta_degC, g, km, rad, s
+from unyt import A, K, cm, delta_degC, g, km, rad, s
 from unyt._array_functions import (
     _HANDLED_FUNCTIONS as HANDLED_FUNCTIONS,
     _UNSUPPORTED_FUNCTIONS as UNSUPPORTED_FUNCTIONS,
 )
 from unyt.array import unyt_array, unyt_quantity
 from unyt.exceptions import (
-    InvalidUnitOperation,
     UnitConversionError,
     UnitInconsistencyError,
     UnytError,
@@ -1244,28 +1243,6 @@ def test_deltas(func, input_units, output_units):
     assert res.units == output_units
 
 
-@pytest.mark.parametrize(
-    "func",
-    [
-        np.diff,
-        np.ediff1d,
-        np.gradient,
-        np.ptp,
-    ],
-)
-@pytest.mark.parametrize("input_units", [degC, degF])
-def test_invalid_deltas(func, input_units):
-    x = np.arange(0, 4) * input_units
-    with pytest.raises(
-        InvalidUnitOperation,
-        match=re.escape(
-            "Quantities with units of Fahrenheit or Celsius "
-            "cannot be multiplied, divided, subtracted or added."
-        ),
-    ):
-        func(x)
-
-
 @pytest.mark.parametrize(
     "func",
     [
diff --git a/unyt/tests/test_unyt_array.py b/unyt/tests/test_unyt_array.py
index 76c9a93..e61db46 100644
--- a/unyt/tests/test_unyt_array.py
+++ b/unyt/tests/test_unyt_array.py
@@ -2493,6 +2493,26 @@ def test_delta_degF():
     assert 2 * delta_degF == unyt_quantity(2, "delta_degF")
 
 
+@pytest.mark.parametrize(
+    ("u0", "u1", "uout"),
+    [
+        (K, K, K),
+        (R, R, R),
+        (degC, degC, delta_degC),
+        (degF, degF, delta_degF),
+        (degC, delta_degC, degC),
+        (delta_degC, degC, degC),
+        (degF, delta_degF, degF),
+        (delta_degF, degF, degF),
+    ],
+)
+def test_delta_temperature_diff(u0, u1, uout):
+    # using repr comparison because
+    # 1) we don't care that Unit instances might not be the same
+    # 2) some temperature units will compare as equal even when they are not identical (e.g. K and delta_degC)
+    assert repr((2 * u0 - 1 * u1).units) == repr(uout)
+
+
 def test_mil():
     assert_allclose_units(unyt_quantity(1, "mil"), unyt_quantity(0.001, "inch"))