"""Newton-Raphson calculator"""
import os
from typing import Any
from collections.abc import Callable
if os.environ.get('__USE_CUDA', '0') == '1':
import cupy as np
else:
import numpy as np
from scipy import optimize
from g3d.defs import G3DCode
from g3d.defs import G3DDisplay
from g3d.defs import CFlowSolverMode
[docs]
class G3DNewtonRaphsonSolver(G3DCode):
"""Newton-Raphson solver"""
def __init__(self, *args: Any, **kwargs: Any) -> None:
"""Initialize Newton-Raphson solver parameters"""
# default values
self.__default_solver_mode = CFlowSolverMode.NATIVE_NEWTON
self.__default_end_value = 0.1
self.__default_target_value = None
self.__default_start_value = None
self.__default_delta = 0.0001
self.__default_tolerance = 1.e-6
self.__default_max_iter = 100
self.__default_coarse_delta = 0.01
self.__default_verbose = False
self.__default_lower_bound = 1.
self.__default_upper_bound = 1.
self.__default_check_upper_bound = False
self.__default_check_lower_bound = False
# initialize parameters
self._solver_mode = self.__default_solver_mode
self._end_value = self.__default_end_value
self._target_value = self.__default_target_value
self._start_value = self.__default_start_value
self._delta = self.__default_delta
self._tolerance = self.__default_tolerance
self._max_iter = self.__default_max_iter
self._coarse_delta = self.__default_coarse_delta
self._verbose = self.__default_verbose
self._lower_bound = self.__default_lower_bound
self._upper_bound = self.__default_upper_bound
self._check_upper_bound = self.__default_check_upper_bound
self._check_lower_bound = self.__default_check_lower_bound
self._scipy_kwargs = {}
super().__init__(*args, **kwargs)
def _set_default(self) -> None:
super()._set_default()
self._tag = 'g3d_newton_raphson_solver'
self.reset()
[docs]
def reset(self) -> None:
"""Resets Newton-Raphson solver parameters to their default values"""
self._solver_mode = self.__default_solver_mode
self._end_value = self.__default_end_value
self._target_value = self.__default_target_value
self._start_value = self.__default_start_value
self._delta = self.__default_delta
self._tolerance = self.__default_tolerance
self._max_iter = self.__default_max_iter
self._coarse_delta = self.__default_coarse_delta
self._verbose = self.__default_verbose
self._lower_bound = self.__default_lower_bound
self._upper_bound = self.__default_upper_bound
self._check_upper_bound = self.__default_check_upper_bound
self._check_lower_bound = self.__default_check_lower_bound
self.__reset_scipy_kwargs()
def __reset_scipy_kwargs(self, kwargs_dict: dict[Any, Any] = None) -> None:
"""Resets scipy parameter dict"""
self._scipy_kwargs = {}
if kwargs_dict is not None:
self.__add_scipy_kwargs(kwargs_dict)
def __set_scipy_kwargs(self, kwargs_dict: dict[Any, Any]) -> None:
"""Updates scipy kwarg dict"""
self.__reset_scipy_kwargs(kwargs_dict)
def __add_scipy_kwargs(self, kwargs_dict: dict[Any, Any]) -> None:
for key, value in kwargs_dict.items():
self._scipy_kwargs[key] = value
[docs]
def update_settings(self, **kwargs: Any) -> None:
"""Updates Newton-Raphson solver parameters"""
for key, value in kwargs.items():
self.set_from_key_value(key, value)
[docs]
def get_settings_from_solver(self, newton: 'G3DNewtonRaphsonSolver') -> None:
"""
copy settings from provided Newton-Raphson solver object
input: Newton-Raphson solver object
"""
self._solver_mode = newton.solver_mode
self._end_value = newton.end_value
self._target_value = newton.target_value
self._start_value = newton.start_value
self._delta = newton.delta
self._tolerance = newton.tolerance
self._max_iter = newton.max_iter
self._coarse_delta = newton.coarse_delta
self._verbose = newton.verbose
self._lower_bound = newton.lower_bound
self._upper_bound = newton.upper_bound
self._check_lower_bound = newton.check_lower_bound
self._check_upper_bound = newton.check_upper_bound
self.__reset_scipy_kwargs(newton.scipy_kwargs)
@property
def solver_mode(self) -> Any:
"""Newton-Raphson solver mode: native (default) or scipy"""
return self._solver_mode
@solver_mode.setter
def solver_mode(self, val: Any) -> None:
self._solver_mode = val
@property
def end_value(self) -> float:
"""Search end value"""
return self._end_value
@end_value.setter
def end_value(self, val: float) -> None:
self._end_value = val
@property
def target_value(self) -> float:
"""Function target value"""
return self._target_value
@target_value.setter
def target_value(self, val: float) -> None:
self._target_value = val
@property
def start_value(self) -> float:
"""Newton search start value"""
return self._start_value
@start_value.setter
def start_value(self, val: float) -> None:
self._start_value = val
@property
def delta(self) -> float:
"""Delta used for numerical differentiation (default 1.e-4)"""
return self._delta
@delta.setter
def delta(self, val: float) -> None:
self._delta = val
@property
def tolerance(self) -> float:
"""Newton-Raphson tolerance = abs(value - target_value), (default 1.e-6)"""
return self._tolerance
@tolerance.setter
def tolerance(self, val: float) -> None:
self._tolerance = val
@property
def max_iter(self) -> int:
"""Maximum number of Newton-Raphson iterations (default 100)"""
return self._max_iter
@max_iter.setter
def max_iter(self, val: int) -> None:
self._max_iter = val
@property
def coarse_delta(self) -> float:
"""Delta used in coarse search (default 1.e-2)"""
return self._coarse_delta
@coarse_delta.setter
def coarse_delta(self, val: float) -> None:
self._coarse_delta = val
@property
def verbose(self) -> bool:
"""Verbose mode => prints iteration information (default False)"""
return self._verbose
@verbose.setter
def verbose(self, val: bool) -> None:
self._verbose = val
@property
def lower_bound(self) -> float:
"""Lower function bound, one-sided numerical differentiation used if needed"""
return self._lower_bound
@lower_bound.setter
def lower_bound(self, val: float) -> None:
self._lower_bound = val
self._check_lower_bound = True
@property
def upper_bound(self) -> float:
"""Upper function bound, one-sided numerical differentiation used if needed"""
return self._upper_bound
@upper_bound.setter
def upper_bound(self, val: float) -> None:
self._upper_bound = val
self._check_upper_bound = True
@property
def check_lower_bound(self) -> bool:
"""Check lower bound (bool, default value: False)"""
return self._check_lower_bound
@check_lower_bound.setter
def check_lower_bound(self, val: bool) -> None:
self._check_lower_bound = val
@property
def check_upper_bound(self) -> bool:
"""Check upper bound (bool, default value: False)"""
return self._check_upper_bound
@check_upper_bound.setter
def check_upper_bound(self, val: bool) -> None:
self._check_upper_bound = val
@property
def scipy_kwargs(self) -> dict[Any, Any]:
"""Scipy Newton-Raphson parameters"""
return self._scipy_kwargs
@scipy_kwargs.setter
def scipy_kwargs(self, kwargs_dict: dict[Any, Any]) -> None:
self.__set_scipy_kwargs(kwargs_dict)
def __set_from_kwargs(self, **kwargs: Any) -> None:
for key, value in kwargs.items():
self.__set_from_key_value(key, value)
[docs]
def set_from_key_value(self, key: str, value: Any) -> bool:
"""Set parameters from provided key => value input"""
if super().set_from_key_value(key, value):
return True
return self.__set_from_key_value(key, value)
def __set_from_key_value(self, key: str, value: Any) -> bool:
if key == 'nrs':
if type(value) is G3DNewtonRaphsonSolver:
self.get_settings_from_solver(value)
return True
elif key == 'solver_mode':
if type(value) is str:
if value == 'native_newton' or value == 'native':
self._solver_mode = CFlowSolverMode.NATIVE_NEWTON
if value == 'scipy_newton' or value == 'scipy':
self._solver_mode = CFlowSolverMode.SCIPY_NEWTON
else:
self._solver_mode = value
return True
elif key == 'end_value' or key == 'end':
self._end_value = value
return True
elif key == 'target_value' or key == 'target':
self._target_value = value
return True
elif key == 'start_value' or key == 'start':
self._start_value = value
self._scipy_kwargs['x0'] = value
return True
elif key == 'delta' or key == 'h':
self._delta = value
return True
elif key == 'tolerance' or key == 'tol':
self._tolerance = value
self._scipy_kwargs['tol'] = value
return True
elif key == 'max_iter' or key == 'maxiter':
self._max_iter = value
self._scipy_kwargs['maxiter'] = value
return True
elif key == 'coarse_delta' or key == 'coarse':
self._coarse_delta = value
return True
elif key == 'verbose_newton':
self._verbose = value
return True
elif key == 'lower_bound':
self._lower_bound = value
self._check_lower_bound = True
return True
elif key == 'upper_bound':
self._upper_bound = value
self._check_upper_bound = True
return True
elif key in ['fprime', 'fprime2', 'rtol', 'x1', 'full_output', 'disp']:
self._scipy_kwargs[key] = value
return True
return False
[docs]
def print(self) -> None:
"""Prints solver settings"""
if self._start_value is not None:
print("start value: %f" % float(self._start_value))
if self._target_value is not None:
print("target value: %f" % float(self._target_value))
print("delta: %e" % self._delta)
print("tolerance: %e" % self._tolerance)
if self._check_lower_bound:
print("lower_bound: %e" % self._lower_bound)
if self._check_upper_bound:
print("upper_bound: %e" % self._upper_bound)
print("max iterations: %d" % self._max_iter)
[docs]
def find_start_value(self, func: Callable[[float], float], *args: Any, **kwargs: Any) -> float:
"""Finds start value using coarse search for provided function"""
self.set(*args, **kwargs)
# find a start value by doing a coarse search
if self._end_value > self._start_value:
while self._start_value < self._end_value:
if func(self._start_value) > self._target_value:
break
self._start_value += self._coarse_delta
self._start_value -= self._coarse_delta
else:
while self._start_value > self._end_value:
if func(self._start_value) > self._target_value:
break
self._start_value -= self._coarse_delta
self._start_value += self._coarse_delta
return self._start_value
[docs]
def solve(self, func: Callable[[float], float], *args: Any, **kwargs: Any) -> float:
"""Executes scipy Newton-Raphson for provided function and parameters"""
self.__set_from_kwargs(**kwargs)
if self._solver_mode == CFlowSolverMode.NATIVE_NEWTON:
return self.native_solve(func, *args, **kwargs)
if self._solver_mode == CFlowSolverMode.SCIPY_NEWTON:
self._scipy_kwargs['x0'] = self._start_value
def n_func(x):
"""
:param x:
:return:
"""
return func(x) - self._target_value
return optimize.newton(n_func, **self._scipy_kwargs)
return -1
[docs]
def native_solve(self, func: Callable[[float], float], *args: Any, **kwargs: Any) -> float:
"""Executes native Newton-Raphson for provided function and parameters"""
self.set(*args, **kwargs)
g3d_display = G3DDisplay()
x = self._start_value
delta = self._delta
i = 1
while True:
F = func(x) - self._target_value
dm = delta
dp = delta
if self._check_lower_bound:
dm = min([dm, abs(x - self._lower_bound)])
if self._check_upper_bound:
dp = min([dp, abs(self._upper_bound - x)])
dF = (func(x + dp) - func(x - dm)) / (dm + dp)
if dF == 0.:
# dF = 1.e-20
delta *= 2
else:
x = x - (F / dF)
if g3d_display.verbose or self._verbose:
print(' iteration: %7d, value: %12.7E, (F= %12.7E, dF= %12.7E, delta= %10.5E)' %
(i, x, F, dF, dm + dp))
if np.abs(F) < self._tolerance or i >= self._max_iter:
break
i += 1
return x