wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-08-03 13:10:16 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'master' of https://github.com/simpeg/simpeg into cylClean

Browse Source 
Conflicts:
	SimPEG/Mesh/View.py
This commit is contained in:
rowanc1
2014-03-17 09:20:36 -07:00
parent 5bcdc5d46a 0c7eb57d46
commit 8276759c3e
12 changed files with 401 additions and 377 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Inversion.py
-65
View File
@@ -32,24 +32,6 @@ class BaseInversion(object):
self.opt.printers.insert(2,IterationPrinters.phi_d)
self.opt.printers.insert(3,IterationPrinters.phi_m)
#TODO: Move this to the data class?
@property
def phi_d_target(self):
"""
target for phi_d
By default this is the number of data.
Note that we do not set the target if it is None, but we return the default value.
"""
if getattr(self, '_phi_d_target', None) is None:
return self.data.dobs.size #
return self._phi_d_target
@phi_d_target.setter
def phi_d_target(self, value):
self._phi_d_target = value
@Utils.timeIt
def run(self, m0):
"""run(m0)
@@ -70,50 +52,3 @@ class BaseInversion(object):
**finish** is called at the end of the optimization.
"""
pass
def save(self, group):
group.attrs['phi_d'] = self.phi_d
group.attrs['phi_m'] = self.phi_m
group.setArray('m', self.m)
group.setArray('dpred', self.dpred)
# class Inversion(Cooling, Remember, BaseInversion):
# maxIter = 10
# name = "SimPEG Inversion"
# def __init__(self, prob, reg, opt, data, **kwargs):
# BaseInversion.__init__(self, prob, reg, opt, data, **kwargs)
# self.stoppers.append(StoppingCriteria.phi_d_target_Inversion)
# if StoppingCriteria.phi_d_target_Minimize not in self.opt.stoppers:
# self.opt.stoppers.append(StoppingCriteria.phi_d_target_Minimize)
# class TimeSteppingInversion(Remember, BaseInversion):
# """
# A slightly different view on regularization parameters,
# let Beta be viewed as 1/dt, and timestep by updating the
# reference model every optimization iteration.
# """
# maxIter = 1
# name = "Time-Stepping SimPEG Inversion"
# def __init__(self, prob, reg, opt, data, **kwargs):
# BaseInversion.__init__(self, prob, reg, opt, data, **kwargs)
# self.stoppers.append(StoppingCriteria.phi_d_target_Inversion)
# if StoppingCriteria.phi_d_target_Minimize not in self.opt.stoppers:
# self.opt.stoppers.append(StoppingCriteria.phi_d_target_Minimize)
# def _startup_TimeSteppingInversion(self, m0):
# def _doEndIteration_updateMref(self, xt):
# if self.debug: 'Updating the reference model.'
# self.parent.reg.mref = self.xc
# self.opt.hook(_doEndIteration_updateMref, overwrite=True)
SimPEG/Mesh/View.py
+228 -224
View File
@@ -14,229 +14,153 @@ class TensorView(object):
def __init__(self):
pass
def plotImage(self, I, imageType='CC', figNum=1,ax=None,direction='z',numbering=True,annotationColor='w',showIt=False,clim=None):
# def components(self):
# plotAll = len(imageType) == 1
# options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":False}
# fig = plt.figure(figNum)
# # Determine the subplot number: 131, 121
# numPlots = 130 if plotAll else len(imageType)/2*10+100
# pltNum = 1
# fxyz = self.r(I,'F','F','M')
# if plotAll or 'Fx' in imageType:
# ax_x = plt.subplot(numPlots+pltNum)
# self.plotImage(fxyz[0], imageType='Fx', ax=ax_x, **options)
# pltNum +=1
# if plotAll or 'Fy' in imageType:
# ax_y = plt.subplot(numPlots+pltNum)
# self.plotImage(fxyz[1], imageType='Fy', ax=ax_y, **options)
# pltNum +=1
# if plotAll or 'Fz' in imageType:
# ax_z = plt.subplot(numPlots+pltNum)
# self.plotImage(fxyz[2], imageType='Fz', ax=ax_z, **options)
# pltNum +=1
# if showIt: plt.show()
def plotImage(self, v, vType='CC', grid=False, view='real',
ax=None, clim=None, showIt=False,
pcolorOpts={},
streamOpts={'color':'k'},
gridOpts={'color':'k'},
numbering=True, annotationColor='w'
):
"""
Mesh.plotImage(I)
Mesh.plotImage(v)
Plots scalar fields on the given mesh.
Input:
:param numpy.array I: scalar field
:param numpy.array v: vector
Optional Input:
Optional Inputs:
:param str imageType: type of image ('CC','N','F','Fx','Fy','Fz','E','Ex','Ey','Ez') or combinations, e.g. ExEy or FxFz
:param int figNum: number of figure to plot to
:param str vType: type of vector ('CC','N','F','Fx','Fy','Fz','E','Ex','Ey','Ez')
:param matplotlib.axes.Axes ax: axis to plot to
:param str direction: slice dimensions, 3D only ('x', 'y', 'z')
:param bool showIt: call plt.show()
3D Inputs:
:param bool numbering: show numbering of slices, 3D only
:param str annotationColor: color of annotation, e.g. 'w', 'k', 'b'
:param bool showIt: call plt.show()
.. plot::
:include-source:
from SimPEG import Mesh, np
M = Mesh.TensorMesh([20, 20])
I = np.sin(M.gridCC[:,0]*2*np.pi)*np.sin(M.gridCC[:,1]*2*np.pi)
M.plotImage(I, showIt=True)
v = np.sin(M.gridCC[:,0]*2*np.pi)*np.sin(M.gridCC[:,1]*2*np.pi)
M.plotImage(v, showIt=True)
.. plot::
:include-source:
from SimPEG import Mesh, np
M = Mesh.TensorMesh([20,20,20])
I = np.sin(M.gridCC[:,0]*2*np.pi)*np.sin(M.gridCC[:,1]*2*np.pi)*np.sin(M.gridCC[:,2]*2*np.pi)
M.plotImage(I, annotationColor='k', showIt=True)
v = np.sin(M.gridCC[:,0]*2*np.pi)*np.sin(M.gridCC[:,1]*2*np.pi)*np.sin(M.gridCC[:,2]*2*np.pi)
M.plotImage(v, annotationColor='k', showIt=True)
"""
assert type(I) == np.ndarray, "I must be a numpy array"
assert type(numbering) == bool, "numbering must be a bool"
assert direction in ["x", "y","z"], "direction must be either x,y, or z"
if imageType == 'CC':
assert I.size == self.nC, "Incorrect dimensions for CC."
elif imageType == 'N':
assert I.size == self.nN, "Incorrect dimensions for N."
elif imageType == 'Fx':
if I.size != np.prod(self.vnFx): I, fy, fz = self.r(I,'F','F','M')
elif imageType == 'Fy':
if I.size != np.prod(self.vnFy): fx, I, fz = self.r(I,'F','F','M')
elif imageType == 'Fz':
if I.size != np.prod(self.vnFz): fx, fy, I = self.r(I,'F','F','M')
elif imageType == 'Ex':
if I.size != np.prod(self.vnEx): I, ey, ez = self.r(I,'E','E','M')
elif imageType == 'Ey':
if I.size != np.prod(self.vnEy): ex, I, ez = self.r(I,'E','E','M')
elif imageType == 'Ez':
if I.size != np.prod(self.vnEz): ex, ey, I = self.r(I,'E','E','M')
elif imageType[0] == 'E':
plotAll = len(imageType) == 1
options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":False}
fig = plt.figure(figNum)
# Determine the subplot number: 131, 121
numPlots = 130 if plotAll else len(imageType)/2*10+100
pltNum = 1
ex, ey, ez = self.r(I,'E','E','M')
if plotAll or 'Ex' in imageType:
ax_x = plt.subplot(numPlots+pltNum)
self.plotImage(ex, imageType='Ex', ax=ax_x, **options)
pltNum +=1
if plotAll or 'Ey' in imageType:
ax_y = plt.subplot(numPlots+pltNum)
self.plotImage(ey, imageType='Ey', ax=ax_y, **options)
pltNum +=1
if plotAll or 'Ez' in imageType:
ax_z = plt.subplot(numPlots+pltNum)
self.plotImage(ez, imageType='Ez', ax=ax_z, **options)
pltNum +=1
if showIt: plt.show()
return
elif imageType[0] == 'F':
plotAll = len(imageType) == 1
options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":False}
fig = plt.figure(figNum)
# Determine the subplot number: 131, 121
numPlots = 130 if plotAll else len(imageType)/2*10+100
pltNum = 1
fxyz = self.r(I,'F','F','M')
if plotAll or 'Fx' in imageType:
ax_x = plt.subplot(numPlots+pltNum)
self.plotImage(fxyz[0], imageType='Fx', ax=ax_x, **options)
pltNum +=1
if plotAll or 'Fy' in imageType:
ax_y = plt.subplot(numPlots+pltNum)
self.plotImage(fxyz[1], imageType='Fy', ax=ax_y, **options)
pltNum +=1
if plotAll or 'Fz' in imageType:
ax_z = plt.subplot(numPlots+pltNum)
self.plotImage(fxyz[2], imageType='Fz', ax=ax_z, **options)
pltNum +=1
if showIt: plt.show()
return
else:
raise Exception("imageType must be 'CC', 'N','Fx','Fy','Fz','Ex','Ey','Ez'")
if ax is None:
fig = plt.figure(figNum)
fig.clf()
fig = plt.figure()
ax = plt.subplot(111)
else:
assert isinstance(ax,matplotlib.axes.Axes), "ax must be an Axes!"
fig = ax.figure
if self.dim == 1:
if imageType == 'CC':
ph = ax.plot(self.vectorCCx, I, '-ro')
elif imageType == 'N':
ph = ax.plot(self.vectorNx, I, '-bs')
if vType == 'CC':
ph = ax.plot(self.vectorCCx, v, '-ro')
elif vType == 'N':
ph = ax.plot(self.vectorNx, v, '-bs')
ax.set_xlabel("x")
ax.axis('tight')
elif self.dim == 2:
if imageType == 'CC':
C = I[:].reshape(self.vnC, order='F')
elif imageType == 'N':
C = I[:].reshape(self.vnN, order='F')
C = 0.25*(C[:-1, :-1] + C[1:, :-1] + C[:-1, 1:] + C[1:, 1:])
elif imageType == 'Fx':
C = I[:].reshape(self.vnFx, order='F')
C = 0.5*(C[:-1, :] + C[1:, :] )
elif imageType == 'Fy':
C = I[:].reshape(self.vnFy, order='F')
C = 0.5*(C[:, :-1] + C[:, 1:] )
elif imageType == 'Ex':
C = I[:].reshape(self.vnEx, order='F')
C = 0.5*(C[:,:-1] + C[:,1:] )
elif imageType == 'Ey':
C = I[:].reshape(self.vnEy, order='F')
C = 0.5*(C[:-1,:] + C[1:,:] )
return self._plotImage2D(v, vType=vType, grid=grid, view=view,
ax=ax, clim=clim, showIt=showIt,
pcolorOpts=pcolorOpts, streamOpts=streamOpts,
gridOpts=gridOpts)
elif self.dim == 3:
# get copy of image and average to cell-centers is necessary
if vType == 'CC':
vc = v.reshape(self.vnC, order='F')
elif vType == 'N':
vc = (self.aveN2CC*v).reshape(self.vnC, order='F')
elif vType in ['Fx', 'Fy', 'Fz', 'Ex', 'Ey', 'Ez']:
aveOp = 'ave' + vType[0] + '2CCV'
v = getattr(self,aveOp)*v # average to cell centers
ind_xyz = {'x':0,'y':1,'z':2}[vType[1]]
vc = self.r(v.reshape((self.nC,-1),order='F'), 'CC','CC','M')[ind_xyz]
# determine number oE slices in x and y dimension
nX = np.ceil(np.sqrt(self.nCz))
nY = np.ceil(self.nCz/nX)
# allocate space for montage
nCx = self.nCx
nCy = self.nCy
C = np.zeros((nX*nCx,nY*nCy))
for iy in range(int(nY)):
for ix in range(int(nX)):
iz = ix + iy*nX
if iz < self.nCz:
C[ix*nCx:(ix+1)*nCx, iy*nCy:(iy+1)*nCy] = vc[:, :, iz]
else:
C[ix*nCx:(ix+1)*nCx, iy*nCy:(iy+1)*nCy] = np.nan
C = np.ma.masked_where(np.isnan(C), C)
xx = np.r_[0, np.cumsum(np.kron(np.ones((nX, 1)), self.hx).ravel())]
yy = np.r_[0, np.cumsum(np.kron(np.ones((nY, 1)), self.hy).ravel())]
# Plot the mesh
if clim is None:
clim = [C.min(),C.max()]
ph = ax.pcolormesh(self.vectorNx, self.vectorNy, C.T, vmin=clim[0], vmax=clim[1])
ph = ax.pcolormesh(xx, yy, C.T, vmin=clim[0], vmax=clim[1])
# Plot the lines
gx = np.arange(nX+1)*(self.vectorNx[-1]-self.x0[0])
gy = np.arange(nY+1)*(self.vectorNy[-1]-self.x0[1])
# Repeat and seperate with NaN
gxX = np.c_[gx, gx, gx+np.nan].ravel()
gxY = np.kron(np.ones((nX+1, 1)), np.array([0, sum(self.hy)*nY, np.nan])).ravel()
gyX = np.kron(np.ones((nY+1, 1)), np.array([0, sum(self.hx)*nX, np.nan])).ravel()
gyY = np.c_[gy, gy, gy+np.nan].ravel()
ax.plot(gxX, gxY, annotationColor+'-', linewidth=2)
ax.plot(gyX, gyY, annotationColor+'-', linewidth=2)
ax.axis('tight')
ax.set_xlabel("x")
ax.set_ylabel("y")
elif self.dim == 3:
if direction == 'z':
# get copy of image and average to cell-centres is necessary
if imageType == 'CC':
Ic = I[:].reshape(self.vnC, order='F')
elif imageType == 'N':
Ic = I[:].reshape(self.vnN, order='F')
Ic = .125*(Ic[:-1,:-1,:-1]+Ic[1:,:-1,:-1] + Ic[:-1,1:,:-1]+ Ic[1:,1:,:-1]+ Ic[:-1,:-1,1:]+Ic[1:,:-1,1:] + Ic[:-1,1:,1:]+ Ic[1:,1:,1:] )
elif imageType == 'Fx':
Ic = I[:].reshape(self.vnFx, order='F')
Ic = .5*(Ic[:-1,:,:]+Ic[1:,:,:])
elif imageType == 'Fy':
Ic = I[:].reshape(self.vnFy, order='F')
Ic = .5*(Ic[:,:-1,:]+Ic[:,1:,:])
elif imageType == 'Fz':
Ic = I[:].reshape(self.vnFz, order='F')
Ic = .5*(Ic[:,:,:-1]+Ic[:,:,1:])
elif imageType == 'Ex':
Ic = I[:].reshape(self.vnEx, order='F')
Ic = .25*(Ic[:,:-1,:-1]+Ic[:,1:,:-1]+Ic[:,:-1,1:]+Ic[:,1:,:1])
elif imageType == 'Ey':
Ic = I[:].reshape(self.vnEy, order='F')
Ic = .25*(Ic[:-1,:,:-1]+Ic[1:,:,:-1]+Ic[:-1,:,1:]+Ic[1:,:,:1])
elif imageType == 'Ez':
Ic = I[:].reshape(self.vnEz, order='F')
Ic = .25*(Ic[:-1,:-1,:]+Ic[1:,:-1,:]+Ic[:-1,1:,:]+Ic[1:,:1,:])
# determine number oE slices in x and y dimension
nX = np.ceil(np.sqrt(self.nCz))
nY = np.ceil(self.nCz/nX)
# allocate space for montage
nCx = self.nCx
nCy = self.nCy
C = np.zeros((nX*nCx,nY*nCy))
if numbering:
pad = np.sum(self.hx)*0.04
for iy in range(int(nY)):
for ix in range(int(nX)):
iz = ix + iy*nX
if iz < self.nCz:
C[ix*nCx:(ix+1)*nCx, iy*nCy:(iy+1)*nCy] = Ic[:, :, iz]
else:
C[ix*nCx:(ix+1)*nCx, iy*nCy:(iy+1)*nCy] = np.nan
ax.text((ix+1)*(self.vectorNx[-1]-self.x0[0])-pad,(iy)*(self.vectorNy[-1]-self.x0[1])+pad,
'#%i'%iz,color=annotationColor,verticalalignment='bottom',horizontalalignment='right',size='x-large')
C = np.ma.masked_where(np.isnan(C), C)
xx = np.r_[0, np.cumsum(np.kron(np.ones((nX, 1)), self.hx).ravel())]
yy = np.r_[0, np.cumsum(np.kron(np.ones((nY, 1)), self.hy).ravel())]
# Plot the mesh
if clim is None:
clim = [C.min(),C.max()]
ph = ax.pcolormesh(xx, yy, C.T, vmin=clim[0], vmax=clim[1])
# Plot the lines
gx = np.arange(nX+1)*(self.vectorNx[-1]-self.x0[0])
gy = np.arange(nY+1)*(self.vectorNy[-1]-self.x0[1])
# Repeat and seperate with NaN
gxX = np.c_[gx, gx, gx+np.nan].ravel()
gxY = np.kron(np.ones((nX+1, 1)), np.array([0, sum(self.hy)*nY, np.nan])).ravel()
gyX = np.kron(np.ones((nY+1, 1)), np.array([0, sum(self.hx)*nX, np.nan])).ravel()
gyY = np.c_[gy, gy, gy+np.nan].ravel()
ax.plot(gxX, gxY, annotationColor+'-', linewidth=2)
ax.plot(gyX, gyY, annotationColor+'-', linewidth=2)
ax.axis('tight')
if numbering:
pad = np.sum(self.hx)*0.04
for iy in range(int(nY)):
for ix in range(int(nX)):
iz = ix + iy*nX
if iz < self.nCz:
ax.text((ix+1)*(self.vectorNx[-1]-self.x0[0])-pad,(iy)*(self.vectorNy[-1]-self.x0[1])+pad,
'#%i'%iz,color=annotationColor,verticalalignment='bottom',horizontalalignment='right',size='x-large')
ax.set_title(imageType)
ax.set_title(vType)
if showIt: plt.show()
return ph
@@ -270,7 +194,7 @@ class TensorView(object):
# Some user error checking
assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts)
assert self.dim == 3, 'Must be a 3D mesh.'
assert self.dim == 3, 'Must be a 3D mesh. Use plotImage.'
assert view in viewOpts, "view must be in ['%s']" % "','".join(viewOpts)
assert normal in normalOpts, "normal must be in ['%s']" % "','".join(normalOpts)
assert type(grid) is bool, 'grid must be a boolean'
@@ -279,14 +203,6 @@ class TensorView(object):
if ind is None: ind = int(szSliceDim/2)
assert type(ind) in [int, long], 'ind must be an integer'
if ax is None:
fig = plt.figure(1)
fig.clf()
ax = plt.subplot(111)
else:
assert isinstance(ax, matplotlib.axes.Axes), "ax must be an matplotlib.axes.Axes"
fig = ax.figure
# The slicing and plotting code!!
def getIndSlice(v):
@@ -299,76 +215,139 @@ class TensorView(object):
if vType == 'CC':
return getIndSlice(self.r(v,'CC','CC','M'))
elif vType == 'CCv':
v = self.r(v.reshape((self.nC,3),order='F'),'CC','CC','M')
assert view == 'vec', 'Other types for CCv not yet supported'
assert view == 'vec', 'Other types for CCv not supported'
else:
# Now just deal with 'F' and 'E'
aveOp = 'ave' + vType + ('2CCV' if view == 'vec' else '2CC')
v = getattr(self,aveOp)*v # average to cell centers (might be a vector)
v = self.r(v.reshape((self.nC,-1),order='F'),'CC','CC','M')
v = self.r(v.reshape((self.nC,-1),order='F'),'CC','CC','M')
if view == 'vec':
outSlice = []
if 'X' not in normal: outSlice.append(getIndSlice(v[0]))
if 'Y' not in normal: outSlice.append(getIndSlice(v[1]))
if 'Z' not in normal: outSlice.append(getIndSlice(v[2]))
return outSlice
return np.r_[mkvc(outSlice[0]), mkvc(outSlice[1])]
else:
return getIndSlice(self.r(v,'CC','CC','M'))
h2d = []
if 'X' not in normal: h2d.append(self.hx)
if 'Y' not in normal: h2d.append(self.hy)
if 'Z' not in normal: h2d.append(self.hz)
tM = self.__class__(h2d) #: Temp Mesh
x2d = []
if 'X' not in normal:
h2d.append(self.hx)
x2d.append(self.x0[0])
if 'Y' not in normal:
h2d.append(self.hy)
x2d.append(self.x0[1])
if 'Z' not in normal:
h2d.append(self.hz)
x2d.append(self.x0[2])
tM = self.__class__(h2d, x2d) #: Temp Mesh
v2d = doSlice(v)
if ax is None:
fig = plt.figure()
ax = plt.subplot(111)
else:
assert isinstance(ax, matplotlib.axes.Axes), "ax must be an matplotlib.axes.Axes"
fig = ax.figure
tM._plotImage2D(v2d, vType=('CCv' if view == 'vec' else 'CC'), grid=grid, view=view,
ax=ax, clim=clim, showIt=showIt,
pcolorOpts=pcolorOpts, streamOpts=streamOpts,
gridOpts=gridOpts)
ax.set_xlabel('y' if normal == 'X' else 'x')
ax.set_ylabel('y' if normal == 'Z' else 'z')
ax.set_title('Slice %d' % ind)
def _plotImage2D(self, v, vType='CC', grid=False, view='real',
ax=None, clim=None, showIt=False,
pcolorOpts={},
streamOpts={'color':'k'},
gridOpts={'color':'k'}
):
vTypeOptsCC = ['N','CC','Fx','Fy','Ex','Ey']
vTypeOptsV = ['CCv','F','E']
vTypeOpts = vTypeOptsCC + vTypeOptsV
if view == 'vec':
assert vType in vTypeOptsV, "vType must be in ['%s'] when view='vec'" % "','".join(vTypeOptsV)
assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts)
viewOpts = ['real','imag','abs','vec']
assert view in viewOpts, "view must be in ['%s']" % "','".join(viewOpts)
if ax is None:
fig = plt.figure()
ax = plt.subplot(111)
else:
assert isinstance(ax, matplotlib.axes.Axes), "ax must be an matplotlib.axes.Axes"
fig = ax.figure
# Reshape to a cell centered variable
if vType == 'CC':
pass
elif vType == 'CCv':
assert view == 'vec', 'Other types for CCv not supported'
elif vType in ['F', 'E', 'N']:
aveOp = 'ave' + vType + ('2CCV' if view == 'vec' else '2CC')
v = getattr(self,aveOp)*v # average to cell centers (might be a vector)
elif vType in ['Fx','Fy','Ex','Ey']:
aveOp = 'ave' + vType[0] + '2CCV'
v = getattr(self,aveOp)*v # average to cell centers (might be a vector)
xORy = {'x':0,'y':1}[vType[1]]
v = v.reshape((self.nC,-1), order='F')[:,xORy]
out = ()
if view in ['real','imag','abs']:
v = self.r(v, 'CC', 'CC', 'M')
v = getattr(np,view)(v) # e.g. np.real(v)
v = doSlice(v)
if clim is None:
clim = [v.min(),v.max()]
out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, v.T, vmin=clim[0], vmax=clim[1], **pcolorOpts),)
out += (ax.pcolormesh(self.vectorNx, self.vectorNy, v.T, vmin=clim[0], vmax=clim[1], **pcolorOpts),)
elif view in ['vec']:
U, V = doSlice(v)
U, V = self.r(v.reshape((self.nC,-1), order='F'), 'CC', 'CC', 'M')
if clim is None:
uv = np.r_[mkvc(U), mkvc(V)]
uv = np.sqrt(uv**2)
uv = np.sqrt(U**2 + V**2)
clim = [uv.min(),uv.max()]
# Matplotlib seems to not support irregular
# spaced vectors at the moment. So we will
# Interpolate down to a regular mesh at the
# smallest mesh size in this 2D slice.
nxi = int(tM.hx.sum()/tM.hx.min())
nyi = int(tM.hy.sum()/tM.hy.min())
tMi = self.__class__([np.ones(nxi)*tM.hx.sum()/nxi,
np.ones(nyi)*tM.hy.sum()/nyi])
P = tM.getInterpolationMat(tMi.gridCC,'CC',zerosOutside=True)
Ui = P*mkvc(U)
Vi = P*mkvc(V)
Ui = tMi.r(Ui, 'CC', 'CC', 'M')
Vi = tMi.r(Vi, 'CC', 'CC', 'M')
nxi = int(self.hx.sum()/self.hx.min())
nyi = int(self.hy.sum()/self.hy.min())
tMi = self.__class__([np.ones(nxi)*self.hx.sum()/nxi,
np.ones(nyi)*self.hy.sum()/nyi], self.x0)
P = self.getInterpolationMat(tMi.gridCC,'CC',zerosOutside=True)
Ui = tMi.r(P*mkvc(U), 'CC', 'CC', 'M')
Vi = tMi.r(P*mkvc(V), 'CC', 'CC', 'M')
# End Interpolation
out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, np.sqrt(U**2+V**2).T, vmin=clim[0], vmax=clim[1], **pcolorOpts),)
out += (ax.pcolormesh(self.vectorNx, self.vectorNy, np.sqrt(U**2+V**2).T, vmin=clim[0], vmax=clim[1], **pcolorOpts),)
out += (ax.streamplot(tMi.vectorCCx, tMi.vectorCCy, Ui.T, Vi.T, **streamOpts),)
if grid:
xXGrid = np.c_[tM.vectorNx,tM.vectorNx,np.nan*np.ones(tM.nNx)].flatten()
xYGrid = np.c_[tM.vectorNy[0]*np.ones(tM.nNx),tM.vectorNy[-1]*np.ones(tM.nNx),np.nan*np.ones(tM.nNx)].flatten()
yXGrid = np.c_[tM.vectorNx[0]*np.ones(tM.nNy),tM.vectorNx[-1]*np.ones(tM.nNy),np.nan*np.ones(tM.nNy)].flatten()
yYGrid = np.c_[tM.vectorNy,tM.vectorNy,np.nan*np.ones(tM.nNy)].flatten()
xXGrid = np.c_[self.vectorNx,self.vectorNx,np.nan*np.ones(self.nNx)].flatten()
xYGrid = np.c_[self.vectorNy[0]*np.ones(self.nNx),self.vectorNy[-1]*np.ones(self.nNx),np.nan*np.ones(self.nNx)].flatten()
yXGrid = np.c_[self.vectorNx[0]*np.ones(self.nNy),self.vectorNx[-1]*np.ones(self.nNy),np.nan*np.ones(self.nNy)].flatten()
yYGrid = np.c_[self.vectorNy,self.vectorNy,np.nan*np.ones(self.nNy)].flatten()
out += (ax.plot(np.r_[xXGrid,yXGrid],np.r_[xYGrid,yYGrid],**gridOpts)[0],)
ax.set_xlabel('y' if normal == 'X' else 'x')
ax.set_ylabel('y' if normal == 'Z' else 'z')
ax.set_title('Slice %d' % ind)
ax.set_xlim(*tM.vectorNx[[0,-1]])
ax.set_ylim(*tM.vectorNy[[0,-1]])
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_xlim(*self.vectorNx[[0,-1]])
ax.set_ylim(*self.vectorNy[[0,-1]])
if showIt: plt.show()
return out
def plotGrid(self, ax=None, nodes=False, faces=False, centers=False, edges=False, lines=True, showIt=False):
"""Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions.
@@ -401,15 +380,20 @@ class TensorView(object):
"""
axOpts = {'projection':'3d'} if self.dim == 3 else {}
if ax is None: ax = plt.subplot(111, **axOpts)
if ax is None:
fig = plt.figure()
ax = plt.subplot(111, **axOpts)
else:
assert isinstance(ax, matplotlib.axes.Axes), "ax must be an matplotlib.axes.Axes"
fig = ax.figure
if self.dim == 1:
if nodes:
ax.plot(xn, np.ones(self.nN), 'bs')
ax.plot(self.gridN, np.ones(self.nN), 'bs')
if centers:
ax.plot(xc, np.ones(self.nC), 'ro')
ax.plot(self.gridCC, np.ones(self.nC), 'ro')
if lines:
ax.plot(xn, np.ones(self.nN), 'b-')
ax.plot(self.gridN, np.ones(self.nN), 'b.-')
ax.set_xlabel('x1')
elif self.dim == 2:
if nodes:
@@ -636,11 +620,31 @@ class LomView(object):
if __name__ == '__main__':
from SimPEG import *
mT = Utils.meshTensors(((2,5),(4,2),(2,5)),((2,2),(6,2),(2,2)),((2,2),(6,2),(2,2)))
M = Mesh.TensorMesh(mT)
M = Mesh.TensorMesh(mT, x0=[10,20,14])
q = np.zeros(M.vnC)
q[[4,4],[4,4],[2,6]]=[-1,1]
q = Utils.mkvc(q)
A = M.faceDiv*M.cellGrad
b = Solver(A).solve(q)
M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, showIt=True, pcolorOpts={'alpha':0.8})
M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, pcolorOpts={'alpha':0.8})
M2 = Mesh.TensorMesh([10,20],x0=[10,5])
f = np.r_[np.sin(M2.gridFx[:,0]*2*np.pi), np.sin(M2.gridFy[:,1]*2*np.pi)]
M2.plotImage(f, 'F', view='vec', grid=True, pcolorOpts={'alpha':0.8})
M2.plotImage(f, 'Fx')
f = np.r_[np.sin(M2.gridEx[:,0]*2*np.pi), np.sin(M2.gridEy[:,1]*2*np.pi)]
M2.plotImage(f, 'E', view='vec', grid=True, pcolorOpts={'alpha':0.8})
c = np.r_[np.sin(M2.gridCC[:,0]*2*np.pi)]
M2.plotImage(c, 'CC', view='real')
from SimPEG import Mesh, np
M = Mesh.TensorMesh([20,20,20])
v = np.sin(M.gridCC[:,0]*2*np.pi)*np.sin(M.gridCC[:,1]*2*np.pi)*np.sin(M.gridCC[:,2]*2*np.pi)
M.plotImage(v, annotationColor='k')
Mesh.TensorMesh([10]).plotGrid()
plt.show()
SimPEG/ObjFunction.py
+23 -19
View File
@@ -1,7 +1,7 @@
import Utils, Parameters, numpy as np, scipy.sparse as sp
import Utils, Parameters, Survey, Problem, numpy as np, scipy.sparse as sp
class BaseObjFunction(object):
"""BaseObjFunction(data, reg, **kwargs)"""
"""BaseObjFunction(forward, reg, **kwargs)"""
__metaclass__ = Utils.SimPEGMetaClass
@@ -10,6 +10,9 @@ class BaseObjFunction(object):
debug = False #: Print debugging information
counter = None #: Set this to a SimPEG.Utils.Counter() if you want to count things
surveyPair = Survey.BaseSurvey
problemPair = Problem.BaseProblem
name = 'Base Objective Function' #: Name of the objective function
u_current = None #: The most current evaluated field
@@ -31,18 +34,19 @@ class BaseObjFunction(object):
def objFunc(self): return self
@property
def opt(self): return getattr(self.parent,'opt',None)
@property
def prob(self): return self.data.prob
@property
def mesh(self): return self.data.prob.mesh
@property
def model(self): return self.data.prob.model
def __init__(self, data, reg, **kwargs):
def __init__(self, forward, reg, **kwargs):
Utils.setKwargs(self, **kwargs)
self.data = data
assert forward.ispaired, 'The forward problem and survey must be paired.'
if isinstance(forward, self.surveyPair):
self.survey = forward
self.prob = forward.prob
elif isinstance(forward, self.problemPair):
self.prob = forward
self.survey = forward.survey
self.reg = reg
self.reg.parent = self
@@ -73,13 +77,13 @@ class BaseObjFunction(object):
self.u_current = None
self.m_current = m
u = self.data.prob.fields(m)
u = self.prob.fields(m)
self.u_current = u
phi_d = self.dataObj(m, u=u)
phi_m = self.reg.modelObj(m)
self.dpred = self.data.dpred(m, u=u) # This is a cheap matrix vector calculation.
self.dpred = self.survey.dpred(m, u=u) # This is a cheap matrix vector calculation.
self.phi_d, self.phi_d_last = phi_d, self.phi_d
self.phi_m, self.phi_m_last = phi_m, self.phi_m
@@ -124,7 +128,7 @@ class BaseObjFunction(object):
u is the field of interest; d_obs is the observed data; and W is the weighting matrix.
"""
# TODO: ensure that this is a data is vector and Wd is a matrix.
R = self.data.residualWeighted(m, u=u)
R = self.survey.residualWeighted(m, u=u)
return 0.5*np.vdot(R, R)
@Utils.timeIt
@@ -160,11 +164,11 @@ class BaseObjFunction(object):
\\frac{\partial \mu_\\text{data}}{\partial \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ R}
"""
if u is None: u = self.data.prob.fields(m)
if u is None: u = self.prob.fields(m)
R = self.data.residualWeighted(m, u=u)
R = self.survey.residualWeighted(m, u=u)
dmisfit = self.data.prob.Jtvec(m, self.data.Wd * R, u=u)
dmisfit = self.prob.Jtvec(m, self.survey.Wd * R, u=u)
return dmisfit
@@ -204,12 +208,12 @@ class BaseObjFunction(object):
\\frac{\partial^2 \mu_\\text{data}}{\partial^2 \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ W J}
"""
if u is None: u = self.data.prob.fields(m)
if u is None: u = self.prob.fields(m)
R = self.data.residualWeighted(m, u=u)
R = self.survey.residualWeighted(m, u=u)
# TODO: abstract to different norms a little cleaner.
# \/ it goes here. in l2 it is the identity.
dmisfit = self.data.prob.Jtvec_approx(m, self.data.Wd * self.data.Wd * self.data.prob.Jvec_approx(m, v, u=u), u=u)
dmisfit = self.prob.Jtvec_approx(m, self.survey.Wd * self.survey.Wd * self.prob.Jvec_approx(m, v, u=u), u=u)
return dmisfit
SimPEG/Parameters.py
+3 -3
View File
@@ -41,7 +41,7 @@ class Parameter(object):
@property
def reg(self): return self.parent.reg
@property
def data(self): return self.parent.data
def survey(self): return self.parent.survey
@property
def prob(self): return self.parent.prob
@property
@@ -131,13 +131,13 @@ class BetaEstimate(Parameter):
:return: beta0
"""
objFunc = self.parent
data = objFunc.data
survey = objFunc.survey
m = objFunc.m_current
u = objFunc.u_current
if u is None:
u = data.prob.fields(m)
u = survey.prob.fields(m)
x0 = np.random.rand(*m.shape)
t = x0.dot(objFunc.dataObj2Deriv(m,x0,u=u))
SimPEG/Problem.py
+29 -26
View File
@@ -1,4 +1,4 @@
import Utils, Data, numpy as np, scipy.sparse as sp
import Utils, Survey, numpy as np, scipy.sparse as sp
import Model
class BaseProblem(object):
@@ -38,40 +38,42 @@ class BaseProblem(object):
counter = None #: A SimPEG.Utils.Counter object
dataPair = Data.BaseData
surveyPair = Survey.BaseSurvey
modelPair = Model.BaseModel
def __init__(self, mesh, model, **kwargs):
def __init__(self, model, **kwargs):
Utils.setKwargs(self, **kwargs)
self.mesh = mesh
assert (isinstance(model, self.modelPair) or
isinstance(model, Model.ComboModel) and isinstance(model.models[0], self.modelPair)
), "Model object must be an instance of a %s class."%(self.modelPair.__name__)
self.model = model
@property
def data(self):
def mesh(self): return self.model.mesh
@property
def survey(self):
"""
The data object for this problem.
The survey object for this problem.
"""
return getattr(self, '_data', None)
return getattr(self, '_survey', None)
def pair(self, d):
"""Bind a data to this problem instance using pointers."""
assert isinstance(d, self.dataPair), "Data object must be an instance of a %s class."%(self.dataPair.__name__)
"""Bind a survey to this problem instance using pointers."""
assert isinstance(d, self.surveyPair), "Data object must be an instance of a %s class."%(self.surveyPair.__name__)
if d.ispaired:
raise Exception("The data object is already paired to a problem. Use data.unpair()")
self._data = d
raise Exception("The survey object is already paired to a problem. Use survey.unpair()")
self._survey = d
d._prob = self
def unpair(self):
"""Unbind a data from this problem instance."""
"""Unbind a survey from this problem instance."""
if not self.ispaired: return
self.data._prob = None
self._data = None
self.survey._prob = None
self._survey = None
@property
def ispaired(self): return self.data is not None
def ispaired(self): return self.survey is not None
@Utils.timeIt
def Jvec(self, m, v, u=None):
@@ -156,25 +158,26 @@ class BaseProblem(object):
"""
pass
def createSyntheticData(self, m, std=0.05, u=None, **geometry_kwargs):
#TODO: Rename and refactor to createSyntheticData
def createSyntheticSurvey(self, m, std=0.05, u=None, **geometry_kwargs):
"""
Create synthetic data given a model, and a standard deviation.
Create synthetic survey given a model, and a standard deviation.
:param numpy.array m: geophysical model
:param numpy.array std: standard deviation
:rtype: numpy.array, numpy.array
:return: dobs, Wd
:rtype: SurveyObject
:return: survey
Returns the observed data with random Gaussian noise
and Wd which is the same size as data, and can be used to weight the inversion.
"""
data = self.dataPair(mtrue=m, **geometry_kwargs)
data.pair(self)
data.dtrue = data.dpred(m, u=u)
noise = std*abs(data.dtrue)*np.random.randn(*data.dtrue.shape)
data.dobs = data.dtrue+noise
data.std = data.dobs*0 + std
return data
survey = self.surveyPair(mtrue=m, **geometry_kwargs)
survey.pair(self)
survey.dtrue = survey.dpred(m, u=u)
noise = std*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
survey.dobs = survey.dtrue+noise
survey.std = survey.dobs*0 + std
return survey
SimPEG/Regularization.py
+1 -1
View File
@@ -46,7 +46,7 @@ class BaseRegularization(object):
@property
def prob(self): return self.parent.prob
@property
def data(self): return self.parent.data
def survey(self): return self.parent.survey
@property
def mesh(self): return self.model.mesh
SimPEG/Data.py → SimPEG/Survey.py
+93 -24
View File
@@ -1,8 +1,8 @@
import Utils, numpy as np
class BaseData(object):
"""Data holds the observed data, and the standard deviations."""
class BaseSurvey(object):
"""Survey holds the observed data, and the standard deviations."""
__metaclass__ = Utils.SimPEGMetaClass
@@ -19,25 +19,32 @@ class BaseData(object):
@property
def prob(self):
"""
The geophysical problem that explains this data, use::
The geophysical problem that explains this survey, use::
data.pair(prob)
survey.pair(prob)
"""
return getattr(self, '_prob', None)
@property
def mesh(self):
"""Mesh of the paired problem."""
if self.ispaired:
return self.prob.mesh
raise Exception('Pair survey to a problem to access the problems mesh.')
def pair(self, p):
"""Bind a problem to this data instance using pointers"""
assert hasattr(p, 'dataPair'), "Problem must have an attribute 'dataPair'."
assert isinstance(self, p.dataPair), "Problem requires data object must be an instance of a %s class."%(p.dataPair.__name__)
"""Bind a problem to this survey instance using pointers"""
assert hasattr(p, 'surveyPair'), "Problem must have an attribute 'surveyPair'."
assert isinstance(self, p.surveyPair), "Problem requires survey object must be an instance of a %s class."%(p.surveyPair.__name__)
if p.ispaired:
raise Exception("The problem object is already paired to a data. Use prob.unpair()")
raise Exception("The problem object is already paired to a survey. Use prob.unpair()")
self._prob = p
p._data = self
p._survey = self
def unpair(self):
"""Unbind a problem from this data instance"""
"""Unbind a problem from this survey instance"""
if not self.ispaired: return
self.prob._data = None
self.prob._survey = None
self._prob = None
@property
@@ -71,7 +78,7 @@ class BaseData(object):
d_\\text{pred} = \mathbf{P} u(m)
"""
return u
raise NotImplemented('projectFields is not yet implemented.')
@Utils.count
def projectFieldsDeriv(self, u):
@@ -108,7 +115,7 @@ class BaseData(object):
"""
Data weighting matrix. This is a covariance matrix used in::
def data.residualWeighted(m,u=None):
def residualWeighted(m,u=None):
return self.Wd*self.residual(m, u=u)
By default, this is based on the norm of the data plus a noise floor.
@@ -139,20 +146,82 @@ class BaseData(object):
"""
return Utils.mkvc(self.Wd*self.residual(m, u=u))
@property
def RHS(self):
"""
Source matrix.
"""
return getattr(self, '_RHS', None)
@RHS.setter
def RHS(self, value):
self._RHS = value
@property
def isSynthetic(self):
"Check if the data is synthetic."
return (self.mtrue is not None)
return self.mtrue is not None
#TODO: Move this to the data class?
# @property
# def phi_d_target(self):
# """
# target for phi_d
# By default this is the number of data.
# Note that we do not set the target if it is None, but we return the default value.
# """
# if getattr(self, '_phi_d_target', None) is None:
# return self.data.dobs.size #
# return self._phi_d_target
# @phi_d_target.setter
# def phi_d_target(self, value):
# self._phi_d_target = value
class BaseRxList(object):
"""SimPEG Receiver List Object"""
locs = None #: Locations (nRx x 3)
knownRxTypes = None #: Set this to a list of strings to ensure that txType is known
def __init__(self, locs, rxType, **kwargs):
self.locs = locs
self.rxType = rxType
Utils.setKwargs(self, **kwargs)
@property
def rxType(self):
"""Receiver Type"""
return getattr(self, '_rxType', None)
@rxType.setter
def rxType(self, value):
known = self.knownRxTypes
if known is not None:
assert value in known, "rxType must be in ['%s']" % ("', '".join(known))
self._rxType = value
class BaseTx(object):
"""SimPEG Transmitter Object"""
loc = None #: Location [x,y,z]
rxList = None #: SimPEG Receiver List
rxListPair = BaseRxList
knownTxTypes = None #: Set this to a list of strings to ensure that txType is known
def __init__(self, loc, txType, rxList, **kwargs):
assert isinstance(rxList, self.rxListPair), 'rxList must be a %s'%self.rxListPair.__name__
self.loc = loc
self.txType = txType
self.rxList = rxList
Utils.setKwargs(self, **kwargs)
@property
def txType(self):
"""Transmitter Type"""
return getattr(self, '_txType', None)
@txType.setter
def txType(self, value):
known = self.knownTxTypes
if known is not None:
assert value in known, "txType must be in ['%s']" % ("', '".join(known))
self._txType = value
if __name__ == '__main__':
d = BaseData()
SimPEG/Utils/__init__.py
+2 -2
View File
@@ -159,11 +159,11 @@ def requires(var):
.. note::
To use data.%s(), SimPEG requires that a problem be bound to the data.
To use survey.%s(), SimPEG requires that a problem be bound to the survey.
If a problem has not been bound, an Exception will be raised.
To bind a problem to the Data object::
data.pair(myProblem)
survey.pair(myProblem)
""" % f.__name__
else:
SimPEG/Utils/matutils.py
+3
View File
@@ -30,6 +30,9 @@ def mkvc(x, numDims=1):
if type(x) == np.matrix:
x = np.array(x)
if hasattr(x, 'tovec'):
x = x.tovec()
assert type(x) == np.ndarray, "Vector must be a numpy array"
if numDims == 1:
SimPEG/__init__.py
+1 -1
View File
@@ -5,7 +5,7 @@ from Solver import Solver
import Mesh
import Model
import Problem
import Data
import Survey
import Regularization
import ObjFunction
import Optimization
Tutorials/Linear.py
+18 -12
View File
@@ -1,12 +1,16 @@
from SimPEG import *
import matplotlib.pyplot as plt
class LinearSurvey(Survey.BaseSurvey):
def projectFields(self, u):
return u
class LinearProblem(Problem.BaseProblem):
"""docstring for LinearProblem"""
def __init__(self, mesh, model, G, **kwargs):
Problem.BaseProblem.__init__(self, mesh, model, **kwargs)
surveyPair = LinearSurvey
def __init__(self, model, G, **kwargs):
Problem.BaseProblem.__init__(self, model, **kwargs)
self.G = G
def fields(self, m, u=None):
@@ -20,8 +24,7 @@ class LinearProblem(Problem.BaseProblem):
def example(N):
h = np.ones(N)/N
M = Mesh.TensorMesh([h])
M = Mesh.TensorMesh([N])
nk = 20
jk = np.linspace(1.,20.,nk)
@@ -43,22 +46,25 @@ def example(N):
model = Model.BaseModel(M)
prob = LinearProblem(M, model, G)
data = prob.createSyntheticData(mtrue, std=0.01)
prob = LinearProblem(model, G)
survey = prob.createSyntheticSurvey(mtrue, std=0.01)
return prob, data, model
return prob, survey, model
if __name__ == '__main__':
prob, data, model = example(100)
import matplotlib.pyplot as plt
prob, survey, model = example(100)
M = prob.mesh
reg = Regularization.Tikhonov(model)
objFunc = ObjFunction.BaseObjFunction(data, reg)
beta = Parameters.BetaSchedule()
objFunc = ObjFunction.BaseObjFunction(survey, reg, beta=beta)
opt = Optimization.InexactGaussNewton(maxIter=20)
inv = Inversion.BaseInversion(objFunc, opt)
m0 = np.zeros_like(data.mtrue)
m0 = np.zeros_like(survey.mtrue)
mrec = inv.run(m0)
@@ -68,7 +74,7 @@ if __name__ == '__main__':
plt.figure(2)
plt.plot(M.vectorCCx, data.mtrue, 'b-')
plt.plot(M.vectorCCx, survey.mtrue, 'b-')
plt.plot(M.vectorCCx, mrec, 'r-')
plt.show()
docs/simpeg-framework.png
BIN
View File
Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 16 KiB

After

Width:  |  Height:  |  Size: 20 KiB