wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-08-12 12:30:37 +08:00
Code Issues Packages Projects Releases Wiki Activity

Make interpolation matrices using cython.

Browse Source
This commit is contained in:
rowanc1
2014-04-16 18:37:36 -07:00
parent 8142966920
commit f8204e7751
6 changed files with 8795 additions and 132 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Tests/test_interpolation.py
+11 -10
View File
@@ -20,6 +20,7 @@ cartE2Cyl = lambda M, ex, ey: cart_row2(M.gridEy, ex, ey)
cartF3 = lambda M, fx, fy, fz: np.vstack((cart_row3(M.gridFx, fx, fy, fz), cart_row3(M.gridFy, fx, fy, fz), cart_row3(M.gridFz, fx, fy, fz)))
cartE3 = lambda M, ex, ey, ez: np.vstack((cart_row3(M.gridEx, ex, ey, ez), cart_row3(M.gridEy, ex, ey, ez), cart_row3(M.gridEz, ex, ey, ez)))
TOL = 1e-7
class TestInterpolation1D(OrderTest):
@@ -64,9 +65,9 @@ class TestOutliersInterp1D(unittest.TestCase):
M = Mesh.TensorMesh([4])
Q = M.getInterpolationMat(np.array([[0],[0.126],[0.127]]),'CC',zerosOutside=True)
x = np.arange(4)+1
self.assertTrue(np.all(Q*x == [1,1.004,1.008]))
self.assertTrue(np.linalg.norm(Q*x - np.r_[1,1.004,1.008]) < TOL)
Q = M.getInterpolationMat(np.array([[-1],[0.126],[0.127]]),'CC',zerosOutside=True)
self.assertTrue(np.all(Q*x == [0,1.004,1.008]))
self.assertTrue(np.linalg.norm(Q*x - np.r_[0,1.004,1.008]) < TOL)
class TestInterpolation2d(OrderTest):
name = "Interpolation 2D"
@@ -243,42 +244,42 @@ class TestInterpolation3D(OrderTest):
def test_orderCC(self):
self.type = 'CC'
self.name = 'Interpolation CC'
self.name = 'Interpolation 3D: CC'
self.orderTest()
def test_orderN(self):
self.type = 'N'
self.name = 'Interpolation N'
self.name = 'Interpolation 3D: N'
self.orderTest()
def test_orderFx(self):
self.type = 'Fx'
self.name = 'Interpolation Fx'
self.name = 'Interpolation 3D: Fx'
self.orderTest()
def test_orderFy(self):
self.type = 'Fy'
self.name = 'Interpolation Fy'
self.name = 'Interpolation 3D: Fy'
self.orderTest()
def test_orderFz(self):
self.type = 'Fz'
self.name = 'Interpolation Fz'
self.name = 'Interpolation 3D: Fz'
self.orderTest()
def test_orderEx(self):
self.type = 'Ex'
self.name = 'Interpolation Ex'
self.name = 'Interpolation 3D: Ex'
self.orderTest()
def test_orderEy(self):
self.type = 'Ey'
self.name = 'Interpolation Ey'
self.name = 'Interpolation 3D: Ey'
self.orderTest()
def test_orderEz(self):
self.type = 'Ez'
self.name = 'Interpolation Ez'
self.name = 'Interpolation 3D: Ez'
self.orderTest()
SimPEG/Utils/interputils.py
+117 -118
View File
@@ -2,34 +2,19 @@ import numpy as np
import scipy.sparse as sp
from matutils import mkvc, sub2ind, spzeros
def _interp_point_1D(x, xr_i):
try:
import interputils_cython as pyx
_interp_point_1D = pyx._interp_point_1D
_interpmat1D = pyx._interpmat1D
_interpmat2D = pyx._interpmat2D
_interpmat3D = pyx._interpmat3D
_interpCython = True
except ImportError, e:
print """Efficiency Warning: Interpolation will be slow, use setup.py!
python setup.py build_ext --inplace
"""
given a point, xr_i, this will find which two integers it lies between.
:param numpy.ndarray x: Tensor vector of 1st dimension of grid.
:param float xr_i: Location of a point
:rtype: int,int,float,float
:return: index1, index2, portion1, portion2
"""
im = np.argmin(abs(x-xr_i))
if xr_i - x[im] >= 0: # Point on the left
ind_x1 = im
ind_x2 = im+1
elif xr_i - x[im] < 0: # Point on the right
ind_x1 = im-1
ind_x2 = im
ind_x1 = max(min(ind_x1, x.size-1), 0)
ind_x2 = max(min(ind_x2, x.size-1), 0)
if ind_x1 == ind_x2:
return ind_x1, ind_x1, 0.5, 0.5
hx = x[ind_x2] - x[ind_x1]
wx1 = 1 - (xr_i - x[ind_x1])/hx
wx2 = 1 - (x[ind_x2] - xr_i)/hx
return ind_x1, ind_x2, wx1, wx2
_interpCython = False
def interpmat(locs, x, y=None, z=None):
@@ -60,118 +45,132 @@ def interpmat(locs, x, y=None, z=None):
plt.show()
"""
npts = locs.shape[0]
locs = locs.astype(float)
x = x.astype(float)
if y is None and z is None:
return _interpmat1D(locs, x)
shape = [x.size,]
inds, vals = _interpmat1D(mkvc(locs), x)
elif z is None:
return _interpmat2D(locs, x, y)
y = y.astype(float)
shape = [x.size, y.size]
inds, vals = _interpmat2D(locs, x, y)
else:
return _interpmat3D(locs, x, y, z)
def _interpmat1D(locs, x):
"""Use interpmat with only x component provided."""
nx = x.size
locs = mkvc(locs)
npts = locs.shape[0]
Q = sp.lil_matrix((npts, nx))
for i in range(npts):
ind_x1, ind_x2, wx1, wx2 = _interp_point_1D(x, locs[i])
# dv = (x[ind_x2] - x[ind_x1])
# Get the row in the matrix
inds = [ind_x1, ind_x2]
vals = [wx1,wx2]
for I, V in zip(inds, vals):
Q[i, I] += V
# Q[i, mkvc(inds)] = vals
y = y.astype(float)
z = z.astype(float)
shape = [x.size, y.size, z.size]
inds, vals = _interpmat3D(locs, x, y, z)
I = np.repeat(range(npts),2**len(shape))
J = sub2ind(shape,inds)
Q = sp.csr_matrix((vals,(I, J)),
shape=(npts, np.prod(shape)))
return Q
if not _interpCython:
def _interp_point_1D(x, xr_i):
"""
given a point, xr_i, this will find which two integers it lies between.
:param numpy.ndarray x: Tensor vector of 1st dimension of grid.
:param float xr_i: Location of a point
:rtype: int,int,float,float
:return: index1, index2, portion1, portion2
"""
im = np.argmin(abs(x-xr_i))
if xr_i - x[im] >= 0: # Point on the left
ind_x1 = im
ind_x2 = im+1
elif xr_i - x[im] < 0: # Point on the right
ind_x1 = im-1
ind_x2 = im
ind_x1 = max(min(ind_x1, x.size-1), 0)
ind_x2 = max(min(ind_x2, x.size-1), 0)
if ind_x1 == ind_x2:
return ind_x1, ind_x1, 0.5, 0.5
hx = x[ind_x2] - x[ind_x1]
wx1 = 1 - (xr_i - x[ind_x1])/hx
wx2 = 1 - (x[ind_x2] - xr_i)/hx
return ind_x1, ind_x2, wx1, wx2
def _interpmat1D(locs, x):
"""Use interpmat with only x component provided."""
nx = x.size
npts = locs.shape[0]
inds, vals = [], []
for i in range(npts):
ind_x1, ind_x2, wx1, wx2 = _interp_point_1D(x, locs[i])
inds += [ind_x1, ind_x2]
vals += [wx1,wx2]
return inds, vals
def _interpmat2D(locs, x, y):
"""Use interpmat with only x and y components provided."""
nx = x.size
ny = y.size
npts = locs.shape[0]
def _interpmat2D(locs, x, y):
"""Use interpmat with only x and y components provided."""
nx = x.size
ny = y.size
npts = locs.shape[0]
Q = sp.lil_matrix((npts, nx*ny))
inds, vals = [], []
for i in range(npts):
ind_x1, ind_x2, wx1, wx2 = _interp_point_1D(x, locs[i, 0])
ind_y1, ind_y2, wy1, wy2 = _interp_point_1D(y, locs[i, 1])
for i in range(npts):
ind_x1, ind_x2, wx1, wx2 = _interp_point_1D(x, locs[i, 0])
ind_y1, ind_y2, wy1, wy2 = _interp_point_1D(y, locs[i, 1])
inds += [( ind_x1, ind_y2),
( ind_x1, ind_y1),
( ind_x2, ind_y1),
( ind_x2, ind_y2)]
# dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1])
vals += [wx1*wy2,
wx1*wy1,
wx2*wy1,
wx2*wy2]
# Get the row in the matrix
inds = sub2ind((nx,ny),[
( ind_x1, ind_y2),
( ind_x1, ind_y1),
( ind_x2, ind_y1),
( ind_x2, ind_y2)])
vals = [wx1*wy2,
wx1*wy1,
wx2*wy1,
wx2*wy2]
for I, V in zip(mkvc(inds), vals):
Q[i, I] += V
# Q[i, mkvc(inds)] = vals
return Q
return inds, vals
def _interpmat3D(locs, x, y, z):
"""Use interpmat."""
nx = x.size
ny = y.size
nz = z.size
npts = locs.shape[0]
def _interpmat3D(locs, x, y, z):
"""Use interpmat."""
nx = x.size
ny = y.size
nz = z.size
npts = locs.shape[0]
Q = sp.lil_matrix((npts, nx*ny*nz))
inds, vals = [], []
for i in range(npts):
ind_x1, ind_x2, wx1, wx2 = _interp_point_1D(x, locs[i, 0])
ind_y1, ind_y2, wy1, wy2 = _interp_point_1D(y, locs[i, 1])
ind_z1, ind_z2, wz1, wz2 = _interp_point_1D(z, locs[i, 2])
for i in range(npts):
ind_x1, ind_x2, wx1, wx2 = _interp_point_1D(x, locs[i, 0])
ind_y1, ind_y2, wy1, wy2 = _interp_point_1D(y, locs[i, 1])
ind_z1, ind_z2, wz1, wz2 = _interp_point_1D(z, locs[i, 2])
inds += [( ind_x1, ind_y2, ind_z1),
( ind_x1, ind_y1, ind_z1),
( ind_x2, ind_y1, ind_z1),
( ind_x2, ind_y2, ind_z1),
( ind_x1, ind_y1, ind_z2),
( ind_x1, ind_y2, ind_z2),
( ind_x2, ind_y1, ind_z2),
( ind_x2, ind_y2, ind_z2)]
# dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1]) *(z[ind_z2] - z[ind_z1])
vals += [wx1*wy2*wz1,
wx1*wy1*wz1,
wx2*wy1*wz1,
wx2*wy2*wz1,
wx1*wy1*wz2,
wx1*wy2*wz2,
wx2*wy1*wz2,
wx2*wy2*wz2]
# Get the row in the matrix
inds = sub2ind((nx,ny,nz),[
( ind_x1, ind_y2, ind_z1),
( ind_x1, ind_y1, ind_z1),
( ind_x2, ind_y1, ind_z1),
( ind_x2, ind_y2, ind_z1),
( ind_x1, ind_y1, ind_z2),
( ind_x1, ind_y2, ind_z2),
( ind_x2, ind_y1, ind_z2),
( ind_x2, ind_y2, ind_z2)])
vals = [wx1*wy2*wz1,
wx1*wy1*wz1,
wx2*wy1*wz1,
wx2*wy2*wz1,
wx1*wy1*wz2,
wx1*wy2*wz2,
wx2*wy1*wz2,
wx2*wy2*wz2]
for I, V in zip(mkvc(inds), vals):
Q[i, I] += V
# Q[i, mkvc(inds)] = vals
return Q
return inds, vals
if __name__ == '__main__':
SimPEG/Utils/interputils_cython.c
+8535
View File
File diff suppressed because it is too large Load Diff
SimPEG/Utils/interputils_cython.pyx
+119
View File
@@ -0,0 +1,119 @@
# from __future__ import division
import numpy as np
cimport numpy as np
# from libcpp.vector cimport vector
def _interp_point_1D (np.ndarray[np.float64_t, ndim=1] x, float xr_i):
"""
given a point, xr_i, this will find which two integers it lies between.
:param numpy.ndarray x: Tensor vector of 1st dimension of grid.
:param float xr_i: Location of a point
:rtype: int,int,float,float
:return: index1, index2, portion1, portion2
"""
# TODO: This fails if the point is on the outside of the mesh.
# We may want to replace this by extrapolation?
cdef int im = np.argmin(abs(x-xr_i))
cdef int ind_x1 = 0
cdef int ind_x2 = 0
cdef int xSize = x.shape[0]-1
cdef float wx1 = 0.0
cdef float wx2 = 0.0
cdef float hx = 0.0
if xr_i - x[im] >= 0: # Point on the left
ind_x1 = im
ind_x2 = im+1
elif xr_i - x[im] < 0: # Point on the right
ind_x1 = im-1
ind_x2 = im
ind_x1 = max(min(ind_x1, xSize), 0)
ind_x2 = max(min(ind_x2, xSize), 0)
if ind_x1 == ind_x2:
return ind_x1, ind_x1, 0.5, 0.5
hx = x[ind_x2] - x[ind_x1]
wx1 = 1 - (xr_i - x[ind_x1])/hx
wx2 = 1 - (x[ind_x2] - xr_i)/hx
return ind_x1, ind_x2, wx1, wx2
def _interpmat1D(np.ndarray[np.float64_t, ndim=1] locs,
np.ndarray[np.float64_t, ndim=1] x):
"""Use interpmat with only x component provided."""
cdef int nx = x.size
cdef int npts = locs.shape[0]
inds, vals = [], []
for i in range(npts):
ind_x1, ind_x2, wx1, wx2 = _interp_point_1D(x, locs[i])
inds += [ind_x1, ind_x2]
vals += [wx1,wx2]
return inds, vals
def _interpmat2D(np.ndarray[np.float64_t, ndim=2] locs,
np.ndarray[np.float64_t, ndim=1] x,
np.ndarray[np.float64_t, ndim=1] y):
"""Use interpmat with only x and y components provided."""
cdef int nx = x.size
cdef int ny = y.size
cdef int npts = locs.shape[0]
inds, vals = [], []
for i in range(npts):
ind_x1, ind_x2, wx1, wx2 = _interp_point_1D(x, locs[i, 0])
ind_y1, ind_y2, wy1, wy2 = _interp_point_1D(y, locs[i, 1])
inds += [( ind_x1, ind_y2),
( ind_x1, ind_y1),
( ind_x2, ind_y1),
( ind_x2, ind_y2)]
vals += [wx1*wy2, wx1*wy1, wx2*wy1, wx2*wy2]
return inds, vals
def _interpmat3D(np.ndarray[np.float64_t, ndim=2] locs,
np.ndarray[np.float64_t, ndim=1] x,
np.ndarray[np.float64_t, ndim=1] y,
np.ndarray[np.float64_t, ndim=1] z):
"""Use interpmat."""
cdef int nx = x.size
cdef int ny = y.size
cdef int nz = z.size
cdef int npts = locs.shape[0]
inds, vals = [], []
for i in range(npts):
ind_x1, ind_x2, wx1, wx2 = _interp_point_1D(x, locs[i, 0])
ind_y1, ind_y2, wy1, wy2 = _interp_point_1D(y, locs[i, 1])
ind_z1, ind_z2, wz1, wz2 = _interp_point_1D(z, locs[i, 2])
inds += [( ind_x1, ind_y2, ind_z1),
( ind_x1, ind_y1, ind_z1),
( ind_x2, ind_y1, ind_z1),
( ind_x2, ind_y2, ind_z1),
( ind_x1, ind_y1, ind_z2),
( ind_x1, ind_y2, ind_z2),
( ind_x2, ind_y1, ind_z2),
( ind_x2, ind_y2, ind_z2)]
vals += [wx1*wy2*wz1,
wx1*wy1*wz1,
wx2*wy1*wz1,
wx2*wy2*wz1,
wx1*wy1*wz2,
wx1*wy2*wz2,
wx2*wy1*wz2,
wx2*wy2*wz2]
return inds, vals
SimPEG/Utils/matutils.py
+2
View File
@@ -156,6 +156,8 @@ def ind2sub(shape, inds):
def sub2ind(shape, subs):
"""From the given shape, returns the index of the given subscript"""
if len(shape) == 1:
return subs
if type(subs) is not np.ndarray:
subs = np.array(subs)
if len(subs.shape) == 1:
setup.py
+11 -4
View File
@@ -5,9 +5,14 @@ SimPEG is a python package for simulation and gradient based
parameter estimation in the context of geophysical applications.
"""
import ez_setup
ez_setup.use_setuptools()
from setuptools import setup, find_packages
# import ez_setup
# ez_setup.use_setuptools()
# from setuptools import setup, find_packages
from distutils.core import setup
from setuptools import find_packages
from Cython.Build import cythonize
import numpy as np
CLASSIFIERS = [
'Development Status :: 0.0.1 - Alpha',
@@ -42,5 +47,7 @@ setup(
download_url = "http://github.com/simpeg",
classifiers=CLASSIFIERS,
platforms = ["Windows", "Linux", "Solaris", "Mac OS-X", "Unix"],
use_2to3 = False
use_2to3 = False,
include_dirs=[np.get_include()],
ext_modules = cythonize('SimPEG/Utils/interputils_cython.pyx')
)