Repository: yt-project/unyt (python)

Issue:
Is it possible to get Hz * s to reduce?
* unyt version: 1.0.6
* Python version: 3.6.6
* Operating System: MacOS 10.13.6

### Description

I multiplied Hz and seconds and got a Hz * s quantity.
I expected a result without dimension, as in Km / Km.

### What I Did

```
In [281]: (1*u.Hz) * (1*u.s)
Out[281]: unyt_quantity(1.0, 'Hz*s')

In [282]: (1*u.km) / (1*u.km)
Out[282]: unyt_quantity(1.0, '(dimensionless)')

In [286]: ((1*u.Hz) * (1*u.s)).to(u.dimensionless)
Out[286]: unyt_quantity(1.0, '(dimensionless)')
```


Pull request description:
add automatic unit simplification. closes #49
This adds automatic simplification of units in the form of a new `Unit.simplify()` method. It also makes it so that products and quotients of arrays with units will now automatically be simplified.

Caveats:

* I haven't done any performance benchmarking. In principle calling `simplify()` in the `__array_ufunc__` machinery is "free" because it's behind an LRU cache but I need to actually verify that. There's also an extra multiplication on the array that happens if there's any unit cancellation. In practice though most people would be "manually" doing that simplification and therefore incurring the cost of the extra multiplication anyway. All that said I still need to actually verify that the performance impact is at least acceptable.
* My tests are not *super* extensive and this is an invasive change. It's entirely possible there are untested corner cases.

Patch:
diff --git a/docs/usage.rst b/docs/usage.rst
index c00233d..effa640 100644
--- a/docs/usage.rst
+++ b/docs/usage.rst
@@ -83,16 +83,15 @@ end, days::
 
   >>> period = 2*pi*(semimajor_axis**3/(G*Mjup))**0.5
   >>> period
-  unyt_array([2.64196224e-12, 5.30291672e-12, 1.06837107e-11,
-              2.49212918e-11], 'AU**(3/2)*s/m**(3/2)')
+  unyt_array([ 152867.34547843,  306833.60667034,  618173.2944961 ,
+                1441978.11592457], 's')
   >>> period.to('d')
   unyt_array([ 1.76929798,  3.55131489,  7.1547835 , 16.68956153], 'd')
 
 Note that we haven't added any conversion factors between different units,
-that's all handled internally by :mod:`unyt`. Also note how the intermediate
-result ended up with complicated, ugly units, but the :meth:`unyt_array.to
-<unyt.array.unyt_array.to>` method was able to automagically handle the
-conversion to days.
+that's all handled internally by :mod:`unyt`. Also note how the
+:meth:`unyt_array.to <unyt.array.unyt_array.to>` method was able to
+automagically handle the conversion from seconds to days.
 
 It's also worth emphasizing that :mod:`unyt` represents powers using standard
 python syntax. This means you must use `**` and not `^`, even when writing a
@@ -229,7 +228,9 @@ controlled identically to numpy arrays, using ``numpy.setprintoptions``:
   [1.1235] km
   >>> np.set_printoptions(precision=8)
 
-Print a :math:`\rm{\LaTeX}` representation of a set of units using the :meth:`unyt.unit_object.Unit.latex_representation` function or :attr:`unyt.unit_object.Unit.latex_repr` attribute:
+Print a :math:`\rm{\LaTeX}` representation of a set of units using the
+:meth:`unyt.unit_object.Unit.latex_representation` function or
+:attr:`unyt.unit_object.Unit.latex_repr` attribute:
 
   >>> from unyt import g, cm
   >>> (g/cm**3).units.latex_representation()
