Source code for compmech.interpolate

r"""
Interpolate (:mod:`compmech.interpolate`)
==================================================

.. currentmodule:: compmech.interpolate

This module includes some interpolation utilities that will be used in other
modules.

.. autofunction:: inv_weighted

.. autofunction:: interp

"""
from collections import Iterable

import numpy as np

from compmech.logger import msg, warn


[docs]def inv_weighted(data, mesh, num_sub, col, ncp=5, power_parameter=2): r"""Interpolates the values taken at one group of points into another using an inverse-weighted algorithm In the inverse-weighted algorithm a number of `n_{CP}` measured points of the input parameter ``data`` that are closest to a given node in the input parameter ``mesh`` are found and the imperfection value of this node (represented by the normal displacement `{w_0}_{node}`) is calculated as follows: .. math:: {w_0}_{node} = \frac{\sum_{i}^{n_{CP}}{{w_0}_i\frac{1}{w_i}}} {\sum_{i}^{n_{CP}}{\frac{1}{w_i}}} where `w_i` is the imperfection at each measured point, calculated as: .. math:: w_i = \left[(x_{node}-x_i)^2+(y_{node}-y_i)^2+(z_{node}-y_i)^2 \right]^p with `p` being a power parameter that when increased will increase the relative influence of a closest point. Parameters ---------- data : numpy.ndarray, shape (N, ndim+1) The data or an array containing the imperfection file. The values to be interpolated must be in the last column. mesh : numpy.ndarray, shape (M, ndim) The new coordinates where the values will be interpolated to. num_sub : int The number of sub-sets used during the interpolation. The points are divided in sub-sets to increase the algorithm's efficiency. col : int The index of the column to be used in order to divide the data in sub-sets. Note that the first column index is ``0``. ncp : int, optional Number of closest points used in the inverse-weighted interpolation. power_parameter : float, optional Power of inverse weighted interpolation function. Returns ------- ans : numpy.ndarray A 1-D array with the interpolated values. The size of this array is ``mesh.shape[0]``. """ if mesh.shape[1] != data.shape[1]-1: raise ValueError('Invalid input: mesh.shape[1] != data.shape[1]') msg('Interpolating... ') num_sub = int(num_sub) mesh_size = mesh.shape[0] # memory control mem_limit = 1024*1024*1024*8*2 # 2 GB mem_entries = int(mem_limit / 64) # if float64 is used sec_size = int(mesh_size/num_sub) while sec_size**2*10 > mem_entries: num_sub +=1 sec_size = int(mesh_size/num_sub) if sec_size**2*10 <= mem_entries: warn('New num_sub: {0}'.format(int(mesh_size/float(sec_size)))) break mesh_seq = np.arange(mesh.shape[0]) mesh_argsort = np.argsort(mesh[:, col]) mesh_seq = mesh_seq[mesh_argsort] back_argsort = np.argsort(mesh_seq) mesh = np.asarray(mesh[mesh_argsort], order='F') length = mesh[:, col].max() - mesh[:, col].min() data = np.asarray(data[np.argsort(data[:, col])], order='F') ans = np.zeros(mesh.shape[0], dtype=mesh.dtype) # max_num_limits defines how many times the log will print # "processed ... out of ... entries" max_num_limits = 10 for den in range(max_num_limits, 0, -1): if num_sub % den == 0: limit = int(num_sub/den) break for i in range(num_sub+1): i_inf = sec_size*i i_sup = sec_size*(i+1) if i % limit == 0: msg('\t processed {0:7d} out of {1:7d} entries'.format( min(i_sup, mesh_size), mesh_size)) sub_mesh = mesh[i_inf : i_sup] if not np.any(sub_mesh): continue inf = sub_mesh[:, col].min() sup = sub_mesh[:, col].max() tol = 0.03 if i == 0 or i == num_sub: tol = 0.06 while True: cond1 = data[:, col] >= inf - tol*length cond2 = data[:, col] <= sup + tol*length cond = np.all(np.array((cond1, cond2)), axis=0) sub_data = data[cond] if not np.any(sub_data): tol += 0.01 else: break dist = np.subtract.outer(sub_mesh[:, 0], sub_data[:, 0])**2 for j in range(1, sub_mesh.shape[1]): dist += np.subtract.outer(sub_mesh[:, j], sub_data[:, j])**2 asort = np.argsort(dist, axis=1) lenn = sub_mesh.shape[0] lenp = sub_data.shape[0] asort_mesh = asort + np.meshgrid(np.arange(lenn)*lenp, np.arange(lenp))[0].transpose() # getting the distance of the closest points dist_cp = np.take(dist, asort_mesh[:, :ncp]) # avoiding division by zero dist_cp[(dist_cp==0)] == 1.e-12 # fetching the imperfection of the sub-data imp = sub_data[:, -1] # taking only the imperfection of the closest points imp_cp = np.take(imp, asort[:, :ncp]) # weight calculation total_weight = np.sum(1./(dist_cp**power_parameter), axis=1) weight = 1./(dist_cp**power_parameter) # computing the new imp imp_new = np.sum(imp_cp*weight, axis=1)/total_weight # updating the answer array ans[i_inf : i_sup] = imp_new ans = ans[back_argsort] msg('Interpolation completed!') return ans
[docs]def interp(x, xp, fp, left=None, right=None, period=None): """ One-dimensional linear interpolation Returns the one-dimensional piecewise linear interpolant to a function with given values at discrete data-points. .. note:: This function has been incorporated in NumPy >= 1.10.0 and will be soon removed from here. Parameters ---------- x : array_like The x-coordinates of the interpolated values. xp : 1-D sequence of floats The x-coordinates of the data points, must be increasing if argument ``period`` is not specified. Otherwise, ``xp`` is internally sorted after normalizing the periodic boundaries with ``xp = xp % period``. fp : 1-D sequence of floats The y-coordinates of the data points, same length as ``xp``. left : float, optional Value to return for ``x < xp[0]``, default is ``fp[0]``. right : float, optional Value to return for ``x > xp[-1]``, default is ``fp[-1]``. period : float, optional A period for the x-coordinates. This parameter allows the proper interpolation of angular x-coordinates. Parameters ``left`` and ``right`` are ignored if ``period`` is specified. Returns ------- y : {float, ndarray} The interpolated values, same shape as ``x``. Raises ------ ValueError If ``xp`` and ``fp`` have different length If ``xp`` or ``fp`` are not 1-D sequences If ``period==0`` Notes ----- Does not check that the x-coordinate sequence ``xp`` is increasing. If ``xp`` is not increasing, the results are nonsense. A simple check for increasing is:: np.all(np.diff(xp) > 0) Examples -------- >>> xp = [1, 2, 3] >>> fp = [3, 2, 0] >>> interp(2.5, xp, fp) 1.0 >>> interp([0, 1, 1.5, 2.72, 3.14], xp, fp) array([ 3. , 3. , 2.5 , 0.56, 0. ]) >>> UNDEF = -99.0 >>> interp(3.14, xp, fp, right=UNDEF) -99.0 Plot an interpolant to the sine function: >>> x = np.linspace(0, 2*np.pi, 10) >>> y = np.sin(x) >>> xvals = np.linspace(0, 2*np.pi, 50) >>> yinterp = interp(xvals, x, y) >>> import matplotlib.pyplot as plt >>> plt.plot(x, y, 'o') [<matplotlib.lines.Line2D object at 0x...>] >>> plt.plot(xvals, yinterp, '-x') [<matplotlib.lines.Line2D object at 0x...>] >>> plt.show() Interpolation with periodic x-coordinates: >>> x = [-180, -170, -185, 185, -10, -5, 0, 365] >>> xp = [190, -190, 350, -350] >>> fp = [5, 10, 3, 4] >>> interp(x, xp, fp, period=360) array([7.5, 5., 8.75, 6.25, 3., 3.25, 3.5, 3.75]) """ if period is None: return np.interp(x, xp, fp, left, right) else: if period==0: raise ValueError('Argument `period` must be a non-zero value') period = abs(period) if not isinstance(x, Iterable): x = [x] x = np.asarray(x) xp = np.asarray(xp) fp = np.asarray(fp) if xp.ndim != 1 or fp.ndim != 1: raise ValueError('Data points must be 1-D sequences') if xp.shape[0] != fp.shape[0]: raise ValueError('Inputs `xp` and `fp` must have the same shape') # eliminating discontinuity between periods x = x % period xp = xp % period asort_xp = np.argsort(xp) xp = xp[asort_xp] fp = fp[asort_xp] xp = np.concatenate((xp[-1:]-period, xp, xp[0:1]+period)) fp = np.concatenate((fp[-1:], fp, fp[0:1])) return np.interp(x, xp, fp)
if __name__=='__main__': a = np.array([[1.1, 1.2, 10], [1.2, 1.2, 10], [1.3, 1.3, 10], [1.4, 1.3, 10], [1.5, 1.3, 10], [2.6, 2.3, 5], [2.7, 2.3, 5], [2.6, 2.1, 5], [2.7, 2.1, 5], [2.8, 2.2, 5], [2.8, 2.2, 5], [5.6, 5.3, 20], [5.7, 5.3, 20], [5.6, 5.1, 20], [5.7, 5.1, 20], [5.8, 5.2, 20], [5.8, 5.2, 20]]) b = np.array([[1., 1.], [2., 2.], [4., 4.], [5., 5.]]) print inv_weighted(a, b, num_sub=1, col=1, ncp=10)