Tests

Author: jrmaddison

pyadjoint/adjfloat.py

@@ -20,6 +20,20 @@ def register(cls):
     return register
 
 
+@register_function(np.power)
+class _pyadjoint_power(sp.Function):
+    def fdiff(self, argindex=1):
+        if argindex == 1:
+            return sp.Piecewise(
+                # Let SymPy decide how to handle indeterminate form
+                ((self.args[0] ** self.args[1]).diff(self.args[0]), sp.And(self.args[0] == 0, self.args[1] == 0)),
+                # Otherwise simplify
+                (sp.S.Zero, self.args[1] == 0),
+                (_pyadjoint_power(self.args[0], self.args[1] - 1) * self.args[1], True))
+        elif argindex == 2:
+            return (self.args[0] ** self.args[1]).diff(self.args[1])
+
+
 @register_function(np.hypot)
 class _pyadjoint_hypot(sp.Function):
     def fdiff(self, argindex=1):
@@ -228,11 +242,11 @@ def __sub__(self, other):
     def __rsub__(self, other):
         return super().__rsub__(other)
 
-    @annotate_operator(Operator(operator.pow, 2), operator.pow)
+    @annotate_operator(Operator(_pyadjoint_power, 2), operator.pow)
     def __pow__(self, other):
         return super().__pow__(other)
 
-    @annotate_operator(Operator(roperator(operator.pow), 2), roperator(operator.pow))
+    @annotate_operator(Operator(roperator(_pyadjoint_power), 2), roperator(operator.pow))
     def __rpow__(self, other):
         return super().__rpow__(other)
 
@@ -243,7 +257,7 @@ def __rpow__(self, other):
     divide = register_operator(np.divide, operator.truediv, 2)
     add = register_operator(np.add, operator.add, 2)
     subtract = register_operator(np.subtract, operator.sub, 2)
