Hybrid grad check

fiddy/derivative_check.py

@@ -1,11 +1,12 @@
 import abc
 from typing import Any, Callable, Dict, List, Union
-from itertools import product
+from itertools import chain
 
 from dataclasses import dataclass
 
 import numpy as np
 import pandas as pd
+import math
 
 from .constants import (
     # TYPE_DIMENSION,
@@ -161,3 +162,109 @@ def method(self, *args, **kwargs):
             success=success,
         )
         return derivative_check_result
+
+
+class HybridDerivativeCheck(DerivativeCheck):
+    method_id = "hybrid"
+
+    def method(self, *args, **kwargs):
+        expected_values = []
+        test_values = []
+        success = True
+        for direction_index, directional_derivative in enumerate(
+            self.derivative.directional_derivatives
+        ):
+            test_value = directional_derivative.value
+            test_values.append(test_value)
+
+            expected_value = []
+            for output_index in np.ndindex(self.output_indices):
+                element = self.expectation[output_index][direction_index]
+                expected_value.append(element)
+            expected_value = np.array(expected_value).reshape(test_value.shape)
+            expected_values.append(expected_value)
+
+        # debug
+        assert len(expected_values) == len(
+            test_values
+        ), "Mismatch of step sizes"
+
+        results_all = []
+        directional_derivative_check_results = []
+        for step_size in range(0, len(expected_values)):
+            approxs_for_param = []
+            grads_for_param = []
+            results = []
+            for diff_index, directional_derivative in enumerate(
+                self.derivative.directional_derivatives
+            ):
+                try:
+                    for grad, approx in zip(
+                        expected_values[diff_index - 1][step_size - 1],
+                        test_values[diff_index - 1][step_size - 1],
+                    ):
+                        approxs_for_param.append(approx)
+                        grads_for_param.append(grad)
+                    fd_range = np.percentile(approxs_for_param, [0, 100])
+                    fd_mean = np.mean(approxs_for_param)
+                    grad_mean = np.mean(grads_for_param)
+                    if not (fd_range[0] <= grad_mean <= fd_range[1]):
+                        if np.any(
+                            [
+                                abs(x - y) > kwargs["atol"]
+                                for i, x in enumerate(approxs_for_param)
+                                for j, y in enumerate(approxs_for_param)
+                                if i != j
+                            ]
+                        ):
+                            fd_range = abs(fd_range[1] - fd_range[0])
+                            if (
+                                abs(grad_mean - fd_mean)
+                                / abs(fd_range + np.finfo(float).eps)
+                            ) > kwargs["rtol"]:
+                                results.append(False)
+                            else:
+                                results.append(False)
+                        else:
+                            results.append(
+                                None
+                            )  # can't judge consistency / questionable grad approxs
+                    else:
+                        fd_range = abs(fd_range[1] - fd_range[0])
+                        if math.isinf(
+                            (fd_range)
+                            or math.isnan(fd_range)
+                            or math.isinf(fd_mean)
+                            or math.isnan(fd_mean)
+                        ):
+                            results.append(None)
+                        else:
+                            results.append(True)
+                except (IndexError, TypeError):
+                    # TODO: Fix this, why does this occur?
+                    pass
+
+                directional_derivative_check_result = (
+                    DirectionalDerivativeCheckResult(
+                        direction_id=directional_derivative.id,
+                        method_id=self.method_id,
+                        test=test_value,
+                        expectation=expected_value,
+                        output={"return": results},
+                        success=all(results),
+                    )
+                )
+                directional_derivative_check_results.append(
+                    directional_derivative_check_result
+                )
+                results_all.append(results)
+
+        success = all(chain(*results_all))
+        derivative_check_result = DerivativeCheckResult(
+            method_id=self.method_id,
+            directional_derivative_check_results=directional_derivative_check_results,
+            test=self.derivative.value,
+            expectation=self.expectation,
+            success=success,
+        )
+        return derivative_check_result

fiddy/gradient_check.py