@@ -237,6 +238,43 @@ Print a :math:`\rm{\LaTeX}` representation of a set of units using the :meth:`un
   >>> (g/cm**3).units.latex_repr
   '\\frac{\\rm{g}}{\\rm{cm}^{3}}'
 
+
+Simplifying Units
+-----------------
+
+Unit expressions can often be simplified to cancel pairs of factors with
+compatible dimensions. For example, we can form a unit with dimensions of length
+by dividing a unit with dimensions of length squared by another unit with
+dimensions of length::
+
+  >>> from unyt import m, cm
+  >>> m**2/cm
+  m**2/cm
+
+The :class:`Unit <unyt.unit_object.Unit>` class has a :meth:`simplify()
+<unyt.unit_object.Unit.simplify>` method that we can call to create a new unit
+object to that includes the dimensionless ratio ``m/cm`` as a constant
+coefficient::
+
+  >>> (m**2/cm).simplify()
+  100*m
+
+This will also work for units that are the reciprocals of each other, for example:
+
+  >>> from unyt import s, Hz
+  >>> (s*Hz).simplify()
+  (dimensionless)
+
+Products and quotients of unit objects will not be simplified unless
+``simplify()`` is called explicitly. However, products and quotients of arrays
+and quantities will be simplified to make interactive work more intuitive::
+
+  >>> from unyt import erg, minute, hour
+  >>> power = [20, 40, 80] * erg / minute
+  >>> elapsed_time = 3*hour
+  >>> print(power*elapsed_time)
+  [ 3600.  7200. 14400.] erg
+
 Unit Conversions and Unit Systems
 +++++++++++++++++++++++++++++++++
 
diff --git a/unyt/array.py b/unyt/array.py
index 39cc7c6..661a363 100644
--- a/unyt/array.py
+++ b/unyt/array.py
@@ -166,52 +166,54 @@ def _return_arr(func):
 
 @lru_cache(maxsize=128, typed=False)
 def _sqrt_unit(unit):
-    return unit ** 0.5
+    return 1, unit ** 0.5
 
 
 @lru_cache(maxsize=128, typed=False)
 def _multiply_units(unit1, unit2):
-    return unit1 * unit2
+    ret = (unit1 * unit2).simplify()
+    return ret.as_coeff_unit()
 
 
 def _preserve_units(unit1, unit2=None):
-    return unit1
+    return 1, unit1
 
 
 @lru_cache(maxsize=128, typed=False)
 def _power_unit(unit, power):
-    return unit ** power
+    return 1, unit ** power
 
 
 @lru_cache(maxsize=128, typed=False)
 def _square_unit(unit):
-    return unit * unit
+    return 1, unit * unit
 
 
 @lru_cache(maxsize=128, typed=False)
 def _divide_units(unit1, unit2):
-    return unit1 / unit2
+    ret = (unit1 / unit2).simplify()
+    return ret.as_coeff_unit()
 
 
 @lru_cache(maxsize=128, typed=False)
 def _reciprocal_unit(unit):
-    return unit ** -1
+    return 1, unit ** -1
 
 
 def _passthrough_unit(unit, unit2=None):
-    return unit
+    return 1, unit
 
 
 def _return_without_unit(unit, unit2=None):
-    return None
+    return 1, None
 
 
 def _arctan2_unit(unit1, unit2):
-    return NULL_UNIT
+    return 1, NULL_UNIT
 
 
 def _comparison_unit(unit1, unit2=None):
-    return None
+    return 1, None
 
 
 def _invert_units(unit):
@@ -1092,9 +1094,9 @@ class unyt_array(np.ndarray):
         >>> (1.0*km).has_equivalent('spectral')
         True
         >>> print((1*km).to_equivalent('MHz', equivalence='spectral'))
-        0.29979245800000004 MHz
+        0.299792458 MHz
         >>> print((1*keV).to_equivalent('angstrom', equivalence='spectral'))
-        12.398419315219659 angstrom
+        12.39841931521966 angstrom
         """
         return self.units.has_equivalent(equivalence)
 
@@ -1629,6 +1631,7 @@ class unyt_array(np.ndarray):
             if ufunc in (multiply, divide) and method == "reduce":
                 # a reduction of a multiply or divide corresponds to
                 # a repeated product which we implement as an exponent
+                mul = 1
                 power_sign = POWER_SIGN_MAPPING[ufunc]
                 if "axis" in kwargs and kwargs["axis"] is not None:
                     unit = u ** (power_sign * inp.shape[kwargs["axis"]])
@@ -1636,7 +1639,7 @@ class unyt_array(np.ndarray):
                     unit = u ** (power_sign * inp.size)
             else:
                 # get unit of result
-                unit = self._ufunc_registry[ufunc](u)
+                mul, unit = self._ufunc_registry[ufunc](u)
             # use type(self) here so we can support user-defined
             # subclasses of unyt_array
             ret_class = type(self)
@@ -1707,7 +1710,7 @@ class unyt_array(np.ndarray):
                         )
                     inp1 = np.asarray(inp1) * conv
             # get the unit of the result
-            unit = unit_operator(u0, u1)
+            mul, unit = unit_operator(u0, u1)
             # actually evaluate the ufunc
             out_arr = func(
                 inp0.view(np.ndarray), inp1.view(np.ndarray), out=out_func, **kwargs
@@ -1759,13 +1762,19 @@ class unyt_array(np.ndarray):
             else:
                 out_arr = ret_class(out_arr, unit, bypass_validation=True)
         if out is not None:
+            if mul != 1:
+                multiply(out, mul, out=out)
+                if np.shares_memory(out_arr, out):
+                    mul = 1
             if isinstance(out, unyt_array):
                 try:
                     out.units = out_arr.units
                 except AttributeError:
                     # out_arr is an ndarray
                     out.units = Unit("", registry=self.units.registry)
-        return out_arr
+        if mul == 1:
+            return out_arr
+        return mul * out_arr
 
     def copy(self, order="C"):
         """
diff --git a/unyt/equivalencies.py b/unyt/equivalencies.py
index 802880b..ee6d54e 100644
--- a/unyt/equivalencies.py
+++ b/unyt/equivalencies.py
@@ -505,7 +505,8 @@ class EffectiveTemperatureEquivalence(Equivalence):
             T4 = np.true_divide(
                 x, pc.stefan_boltzmann_constant_mks, out=self._get_out(x)
             )
-            return np.power(T4, 0.25, out=self._get_out(x))
+            ret = np.power(T4, 0.25, out=self._get_out(x))
+            return ret
 
     def __str__(self):
         return "effective_temperature: flux <-> temperature"
diff --git a/unyt/unit_object.py b/unyt/unit_object.py
index 35126cb..2746629 100644
--- a/unyt/unit_object.py
+++ b/unyt/unit_object.py
@@ -14,6 +14,7 @@ A class that represents a unit symbol.
 
 
 import copy
+import itertools
 
 try:
     from functools import lru_cache
@@ -37,6 +38,8 @@ from sympy import (
     Basic,
     Rational,
     sqrt,
+    Mod,
+    floor,
 )
 from sympy.core.numbers import One
 from sympy import sympify, latex
@@ -757,6 +760,94 @@ class Unit(object):
         """
         return self.latex_repr
 
+    def as_coeff_unit(self):
+        """Factor the coefficient multiplying a unit
+
+        For units that are multiplied by a constant dimensionless
+        coefficient, returns a tuple containing the coefficient and
+        a new unit object for the unmultiplied unit.
+
+        Example
+        -------
+
+        >>> import unyt as u
+        >>> unit = (u.m**2/u.cm).simplify()
+        >>> unit
+        100*m
+        >>> unit.as_coeff_unit()
+        (100.0, m)
+        """
+        coeff, mul = self.expr.as_coeff_Mul()
+        coeff = float(coeff)
+        ret = Unit(
+            mul,
+            self.base_value / coeff,
+            self.base_offset,
+            self.dimensions,
+            self.registry,
+        )
+        return coeff, ret
+
+    def simplify(self):
+        """Return a new equivalent unit object with a simplified unit expression
+
+        >>> import unyt as u
+        >>> unit = (u.m**2/u.cm).simplify()
+        >>> unit
+        100*m
+        """
+        expr = self.expr
+        self.expr = _cancel_mul(expr, self.registry)
+        return self
+
+
+def _factor_pairs(expr):
+    factors = expr.as_ordered_factors()
+    expanded_factors = []
+    for f in factors:
+        if f.is_Number:
+            continue
+        base, exp = f.as_base_exp()
+        if exp.q != 1:
+            expanded_factors.append(base ** Mod(exp, 1))
+            exp = floor(exp)
+        if exp >= 0:
+            f = (base,) * exp
+        else:
+            f = (1 / base,) * abs(exp)
+        expanded_factors.extend(f)
+    return list(itertools.combinations(expanded_factors, 2))
+
+
+def _create_unit_from_factor(factor, registry):
+    base, exp = factor.as_base_exp()
+    f = registry[str(base)]
+    return Unit(base, f[0], f[2], f[1], registry, f[3]) ** exp
+
+
+def _cancel_mul(expr, registry):
+    pairs_to_consider = _factor_pairs(expr)
+    uncancelable_pairs = set()
+    while len(pairs_to_consider):
+        pair = pairs_to_consider.pop()
+        if pair in uncancelable_pairs:
+            continue
+        u1 = _create_unit_from_factor(pair[0], registry)
+        u2 = _create_unit_from_factor(pair[1], registry)
+        prod = u1 * u2
+        if prod.dimensions == 1:
+            expr = expr / pair[0]
+            expr = expr / pair[1]
+            value = prod.base_value
+            if value != 1:
+                if value.is_integer():
+                    value = int(value)
+                expr *= value
+        else:
+            uncancelable_pairs.add(pair)
+        pairs_to_consider = _factor_pairs(expr)
+    return expr
+
 
 #
 # Unit manipulation functions
@@ -1084,3 +1175,6 @@ def define_unit(
     if registry is default_unit_registry:
         u = Unit(symbol, registry=registry)
         setattr(unyt, symbol, u)
+
+
+NULL_UNIT = Unit()


Test patch:
diff --git a/unyt/tests/test_units.py b/unyt/tests/test_units.py
index 2f79b21..d796e64 100644
--- a/unyt/tests/test_units.py
+++ b/unyt/tests/test_units.py
@@ -653,3 +653,25 @@ def test_em_unit_base_equivalent():
 def test_symbol_lut_length():
     for v in default_unit_symbol_lut.values():
         assert len(v) == 5
+
+
+def test_simplify():
+    import unyt as u
+
+    answers = {
+        u.Hz * u.s: "dimensionless",
+        u.kg / u.g: "1000",
+        u.Hz * u.s * u.km: "km",
+        u.kHz * u.s: "1000",
+        u.kHz * u.s * u.km: "1000*km",
+        u.kHz * u.s ** 2: "1000*s",
+        u.kHz * u.s ** 2 * u.km: "1000*km*s",
+        u.Hz ** -1 * u.s: "s/Hz",
+        u.Hz ** -1 * u.s * u.km: "km*s/Hz",
+        u.Hz ** 1.5 * u.s ** 1.7: "sqrt(Hz)*s**(7/10)",
+        u.Hz ** 1.5 * u.s ** 1.7 * u.km: "sqrt(Hz)*km*s**(7/10)",
+        u.m ** 2 / u.cm ** 2: "10000",
+    }
+
+    for unit, answer in answers.items():
+        assert str(unit.simplify()) == answer
diff --git a/unyt/tests/test_unyt_array.py b/unyt/tests/test_unyt_array.py
index eb8f466..d59e224 100644
--- a/unyt/tests/test_unyt_array.py
+++ b/unyt/tests/test_unyt_array.py
@@ -1039,7 +1039,7 @@ def binary_ufunc_comparison(ufunc, a, b):
     ret = ufunc(a, b)
 
     if ufunc is np.multiply:
-        assert ret.units == a.units * b.units
+        assert ret.units == (a.units * b.units).simplify().as_coeff_unit()[1]
     elif ufunc in (np.divide, np.true_divide, np.arctan2):
         assert ret.units.dimensions == (a.units / b.units).dimensions
     elif ufunc in (
@@ -1715,11 +1715,12 @@ def test_equivalencies():
     assert_allclose_units(T, F.in_units("K", equivalence="effective_temperature"))
     T.convert_to_units("erg/s/cm**2", "effective_temperature")
     assert_allclose_units(T, F)
-    assert T.units == unyt_quantity(1.0, "erg/cm**2/s").units
+    assert T.units == u.Unit("erg/cm**2/s")
     assert F.units == u.W / u.m ** 2
+    assert_almost_equal(T.in_units("K", "effective_temperature").value, 1e4)
     T.convert_to_units("K", "effective_temperature")
     assert_almost_equal(T.value, 1e4)
-    assert T.units == u.K.units
+    assert T.units == u.K
 
     # to_value test