diff --git a/fiddy/constants.py b/fiddy/constants.py index 36c3297..f017b9b 100644 --- a/fiddy/constants.py +++ b/fiddy/constants.py @@ -62,3 +62,6 @@ class MethodId(str, Enum): # FORWARD = MethodId.FORWARD # HYBRID = MethodId.HYBRID ##>>>>>>> origin/main + + +EPSILON = 1e-5 diff --git a/fiddy/derivative.py b/fiddy/derivative.py index 49d13a7..364b56f 100644 --- a/fiddy/derivative.py +++ b/fiddy/derivative.py @@ -1,5 +1,6 @@ import abc from typing import Any, Callable, Dict, List, Union +import warnings from dataclasses import dataclass @@ -9,6 +10,7 @@ from .constants import ( MethodId, Type, + EPSILON, ) from .analysis import Analysis @@ -133,6 +135,7 @@ def get_derivative( success_checker: Success, *args, analysis_classes: List[Analysis] = None, + relative_sizes: bool = False, directions: Union[List[Type.DIRECTION], Dict[str, Type.DIRECTION]] = None, direction_ids: List[str] = None, direction_indices: List[int] = None, @@ -149,14 +152,22 @@ def get_derivative( The IDs of the directions. directions: List: The directions to step along. Dictionary: keys are direction IDs, values are directions. + relative_sizes: + If `True`, sizes are scaled by the `point`, otherwise not. """ # TODO docs - direction_ids, directions = get_directions( - point=point, - directions=directions, - ids=direction_ids, - indices=direction_indices, - ) + if directions is not None: + direction_ids, directions = get_directions( + directions=directions, + ids=direction_ids, + indices=direction_indices, + ) + else: + direction_ids, directions = get_directions( + point=point, + ids=direction_ids, + indices=direction_indices, + ) if custom_methods is None: custom_methods = {} if analysis_classes is None: @@ -184,6 +195,7 @@ def get_derivative( size=size, method=method, autorun=False, + relative_size=relative_sizes, ) computers.append(computer) directional_derivative = DirectionalDerivative( diff --git a/fiddy/directional_derivative.py b/fiddy/directional_derivative.py index 89a5193..736b23f 100644 --- a/fiddy/directional_derivative.py +++ b/fiddy/directional_derivative.py @@ -1,5 +1,6 @@ import abc from typing import Any, Callable, Dict, List, Union, Tuple +import warnings import numpy as np import pandas as pd @@ -7,6 +8,7 @@ from .constants import ( MethodId, Type, + EPSILON, ) from .step import step @@ -36,6 +38,7 @@ class Computer: completed: bool = False results: List[ComputerResult] = field(default_factory=list) #value: Type.DIRECTIONAL_DERIVATIVE = None + relative_size: bool = False def __post_init__(self): if isinstance(self.method, MethodId): @@ -48,13 +51,30 @@ def __post_init__(self): if self.autorun: self() + def get_size(self): + if not self.relative_size: + return self.size + + # If relative, project point onto direction as scaling factor for size + unit_direction = self.direction / np.linalg.norm(self.direction) + # TODO add some epsilon to size? + size = np.dot(self.point, unit_direction) * self.size + if size == 0: + warnings.warn( + "Point has no component in this direction. " + "Set `Computer.relative_size=False` to avoid this. " + f"Using default small step size `fiddy.EPSILON`: {EPSILON}" + ) + size = EPSILON + return size + def __call__(self): value = self.method( point=self.point, direction=self.direction, - size=self.size, + size=self.get_size(), ) - result = ComputerResult(method_id=self.method.id, value=value, metadata={'size': self.size}) + result = ComputerResult(method_id=self.method.id, value=value, metadata={'size': self.get_size(), 'size_absolute': self.size}) self.results.append(result) self.completed = True diff --git a/fiddy/step.py b/fiddy/step.py index 362b344..f71c67a 100644 --- a/fiddy/step.py +++ b/fiddy/step.py @@ -22,4 +22,4 @@ def step( Returns: The step. """ - return direction * size / np.linalg.norm(direction) + return direction / np.linalg.norm(direction) * size diff --git a/tests/test_derivative.py b/tests/test_derivative.py index 62f2586..8585d0f 100644 --- a/tests/test_derivative.py +++ b/tests/test_derivative.py @@ -112,3 +112,64 @@ def test_get_derivative(point, sizes, output_shape): ) result = check(rtol=1e-2) assert result.success + + +def test_get_derivative_relative(): + point = np.array((3,4,0)) + size = 1e-1 + output_shape = (1,) + + direction = np.array([1,0,0]) + + directions = [direction] + success_checker = Consistency(atol=1e-2) + + function = partial(rosenbrock, output_shape=output_shape) + + # Expected finite difference derivatives + f_0 = function(point) + f_a = function(point + direction*size) + f_r = function(point + point*direction*size) # cardinal direction, simplifies to this, but usually need dot product + + g_a = (f_a-f_0)/size + g_r = (f_r-f_0)/(point*direction*size).sum() # cardinal direction, simplifies to this, but usually need dot product + + # Fiddy finite difference derivatives + kwargs = { + 'function': function, + 'point': point, + 'sizes': [size], + 'method_ids': [MethodId.FORWARD], + 'directions': [direction], + 'success_checker': success_checker, + } + fiddy_r = float(np.squeeze(get_derivative(**kwargs, relative_sizes=True).value)) + fiddy_a = float(np.squeeze(get_derivative(**kwargs, relative_sizes=False).value)) + + # Relative step sizes work + assert np.isclose(fiddy_r, g_r) + assert np.isclose(fiddy_a, g_a) + + # Same thing, now with non-cardinal direction + function = lambda x: (x[0]-2)**2 + (x[1]+3)**2 + point = np.array([3,4]) + direction = np.array([1,1]) + unit_direction = direction / np.linalg.norm(direction) + kwargs['function'] = function + kwargs['directions'] = [direction] + kwargs['point'] = point + + size_r = size * np.dot(point, unit_direction) + + f_0 = function(point) + f_a = function(point + unit_direction * size) + f_r = function(point + unit_direction * size_r) + + g_a = (f_a-f_0)/size + g_r = (f_r-f_0)/size_r + + fiddy_r = float(np.squeeze(get_derivative(**kwargs, relative_sizes=True).value)) + fiddy_a = float(np.squeeze(get_derivative(**kwargs, relative_sizes=False).value)) + + assert np.isclose(fiddy_r, g_r) + assert np.isclose(fiddy_a, g_a)