Source code for openmdao.surrogate_models.nn_interpolators.linear_interpolator

"""Define the LinearInterpolator class."""

import numpy as np

from openmdao.surrogate_models.nn_interpolators.nn_base import NNBase
from six.moves import range


[docs]class LinearInterpolator(NNBase): """ Interpolate values by forming a hyperplane between the points closest to the prescribed inputs. """ def _find_hyperplane(self, neighbor_idx): """ Find hyperplane. Parameters ---------- neighbor_idx : ndarray Location indices for nearest neighbors. Returns ------- ndarray Normal vector for the hyperplane. ndarray Constant of the hyperplane. """ # Extra Inputs for Finding the normal are found below # Number of row vectors needed always dimensions - 1 indep_dims = self._indep_dims dep_dims = self._dep_dims # Number of Prediction Points nppts = neighbor_idx.shape[0] # Preallocate storage pc = np.empty((nppts, dep_dims), dtype='float') normal = np.empty((nppts, indep_dims + 1, dep_dims), dtype='float') nvect = np.empty((nppts, indep_dims, indep_dims + 1), dtype='float') trnd = np.concatenate((self._tp[neighbor_idx, :], self._tv[neighbor_idx, 0].reshape(nppts, indep_dims + 1, 1)), axis=2) nvect[:, :, :-1] = trnd[:, 1:, :-1] - trnd[:, :-1, :-1] for i in range(dep_dims): # Planar vectors need both dep and ind dimensions trnd[:, :, -1] = self._tv[neighbor_idx, i] # Go through each neighbor # Creates array[neighbor, dimension] from NN results nvect[:, :, -1] = trnd[:, 1:, -1] - trnd[:, :-1, -1] # Normal vector is in the null space of nvect. # Since nvect is of size indep x (indep + 1), # the normal vector will be the last entry in # V in the U, Sigma, V = svd(nvect). normal[:, :, i] = np.linalg.svd(nvect)[2][:, -1, :] # Use the point of the closest neighbor to # solve for pc - the constant of the n-dimensional plane. pc[:, i] = np.einsum('ij,ij->i', trnd[:, 0, :], normal[:, :, i]) return normal, pc def __call__(self, prediction_points): """ Do linear interpolation by defining plane with set number of nearest neighbors to predicted. Parameters ---------- prediction_points : ndarray Points at which predictions are computed. Returns ------- ndarray Predicted values at requested points. """ if len(prediction_points.shape) == 1: # Reshape vector to n x 1 array prediction_points.shape = (1, prediction_points.shape[0]) normalized_pts = (prediction_points - self._tpm) / self._tpr # Linear interp only uses as many neighbors as it has dimensions points_needed = self._indep_dims + 1 # KData query takes (data, #ofneighbors) to determine closest # training points to predicted data ndist, nloc = self._KData.query(normalized_pts.real, points_needed) normal, pc = self._find_hyperplane(nloc) # Set all predictions from values on plane predictions = np.einsum('ij,ijk->ik', normalized_pts, normal[:, :self._indep_dims, :]) - pc # Check to see if there are any collinear points and replace them n0 = np.where(normal[:, -1, :] == 0) predictions[n0, :] = self._tv[nloc[0, n0], :] # Finish computation for the good normals n = np.where(normal[:, -1, :] != 0) predictions[n] /= -normal[:, -1, :][n] # Rescale to original units predictions = (predictions * self._tvr) + self._tvm self._pt_cache = (normalized_pts, ndist, nloc) return predictions
[docs] def gradient(self, prediciton_points): """ Find the gradient at each location of a set of supplied predicted points. Parameters ---------- prediciton_points : ndarray Prediction points at which the gradient is calculated. Returns ------- ndarray Gradient at the prediction points. """ if len(prediciton_points.shape) == 1: # Reshape vector to n x 1 array prediciton_points.shape = (1, prediciton_points.shape[0]) normPredPts = (prediciton_points - self._tpm) / self._tpr nppts = normPredPts.shape[0] gradient = np.zeros((nppts, self._dep_dims, self._indep_dims), dtype="float") # Linear interp only uses as many neighbors as it has dimensions dims = self._indep_dims + 1 # Find the neighbors if self._pt_cache is not None and \ np.allclose(self._pt_cache[0], normPredPts): ndist, nloc = self._pt_cache[1:] else: ndist, nloc = self._KData.query(normPredPts.real, dims) normal, pc = self._find_hyperplane(nloc) if np.any(normal[:, -1, :]) == 0: return gradient gradient[:] = (-normal[:, :-1, :] / normal[:, -1, :]).squeeze().T grad = gradient * (self._tvr[:, np.newaxis] / self._tpr) return grad