import numdifftools as nd
import numpy as np
from astromodels import SettingOutOfBounds
[docs]
class ParameterOnBoundary(RuntimeError):
pass
[docs]
class CannotComputeHessian(RuntimeError):
pass
def _get_wrapper(function, point, minima, maxima):
point = np.array(point, ndmin=1, dtype=float)
minima = np.array(minima, ndmin=1, dtype=float)
maxima = np.array(maxima, ndmin=1, dtype=float)
n_dim = point.shape[0]
# Find order of magnitude of each coordinate. If one of the coordinates is exactly zero we need
# to treat it differently
idx = point == 0.0
orders_of_magnitude = np.zeros_like(point)
orders_of_magnitude[idx] = 1.0
orders_of_magnitude[~idx] = 10 ** np.ceil(
np.log10(np.abs(point[~idx]))
) # type: np.ndarray
scaled_point = point / orders_of_magnitude
scaled_minima = minima / orders_of_magnitude
scaled_maxima = maxima / orders_of_magnitude
# Decide a delta for the finite differentiation
# The algorithm implemented in numdifftools is robust with respect to the choice
# of delta, as long as we are not going beyond the boundaries (which would cause
# the procedure to fail)
scaled_deltas = np.zeros_like(scaled_point)
for i in range(n_dim):
scaled_value = scaled_point[i]
scaled_min_value, scaled_max_value = (scaled_minima[i], scaled_maxima[i])
if scaled_value == scaled_min_value or scaled_value == scaled_max_value:
raise ParameterOnBoundary(
"Value for parameter number %s is on the boundary" % i
)
if not np.isnan(scaled_min_value):
# Parameter with low bound
distance_to_min = scaled_value - scaled_min_value
else:
# No defined minimum
distance_to_min = np.inf
if not np.isnan(scaled_max_value):
# Parameter with hi bound
distance_to_max = scaled_max_value - scaled_value
else:
# No defined maximum
distance_to_max = np.inf
# Delta is the minimum between 0.03% of the value, and 1/2.5 times the minimum
# distance to either boundary. 1/2 of that factor is due to the fact that numdifftools uses
# twice the delta to compute the differential, and the 0.5 is due to the fact that we don't want
# to go exactly equal to the boundary
if scaled_point[i] == 0.0:
scaled_deltas[i] = min([1e-5, distance_to_max / 2.5, distance_to_min / 2.5])
else:
scaled_deltas[i] = min(
[
0.003 * abs(scaled_point[i]),
distance_to_max / 2.5,
distance_to_min / 2.5,
]
)
def wrapper(x):
scaled_back_x = x * orders_of_magnitude # type: np.ndarray
try:
result = function(*scaled_back_x)
except SettingOutOfBounds:
raise CannotComputeHessian(
"Cannot compute Hessian, parameters out of bounds at %s" % scaled_back_x
)
else:
return result
return wrapper, scaled_deltas, scaled_point, orders_of_magnitude, n_dim
[docs]
def get_jacobian(function, point, minima, maxima):
wrapper, scaled_deltas, scaled_point, orders_of_magnitude, n_dim = _get_wrapper(
function, point, minima, maxima
)
# Compute the Jacobian matrix at best_fit_values
jacobian_vector = nd.Jacobian(wrapper, scaled_deltas, method="central")(
scaled_point
)
# Transform it to numpy matrix
jacobian_vector = np.array(jacobian_vector)
# Now correct back the Jacobian for the scales
jacobian_vector /= orders_of_magnitude
return jacobian_vector[0]
[docs]
def get_hessian(function, point, minima, maxima):
wrapper, scaled_deltas, scaled_point, orders_of_magnitude, n_dim = _get_wrapper(
function, point, minima, maxima
)
# Compute the Hessian matrix at best_fit_values
hessian_matrix_ = nd.Hessian(wrapper, scaled_deltas)(scaled_point)
# Transform it to numpy matrix
hessian_matrix = np.array(hessian_matrix_)
# Now correct back the Hessian for the scales
for i in range(n_dim):
for j in range(n_dim):
hessian_matrix[i, j] /= orders_of_magnitude[i] * orders_of_magnitude[j]
return hessian_matrix