Repository: lmfit/uncertainties

Issue:
Clean up power derivative calculations
`uncertainties` calculates the `x` and `y` derivatives of `pow(x, y) = x**y` using the `pow_deriv_0` and `pow_deriv_1` local functions defined in `get_ops_with_reflection` in `ops.py`.

https://github.com/lmfit/uncertainties/blob/19865cc0e87558d1d105c9d674ebace501e41395/uncertainties/ops.py#L131-L143

My issues with `pow_deriv_0`:

- It returns complex results for `x<0` with non-integer `y` even though `uncertainties` only supports real values. I would prefer it to return `nan` or raise a `ValueError` (weak preference for returning `nan`).
- For `x=0`, `y<0` and `y` an integer it raises a `ZeroDivisionError`, but for `x=0`, `y<0` and `y` a non-integer it returns `nan`. I think for `x=0` and `y<0` this function should always return `nan` or should always raise an exception (weak preference for returning `nan`).

My issues with `pow_deriv_1`:

- It raises a `ValueError` for all `x<0`. This is different than `pow_deriv_0` which returns complex results for this case. It should have the same behavior as `pow_deriv_0`, either return `nan` (my weak preference) or raise a `ValueError`.
- For `x=0` and `y<=0` it always raises `ValueError`. This is ok but again it's not consistent with the `pow_deriv_0` behavior. The two should be consistent, my preference would be to return `nan` or raise `ValueError`.

Other issues with both:
 
- But one big issue with both of these is that when the `ValueError` is raised, it is only raised as a consequence of the derivative calculation and is raised by built-in functions. It would be preferable to explicitly slice up the domain and handle all cases explicitly rather than relying on implicit behavior of the functions involved.
- They are both defined locally in the scope of another function `get_ops_with_reflection` in `ops.py`. They should be defined top-level in the module. This will allow direct testing to help clarify these edge cases.

Discussion:

Solution patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index cee272e..529ac3c 100644
--- a/CHANGES.rst
+++ b/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:
 
diff --git a/doc/user_guide.rst b/doc/user_guide.rst
index 9226bd9..8fee187 100644
--- a/doc/user_guide.rst
+++ b/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
 ======================
 
diff --git a/uncertainties/ops.py b/uncertainties/ops.py
index d390459..53e64c8 100644
--- a/uncertainties/ops.py
+++ b/uncertainties/ops.py
@@ -64,6 +64,28 @@ def nan_if_exception(f):
     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),


Test patch:
diff --git a/tests/test_power.py b/tests/test_power.py
index f89f38e..82d00e9 100644
--- a/tests/test_power.py
+++ b/tests/test_power.py
@@ -3,11 +3,40 @@ from math import pow as math_pow
 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)