-    power = register_operator(np.power, operator.pow, 2)
+    power = register_operator(np.power, _pyadjoint_power, 2)
     minimum = register_operator(
         np.minimum,
         lambda self, other: sp.Piecewise((self, self <= other),

tests/pyadjoint/test_floats.py

@@ -2,6 +2,7 @@
 import math
 import numpy as np
 import operator
+import sympy as sp
 from numpy.testing import assert_approx_equal
 from numpy.random import rand
 from pyadjoint import *
@@ -165,6 +166,8 @@ def test_float_abs(v, abs):
     assert_approx_equal(b, abs(v))
 
     rf = ReducedFunctional(b, Control(a))
+    assert_approx_equal(rf(v), abs(v))
+    assert_approx_equal(rf.tlm(1.0), 1.0 if v >= 0 else -1.0)
     assert_approx_equal(rf.derivative(), 1.0 if v >= 0 else -1.0)
     assert_approx_equal(rf.hessian(1.0), 0.0)
 
@@ -177,6 +180,8 @@ def test_float_exp(v, exp):
     assert_approx_equal(b, math.exp(v))
 
     rf = ReducedFunctional(b, Control(a))
+    assert_approx_equal(rf(v), math.exp(v))
+    assert_approx_equal(rf.tlm(1.0), math.exp(v))
     assert_approx_equal(rf.derivative(), math.exp(v))
     assert_approx_equal(rf.hessian(1.0), math.exp(v))
 
@@ -189,10 +194,109 @@ def test_float_loglog(log):
     assert_approx_equal(b, math.log(math.log(v)))
 
     rf = ReducedFunctional(b, Control(a))
+    assert_approx_equal(rf(v), math.log(math.log(v)))
+    assert_approx_equal(rf.tlm(1.0), 1.0 / (v * math.log(v)))
     assert_approx_equal(rf.derivative(), 1.0 / (v * math.log(v)))
     assert_approx_equal(rf.hessian(1.0), -(1.0 + math.log(v)) / ((v * math.log(v)) ** 2))
 
 
+def compose(*args):
+    def fn(x):
+        for arg in reversed(args):
+            x = arg(x)
+        return x
+    return fn
+
+
+def sq(x):
+    return x ** 2
+
+
+def hypotsq(x):
+    return 1 + x ** 2
+
+
+@pytest.mark.parametrize("np_operator, sp_operator",
+                         (
+                             (operator.abs, lambda x: sp.Piecewise((x, x >= 0), (-x, True))),
+                             (operator.pos, operator.pos),
+                             (operator.neg, operator.neg),
+                             (lambda x: 1 + x, lambda x: 1 + x),
+                             (lambda x: x + 1, lambda x: 1 + x),
+                             (lambda x: 1 - x, lambda x: 1 - x),
+                             (lambda x: x - 1, lambda x: x - 1),
+                             (lambda x: 2 * x, lambda x: 2 * x),
+                             (lambda x: x * 2, lambda x: 2 * x),
+                             (compose(lambda x: 2 / x, hypotsq), compose(lambda x: 2 / x, hypotsq)),
+                             (lambda x: x / 2, lambda x: x / 2),
+                             (lambda x: 2 ** x, lambda x: 2 ** x),
+                             (lambda x: x ** 2, sq),
+                             (np.absolute, lambda x: sp.Piecewise((x, x >= 0), (-x, True))),
+                             (np.positive, operator.pos),
+                             (np.negative, operator.neg),
+                             (lambda x: np.add(x, 1), lambda x: 1 + x),
+                             (lambda x: np.add(1, x), lambda x: 1 + x),
+                             (lambda x: np.subtract(x, 1), lambda x: x - 1),
+                             (lambda x: np.subtract(1, x), lambda x: 1 - x),
+                             (lambda x: np.multiply(x, 2), lambda x: 2 * x),
+                             (lambda x: np.multiply(2, x), lambda x: 2 * x),
+                             (lambda x: np.divide(x, 2), lambda x: x / 2),
+                             (compose(lambda x: np.divide(2, x), hypotsq), compose(lambda x: 2 / x, hypotsq)),
+                             (lambda x: np.power(x, 2), sq),
+                             (lambda x: np.power(2, x), lambda x: 2 ** x),
+                             (np.sin, sp.sin),
+                             (np.cos, sp.cos),
+                             (np.tan, sp.tan),
+                             (compose(np.arcsin, np.tanh), compose(sp.asin, sp.tanh)),
+                             (compose(np.arccos, np.tanh), compose(sp.acos, sp.tanh)),
+                             (np.arctan, sp.atan),
+                             (lambda x: np.arctan2(x, 1), sp.atan),
+                             (lambda x: np.arctan2(1, x), lambda x: sp.atan2(1, x)),
+                             (compose(sq, lambda x: np.hypot(1, x)), hypotsq),
+                             (compose(sq, lambda x: np.hypot(x, 1)), hypotsq),
+                             (np.sinh, sp.sinh),
+                             (np.cosh, sp.cosh),
+                             (np.tanh, sp.tanh),
+                             (np.arcsinh, sp.asinh),
+                             (compose(np.arccosh, hypotsq, hypotsq), (compose(sp.acosh, hypotsq, hypotsq))),
+                             (compose(np.arctanh, np.sin), compose(sp.atanh, sp.sin)),
+                             (np.exp, sp.exp),
+                             (np.exp2, lambda x: 2 ** x),
+                             (np.expm1, lambda x: sp.exp(x) - 1),
+                             (compose(np.log, hypotsq), compose(sp.log, hypotsq)),
+                             (compose(np.log2, hypotsq), compose(lambda x: sp.log(x, 2), hypotsq)),
+                             (compose(np.log10, hypotsq), compose(lambda x: sp.log(x, 10), hypotsq)),
+                             (compose(np.log1p, sq), compose(sp.log, hypotsq)),
+                             (compose(np.sqrt, hypotsq), compose(sp.sqrt, hypotsq)),
+                             (np.square, sq),
+                             (compose(np.cbrt, hypotsq), compose(lambda x: x ** sp.Rational(1, 3), hypotsq)),
+                             (compose(np.reciprocal, hypotsq), compose(lambda x: 1 / x, hypotsq)),
+                             (lambda x: np.minimum(x, 1), lambda x: sp.Piecewise((x, x <= 1), (sp.S.One, True))),
+                             (lambda x: np.minimum(1, x), lambda x: sp.Piecewise((x, x <= 1), (sp.S.One, True))),
+                             (lambda x: np.maximum(x, 1), lambda x: sp.Piecewise((x, x >= 1), (sp.S.One, True))),
+                             (lambda x: np.maximum(1, x), lambda x: sp.Piecewise((x, x >= 1), (sp.S.One, True))),
+                         ))
+@pytest.mark.parametrize("v", (-np.sqrt(np.pi), 0, np.exp(0.5)))
+def test_float_operators(np_operator, sp_operator, v):
+    np_operator = compose(np.exp, np_operator)
+    sp_operator = compose(sp.exp, sp_operator)
+
+    a = AdjFloat(v)
+    b = np_operator(a)
+
+    x = sp.Symbol("x", real=True)
+    op_ref = sp.lambdify((x,), sp_operator(x), modules=["numpy"])(v)
+    dop_ref = sp.lambdify((x,), sp_operator(x).diff(x), modules=["numpy"])(v)
+    ddop_ref = sp.lambdify((x,), sp_operator(x).diff(x).diff(x), modules=["numpy"])(v)
+
+    assert_approx_equal(b, op_ref)
+    rf = ReducedFunctional(b, Control(a))
+    assert_approx_equal(rf(v), op_ref)
+    assert_approx_equal(rf.tlm(1.0), dop_ref)
+    assert_approx_equal(rf.derivative(), dop_ref)
+    assert_approx_equal(rf.hessian(1.0), ddop_ref)
+
+
 def test_float_logexp():
     a = AdjFloat(3.0)
     b = exp(a)
@@ -240,8 +344,6 @@ def test_float_exponentiation():
     # d(a**a)/da = dexp(a log(a))/da = a**a * (log(a) + 1)
     assert_approx_equal(rf.derivative(), 4.0 * (math.log(2.0)+1.0))
 
-    # TODO: __rpow__ is not yet implemented
-
 
 @pytest.mark.parametrize("B", [3,4])
 @pytest.mark.parametrize("E", [6,5])