diff --git a/SimPEG/Inversion.py b/SimPEG/Inversion.py index 297270ef..8dbee17d 100644 --- a/SimPEG/Inversion.py +++ b/SimPEG/Inversion.py @@ -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) diff --git a/SimPEG/Mesh/View.py b/SimPEG/Mesh/View.py index eaa9644c..5b2c1397 100644 --- a/SimPEG/Mesh/View.py +++ b/SimPEG/Mesh/View.py @@ -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() diff --git a/SimPEG/ObjFunction.py b/SimPEG/ObjFunction.py index 2b1912b8..18d9cfc3 100644 --- a/SimPEG/ObjFunction.py +++ b/SimPEG/ObjFunction.py @@ -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 diff --git a/SimPEG/Parameters.py b/SimPEG/Parameters.py index 5a6760e6..417051df 100644 --- a/SimPEG/Parameters.py +++ b/SimPEG/Parameters.py @@ -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)) diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index a509f6cd..c360f02c 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -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 diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 8dc32b79..8f7b9ede 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -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 diff --git a/SimPEG/Data.py b/SimPEG/Survey.py similarity index 53% rename from SimPEG/Data.py rename to SimPEG/Survey.py index 12696406..e754476e 100644 --- a/SimPEG/Data.py +++ b/SimPEG/Survey.py @@ -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() diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index a44cfd74..209d1cff 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -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: diff --git a/SimPEG/Utils/matutils.py b/SimPEG/Utils/matutils.py index f1ea7ec7..ac6c391d 100644 --- a/SimPEG/Utils/matutils.py +++ b/SimPEG/Utils/matutils.py @@ -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: diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 17f64d54..874f4ffe 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -5,7 +5,7 @@ from Solver import Solver import Mesh import Model import Problem -import Data +import Survey import Regularization import ObjFunction import Optimization diff --git a/Tutorials/Linear.py b/Tutorials/Linear.py index 7d3aa805..1bfa7339 100644 --- a/Tutorials/Linear.py +++ b/Tutorials/Linear.py @@ -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() diff --git a/docs/simpeg-framework.png b/docs/simpeg-framework.png index ad8e5c1c..f90f8260 100644 Binary files a/docs/simpeg-framework.png and b/docs/simpeg-framework.png differ