import numpy
import scipy.linalg
[docs]class LinearInterpolationNOGrid:
"""
The class takes a function on a uniform, rectangular grid (with orthogonal or
non-orthogonal lattice vectors) and calculates the function at any point
using linear interpolation. It can be a 1D, 2D or 3D grid.
See http://en.wikipedia.org/wiki/Trilinear_interpolation for the math.
Usage::
func=[[0,1],[2,3]]
start=[0,0]
latticevecs=[[1e-3,0],[0,1e-3]]
interp=LinearInterpolationNOGrid(func,latticevecs)
f=interp(0.3,0.6)
Note: the gridpoints on the upper boundary in each direction
are not accessible by the interpolation function.
TODO: fix that.
Also note: Pay attention to the ordering of the ``func`` array.
"""
#for higher performance, implement bilinear and linear
#interpolation, and extra functions for orthogonal grid
__dim=None
__func=None
__start=None
__latticevecs=None
__transformation_matrix=None
__default=None
def __init__(self,func,grid=1,start=None,default=None):
"""
func: multi-dimensional array containing the function
values (list or numpy.array) in a right-handed coordinate system
(e.g. x to the left, y up), with the indices in the same order
as the grid sizes/lattice vectors (i.e. if the lattice vectors
point in x,y and z direction, func will be called like
func[x,y,z])
start: the start point of the coordinate system given by latticevecs
(list or numpy.array). Default is None, which means that
the grid is starting at the coordinate origin.
grid: either 1) Lattice vectors (list or numpy.array)
or 2) a number, which is the gridsize.
Default is 1.
default: The function value of points outside the area.
Default is None.
"""
#Convert lists to numpy.array
self.__func=numpy.array(func,copy=False)
self.__dim=len(self.__func.shape)
self.__default=default
if start==None:
self.__start=numpy.zeros(self.__dim)
else:
self.__start=numpy.array(start,copy=False)
if isinstance(grid,(int,long,float)): #http://stackoverflow.com/a/3501408/1447622
basisvecs=grid*numpy.eye(self.__dim)
else:
basisvecs=numpy.array(grid,copy=False)
self.__transformation_matrix=self.__calc_transformation_matrix(basisvecs)
def __call__(self,*point):
"""
The interpolation function can be called like::
my_interp(1,2,3)
as well as::
my_interp([1,2,3])
"""
if(len(point)==1):
point=point[0]
grid_element,relative_coordinates=self.__find_position_in_no_grid(point)
f000=f001=f010=f011=f100=f101=f110=f111=0
x=y=z=0
try:
if(self.__dim==1):
f000=self.__func[tuple(grid_element+[0])]
f100=self.__func[tuple(grid_element+[1])]
x,=relative_coordinates
if(self.__dim==2):
f000=self.__func[tuple(grid_element+[0,0])]
f010=self.__func[tuple(grid_element+[0,1])]
f100=self.__func[tuple(grid_element+[1,0])]
f110=self.__func[tuple(grid_element+[1,1])]
x,y=relative_coordinates
if(self.__dim==3):
f000=self.__func[tuple(grid_element+[0,0,0])]
f001=self.__func[tuple(grid_element+[0,0,1])]
f010=self.__func[tuple(grid_element+[0,1,0])]
f011=self.__func[tuple(grid_element+[0,1,1])]
f100=self.__func[tuple(grid_element+[1,0,0])]
f101=self.__func[tuple(grid_element+[1,0,1])]
f110=self.__func[tuple(grid_element+[1,1,0])]
f111=self.__func[tuple(grid_element+[1,1,1])]
x,y,z=relative_coordinates
return self.__trilinear_interpolation(f000, f001, f010, f011,
f100, f101, f110, f111,
x, y, z)
except IndexError:
return self.__default
def __find_position_in_no_grid(self,point):
"""
Finds the grid element the point is in and its position in this
element in terms of the NO basis.
point: The point (numpy.array)
"""
point_in_no_basis=numpy.dot(self.__transformation_matrix,
(point-self.__start))
grid_element=numpy.floor(point_in_no_basis)
relative_coordinates=point_in_no_basis-grid_element
return grid_element,relative_coordinates
def __calc_transformation_matrix(self,basisvecs):
"""
Calculate transformation matrix to find corresponding vector in the
NO basis for a vector in cartesian coordinates.
"""
return scipy.linalg.inv(numpy.transpose(basisvecs))
def __trilinear_interpolation(self,f000, f001, f010, f011, f100, f101,
f110, f111, x, y, z):
return (f111*x*y*z +
f110*x*y*(1 - z) +
f101*x*(1 - y)*z +
f100*x*(1 - y)*(1 - z) +
f011*(1 - x)*y*z +
f010*(1 - x)*y*(1 - z) +
f001*(1 - x)*(1 - y)*z +
f000*(1 - x)*(1 - y)*(1 - z))
[docs] def map_function_to_points(self,points):
"""
Return list of function values at given points.
points: list of points
"""
return [self.__call__(point) for point in points]
[docs] def dim(self):
"""
Return the dimension of the interpolation function.
"""
return self.__dim