@@ -140,7 +140,7 @@ def gradient_check(
 
     if sort:
         results_df.sort_values(
-            by=['dimension', 'size'],
+            by=["dimension", "size"],
             inplace=True,
         )
 
@@ -266,7 +266,7 @@ def keep_lowest_error(
     if not inplace:
         sort_df = minimal_results_df
     sort_df.sort_values(
-        ['success', 'dimension', 'size'],
+        ["success", "dimension", "size"],
         ascending=[False, True, True],
         inplace=True,
     )

fiddy/success.py

@@ -94,8 +94,8 @@ def method(
             success_by_size[size] = np.isclose(
                 values,
                 np.nanmean(values, axis=0),
-                rtol=self.rtol/2,
-                atol=self.atol/2,
+                rtol=self.rtol / 2,
+                atol=self.atol / 2,
                 equal_nan=self.equal_nan,
             ).all()
 
@@ -108,13 +108,16 @@ def method(
         success = False
         value = np.nanmean(np.array(consistent_results), axis=0)
         if consistent_results:
-            success = np.isclose(
-                consistent_results,
-                value,
-                rtol=self.rtol,
-                atol=self.atol,
-                equal_nan=self.equal_nan
-            ).all() and not np.isnan(consistent_results).all()
+            success = (
+                np.isclose(
+                    consistent_results,
+                    value,
+                    rtol=self.rtol,
+                    atol=self.atol,
+                    equal_nan=self.equal_nan,
+                ).all()
+                and not np.isnan(consistent_results).all()
+            )
         value = np.average(np.array(consistent_results), axis=0)
 
         return success, value

tests/test_derivative.py

@@ -10,7 +10,10 @@
 from fiddy.derivative import Computer
 from fiddy.analysis import ApproximateCentral
 from fiddy.success import Consistency
-from fiddy.derivative_check import NumpyIsCloseDerivativeCheck
+from fiddy.derivative_check import (
+    NumpyIsCloseDerivativeCheck,
+    HybridDerivativeCheck,
+)
 
 
 RTOL = 1e-2
@@ -108,7 +111,7 @@ def test_get_derivative(point, sizes, output_shape):
         # FIXME default?
         sizes=[1e-10, 1e-5],
         # FIXME default?
-        method_ids=[MethodId.FORWARD, MethodId.BACKWARD],
+        method_ids=[MethodId.FORWARD, MethodId.BACKWARD, MethodId.CENTRAL],
         # FIXME default?
         analysis_classes=[ApproximateCentral],
         # FIXME default? not just "True" ...
@@ -122,7 +125,55 @@ def test_get_derivative(point, sizes, output_shape):
         expectation=expected_value,
         point=point,
     )
-    result = check(rtol=1e-2)
+    result = check(rtol=1e-2, atol=1e-3)
+    assert result.success
+
+
+@pytest.mark.parametrize(
+    "point, sizes, output_shape",
+    [
+        (np.array(point), sizes, output_shape)
+        for point in [
+            (1, 0, 0),
+            (0.9, 0.1, 0.2, 0.4),
+        ]
+        for sizes in [
+            [1e-10, 1e-5],
+        ]
+        for output_shape in [
+            (1,),
+            (1, 2),
+            (5, 3, 6, 2, 4),
+        ]
+    ],
+)
+def test_get_derivative_hybrid(point, sizes, output_shape):
+    function = partial(rosenbrock, output_shape=output_shape)
+    expected_derivative_function = partial(
+        rosenbrock_der, output_shape=output_shape
+    )
+    derivative = get_derivative(
+        function=function,
+        point=point,
+        # FIXME default?
+        sizes=[1e-10, 1e-5],
+        # FIXME default?
+        method_ids=[MethodId.FORWARD, MethodId.BACKWARD, MethodId.CENTRAL],
+        # FIXME default?
+        analysis_classes=[ApproximateCentral],
+        # FIXME default? not just "True" ...
+        success_checker=Consistency(),
+    )
+    test_value = derivative.value
+    expected_value = expected_derivative_function(point)
+
+    check = HybridDerivativeCheck(
+        derivative=derivative,
+        expectation=expected_value,
+        point=point,
+    )
+    # based on given tolerances, hybrid gradient check should not fail
+    result = check(rtol=1e-2, atol=1e-3)
     assert result.success