Refactor power derivative functions (#294)

- [x] Closes #293
- [x] Executed `pre-commit run --all-files` with no errors
- [x] The change is fully covered by automated unit tests
- [x] Documented in docs/ as appropriate
- [x] Added an entry to the CHANGES file

`uncertainties` allows calculating `x**y` where one or both of `x` and
`y` are `UFloat` objects. For this we need to calculate the `x` and `y`
derivatives of `x**y`. This PR refactors the functions that calculate
those derivatives to more explicitly and clearly handle non-typical
cases. Moving those functions to the module level also allows these
functions to be directly tested.

---------

Co-authored-by: Matt Newville <newville@cars.uchicago.edu>

Author: jagerber48

CHANGES.rst

@@ -12,6 +12,8 @@ Changes
    time. Now such a function can be imported but if the user attempts to
    execute it, a `NotImplementedError` is raised indicating that the
    function can't be used because `numpy` couldn't be imported.
+- Refactored the implementation for the calculation of the derivatives of the power
+   function and improves the corresponding testing.
 
 Adds:
 

doc/user_guide.rst

@@ -271,6 +271,64 @@ To check whether the uncertainty is NaN or Inf, use one of :func:`math.isnan`,
 .. index:: correlations; detailed example
 
 
+Power Function Behavior
+=======================
+
+The value of one :class:`UFloat` raised to the power of another can be calculated in two
+ways:
+
+>>> from uncertainties import umath
+>>>
+>>> x = ufloat(4.5, 0.2)
+>>> y = ufloat(3.4, 0.4)
+>>> print(x**y)
+(1.7+/-1.0)e+02
+>>> print(umath.pow(x, y))
+(1.7+/-1.0)e+02
+
+The function ``x**y`` is defined for all ``x != 0`` and for ``x == 0`` as long as
+``y > 0``.
+There is not a unique definition for ``0**0``, however python takes the convention for
+:class:`float` that ``0**0 == 1``.
+If the power operation is performed on an ``(x, y)`` pair for which ``x**y`` is
+undefined then an exception will be raised:
+
+>>> x = ufloat(0, 0.2)
+>>> y = ufloat(-3.4, 0.4)
+>>> print(x**y)
+Traceback (most recent call last):
+ ...
+ZeroDivisionError: 0.0 cannot be raised to a negative power
+
+On the domain where it is defined, ``x**y`` is always real for ``x >= 0``.
+For ``x < 0`` it is real for all integer values of ``y``.
+If ``x<0`` and ``y`` is not an integer then ``x**y`` has a non-zero imaginary component.
+The :mod:`uncertainties` module does not handle complex values:
+
+>>> x = ufloat(-4.5, 0.2)
+>>> y = ufloat(-3.4, 0.4)
+>>> print(x**y)
+Traceback (most recent call last):
+ ...
+ValueError: The uncertainties module does not handle complex results
+
+The ``x`` derivative is real anywhere ``x**y`` is real except along ``x==0`` for
+non-integer ``y``.
+At these points the ``x`` derivative would be complex so a NaN value is used:
+
+>>> x = ufloat(0, 0.2)
+>>> y=1.5
+>>> print((x**y).error_components())
+{0.0+/-0.2: nan}
+
+The ``y`` derivative is real anywhere ``x**y`` is real as long as ``x>=0``.
+For ``x < 0`` the ``y`` derivative is always complex valued so a NaN value is used:
+
+>>> x = -2
+>>> y = ufloat(1, 0.2)
+>>> print((x**y).error_components())
+{1.0+/-0.2: nan}
+
 Automatic correlations
 ======================
 

tests/test_power.py

@@ -3,11 +3,40 @@
 import pytest
 
 from uncertainties import ufloat
+from uncertainties.ops import pow_deriv_0, pow_deriv_1
 from uncertainties.umath_core import pow as umath_pow
 
 from helpers import nan_close
 
 
+pow_deriv_cases = [
+    (0.5, 2, 1.0, -0.17328679513998632),
+    (0.5, 1.5, 1.0606601717798214, -0.2450645358671368),
+    (0.5, 0, 0.0, -0.6931471805599453),
+    (0.5, -1.5, -8.485281374238571, -1.9605162869370945),
+    (0.5, -2, -16.0, -2.772588722239781),
+    (0, 2, 0, 0),
+    (0, 1.5, float("nan"), 0),
+    (0, 0, 0, float("nan")),
+    (0, -0.5, float("nan"), float("nan")),
+    (0, -2, float("nan"), float("nan")),
+    (-0.5, 2, -1.0, float("nan")),
+    (-0.5, 1.5, float("nan"), float("nan")),
+    (-0.5, 0, -0.0, float("nan")),
+    (-0.5, -1.5, float("nan"), float("nan")),
+    (-0.5, -2, 16.0, float("nan")),
+]
+
+
+@pytest.mark.parametrize("x, y, x_deriv_expected, y_deriv_expected", pow_deriv_cases)
+def test_pow_deriv_0(x, y, x_deriv_expected, y_deriv_expected):
+    x_deriv_actual = pow_deriv_0(x, y)
+    assert nan_close(x_deriv_actual, x_deriv_expected)
+
+    y_deriv_actual = pow_deriv_1(x, y)
+    assert nan_close(y_deriv_actual, y_deriv_expected)
+
+
 zero = ufloat(0, 0.1)
 zero2 = ufloat(0, 0.1)
 one = ufloat(1, 0.1)

uncertainties/ops.py

@@ -64,6 +64,28 @@ def wrapped_f(*args, **kwargs):
     return wrapped_f
 
 
+def pow_deriv_0(x, y):
+    """
+    The formula below works if x is positive or if y is an integer and x is negative
+    of y is an integer, x is zero and y is greater than or equal to 1.
+    """
+    if x > 0 or (y % 1 == 0 and (x < 0 or y >= 1)):
+        return y * x ** (y - 1)
+    elif x == 0 and y == 0:
+        return 0
+    else:
+        return float("nan")
+
+
+def pow_deriv_1(x, y):
+    if x > 0:
+        return log(x) * x**y
+    elif x == 0 and y > 0:
+        return 0
+    else:
+        return float("nan")
+
+
 def get_ops_with_reflection():
     """
     Return operators with a reflection, along with their partial derivatives.
@@ -119,29 +141,6 @@ def get_ops_with_reflection():
             eval("lambda y, x: %s" % expr) for expr in reversed(derivatives)
         ]
 
-    # The derivatives of pow() are more complicated:
-
-    # The case x**y is constant one the line x = 0 and in y = 0;
-    # the corresponding derivatives must be zero in these
-    # cases. If the function is actually not defined (e.g. 0**-3),
-    # then an exception will be raised when the nominal value is
-    # calculated.  These derivatives are transformed to NaN if an
-    # error happens during their calculation:
-
-    def pow_deriv_0(x, y):
-        if y == 0:
-            return 0.0
-        elif x != 0 or y % 1 == 0:
-            return y * x ** (y - 1)
-        else:
-            return float("nan")
-
-    def pow_deriv_1(x, y):
-        if x == 0 and y > 0:
-            return 0.0
-        else:
-            return log(x) * x**y
-
     ops_with_reflection["pow"] = [pow_deriv_0, pow_deriv_1]
     ops_with_reflection["rpow"] = [
         lambda y, x: pow_deriv_1(x, y),