Source code for threeML.utils.differentiation

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 # 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
[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