Source code for bladex.ndinterpolator

import numpy as np
from scipy.spatial.distance import cdist
from scipy.interpolate import splev


[docs]class RBF(object): """ .. _ndinterpolator-label: RBF ndinterpolator ====================== Module focused on the implementation of the Radial Basis Functions interpolation technique. This technique is still based on the use of a set of parameters, the so-called control points, as for FFD, but RBF is interpolatory. Another important key point of RBF strategy relies in the way we can locate the control points: in fact, instead of FFD where control points need to be placed inside a regular lattice, with RBF we havo no more limitations. So we have the possibility to perform localized control points refiniments. :Theoretical Insight: As reference please consult M.D. Buhmann, Radial Basis Functions, volume 12 of Cambridge monographs on applied and computational mathematics. Cambridge University Press, UK, 2003. This implementation follows D. Forti and G. Rozza, Efficient geometrical parametrization techniques of interfaces for reduced order modelling: application to fluid-structure interaction coupling problems, International Journal of Computational Fluid Dynamics. RBF shape parametrization technique is based on the definition of a map :math:`\\mathcal{M}(\\boldsymbol{x}) : \\mathbb{R}^n \\rightarrow \\mathbb{R}^n`, that allows the possibility of transferring data across non-matching grids and facing the dynamic mesh handling. The map introduced is defines as follows .. math:: \\mathcal{M}(\\boldsymbol{x}) = p(\\boldsymbol{x}) + \\sum_{i=1}^{\\mathcal{N}_C} \\gamma_i \\varphi(\\| \\boldsymbol{x} - \\boldsymbol{x_{C_i}} \\|) where :math:`p(\\boldsymbol{x})` is a low_degree polynomial term, :math:`\\gamma_i` is the weight, corresponding to the a-priori selected :math:`\\mathcal{N}_C` control points, associated to the :math:`i`-th basis function, and :math:`\\varphi(\\| \\boldsymbol{x} - \\boldsymbol{x_{C_i}} \\|)` a radial function based on the Euclidean distance between the control points position :math:`\\boldsymbol{x_{C_i}}` and :math:`\\boldsymbol{x}`. A radial basis function, generally, is a real-valued function whose value depends only on the distance from the origin, so that :math:`\\varphi(\\boldsymbol{x}) = \\tilde{\\varphi}( \\|\\boldsymbol{x} \\|)`. The matrix version of the formula above is: .. math:: \\mathcal{M}(\\boldsymbol{x}) = \\boldsymbol{c} + \\boldsymbol{Q}\\boldsymbol{x} + \\boldsymbol{W^T}\\boldsymbol{d}(\\boldsymbol{x}) The idea is that after the computation of the weights and the polynomial terms from the coordinates of the control points before and after the deformation, we can deform all the points of the mesh accordingly. Among the most common used radial basis functions for modelling 2D and 3D shapes, we consider Gaussian splines, Multi- uadratic biharmonic splines, Inverted multi-quadratic biharmonic splines, Thin-plate splines, Beckert and Wendland :math:`C^2` basis all defined and implemented below. :param string basis: RBF basis function :param float radius: cut-off radius """ def __init__(self, basis, radius): self.bases = { 'gaussian_spline': self.gaussian_spline, 'multi_quadratic_biharmonic_spline': self.multi_quadratic_biharmonic_spline, 'inv_multi_quadratic_biharmonic_spline': self.inv_multi_quadratic_biharmonic_spline, 'thin_plate_spline': self.thin_plate_spline, 'beckert_wendland_c2_basis': self.beckert_wendland_c2_basis } if basis in self.bases: self.basis = self.bases[basis] else: raise NameError( """The name of the basis function is not correct. Check the documentation for all the available functions.""") self.radius = radius # The following static methods are the implementations # of the most common radial basis functions @staticmethod
[docs] def gaussian_spline(X, r): """ It implements the following formula: .. math:: \\varphi(\\| \\boldsymbol{x} \\|) = e^{-\\frac{\\| \\boldsymbol{x} \\|^2}{r^2}} :param numpy.ndarray X: l2-norm between given inputs of a function and the locations to perform rbf approximation to that function. :param float r: smoothing length, also called the cut-off radius. :return: result: the result of the formula above. :rtype: float """ return np.exp(-(X * X) / (r * r))
@staticmethod
[docs] def multi_quadratic_biharmonic_spline(X, r): """ It implements the following formula: .. math:: \\varphi(\\| \\boldsymbol{x} \\|) = \\sqrt{\\| \\boldsymbol{x} \\|^2 + r^2} :param numpy.ndarray X: l2-norm between given inputs of a function and the locations to perform rbf approximation to that function. :param float r: smoothing length, also called the cut-off radius. :return: result: the result of the formula above. :rtype: float """ return np.sqrt((X * X) + (r * r))
@staticmethod
[docs] def inv_multi_quadratic_biharmonic_spline(X, r): """ It implements the following formula: .. math:: \\varphi(\\| \\boldsymbol{x} \\|) = (\\| \\boldsymbol{x} \\|^2 + r^2 )^{-\\frac{1}{2}} :param numpy.ndarray X: l2-norm between given inputs of a function and the locations to perform rbf approximation to that function. :param float r: smoothing length, also called the cut-off radius. :return: result: the result of the formula above. :rtype: float """ return 1.0 / (np.sqrt((X * X) + (r * r)))
@staticmethod
[docs] def thin_plate_spline(X, r): """ It implements the following formula: .. math:: \\varphi(\\| \\boldsymbol{x} \\|) = \\left\\| \\frac{\\boldsymbol{x} }{r} \\right\\|^2 \\ln \\left\\| \\frac{\\boldsymbol{x} }{r} \\right\\| :param numpy.ndarray X: l2-norm between given inputs of a function and the locations to perform rbf approximation to that function. :param float r: smoothing length, also called the cut-off radius. :return: result: the result of the formula above. :rtype: float """ arg = X / r result = arg * arg result = np.where(arg > 0, result * np.log(arg), result) return result
@staticmethod
[docs] def beckert_wendland_c2_basis(X, r): """ It implements the following formula: .. math:: \\varphi(\\| \\boldsymbol{x} \\|) = \\left( 1 - \\frac{\\| \\boldsymbol{x} \\|}{r} \\right)^4 + \\left( 4 \\frac{\\| \\boldsymbol{x} \\|}{r} + 1 \\right) :param numpy.ndarray X: l2-norm between given inputs of a function and the locations to perform rbf approximation to that function. :param float r: smoothing length, also called the cut-off radius. :return: result: the result of the formula above. :rtype: float """ arg = X / r first = np.where((1 - arg) > 0, np.power((1 - arg), 4), 0) second = (4 * arg) + 1 return first * second
[docs] def weights_matrix(self, X1, X2): """ This method returns the following matrix: .. math:: \\boldsymbol{D_{ij}} = \\varphi(\\| \\boldsymbol{x_i} - \\boldsymbol{y_j} \\|) :param numpy.ndarray X1: the vector x in the formula above. :param numpy.ndarray X2: the vector y in the formula above. :return: matrix: the matrix D. :rtype: numpy.ndarray """ if X1.size == X1.shape[0]: X1 = X1.reshape(-1, 1) if X2.size == X2.shape[0]: X2 = X2.reshape(-1, 1) return self.basis(cdist(X1, X2), self.radius)
[docs]def reconstruct_f(original_input, original_output, rbf_input, rbf_output, basis, radius): """ Reconstruct a function by using the radial basis function approximations. :param array_like original_input: the original values of function inputs. :param array_like original_output: the original values of function output. :param array_like rbf_input: the input data for RBF approximation. :param array_like rbf_output: the array elements to be updated with the RBF interpolated outputs after the approximation. :param string basis: radial basis function. :param float radius: smoothing length, also called the cut-off radius. :Example: >>> import numpy as np >>> from bladex.ndinterpolator import reconstruct_f >>> x = np.arange(10) >>> y = np.square(x) >>> radius = 10 >>> n_interp = 50 >>> x_rbf = np.linspace(x[0], x[-1], num=n_interp) >>> y_rbf = np.zeros(n_interp) >>> reconstruct_f(original_input=x, original_output=y, rbf_input=x_rbf, rbf_output=y_rbf, radius=radius, basis='beckert_wendland_c2_basis') """ radial = RBF(basis=basis, radius=radius) weights_coeff = np.linalg.solve( radial.weights_matrix(original_input, original_input), original_output) weights_rbf = radial.weights_matrix(rbf_input, original_input) for i in range(rbf_input.shape[0]): for j in range(original_input.shape[0]): rbf_output[i] += weights_coeff[j] * weights_rbf[i][j]
[docs]def scipy_bspline(cv, npoints, degree): """ Construct B-Spline curve via the control points. :param array_like cv: control points vector (spline coefficients). :param int npoints: number of points on the curve. :param int degree: BSpline degree. """ c = cv.shape[0] # calculating BSpline knots kv = np.clip(np.arange(c+degree+1)-degree, 0, c-degree) # u is the array of points at which to return the value of the smoothed spline u = np.linspace(0, c-degree, npoints) # Calculate resulting spline # splev requires u, and a tuple which is a sequence of length 3 containing # the knots, coefficients, and degree of the spline. return np.array(splev(u, (kv, cv.T, degree))).T