diff --git a/SimPEG/DCIP/DCIPUtils.py b/SimPEG/DCIP/DCIPUtils.py index 49f9ba2b..8e1ca5e9 100644 --- a/SimPEG/DCIP/DCIPUtils.py +++ b/SimPEG/DCIP/DCIPUtils.py @@ -169,7 +169,8 @@ def readUBC_DC2DModel(fileName): return model -def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt', clim=None): + +def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt', clim=None, cblabel=True, axlabel = True, colorbar = True, contour = None): """ Read list of 2D tx-rx location and plot a speudo-section of apparent resistivity. @@ -187,9 +188,6 @@ def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt from scipy.interpolate import griddata import pylab as plt - # Set depth to 0 for now - z0 = 0. - # Pre-allocate midx = [] midz = [] @@ -254,38 +252,54 @@ def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt midx = np.hstack([midx, ( Cmid + Pmid )/2 ]) midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + (Tx[0][2] + Tx[1][2])/2 ]) - ax = axs - # Grid points grid_x, grid_z = np.mgrid[np.min(midx):np.max(midx), np.min(midz):np.max(midz)] grid_rho = griddata(np.c_[midx,midz], rho.T, (grid_x, grid_z), method='linear') + # Scale the color scheme if clim == None: vmin, vmax = rho.min(), rho.max() else: vmin, vmax = clim[0], clim[1] + # Plot data grid_rho = np.ma.masked_where(np.isnan(grid_rho), grid_rho) - ph = plt.pcolormesh(grid_x[:,0],grid_z[0,:],grid_rho.T, clim=(vmin, vmax)) - cbar = plt.colorbar(format="$10^{%.1f}$",fraction=0.04,orientation="horizontal") - cmin,cmax = cbar.get_clim() - ticks = np.linspace(cmin,cmax,3) - cbar.set_ticks(ticks) - cbar.ax.tick_params(labelsize=10) - if unitType == 'appConductivity': - cbar.set_label("App.Cond",size=12) - elif unitType == 'appResistivity': - cbar.set_label("App.Res.",size=12) - elif unitType == 'volt': - cbar.set_label("Potential (V)",size=12) - # Plot apparent resistivity - ax.scatter(midx,midz,s=10,c=rho.T, vmin =vmin, vmax = vmax, clim=(vmin, vmax)) + ph = plt.pcolormesh(grid_x[:,0],grid_z[0,:],grid_rho.T, vmin = vmin, vmax = vmax) + plt.gca().tick_params(axis='both', which='major', labelsize=8) - #ax.set_xticklabels([]) - #ax.set_yticklabels([]) + if contour is not None: + plt.contour(grid_x,grid_z,grid_rho,levels = contour,colors = 'r', vmin = vmin, vmax = vmax) + + # Add scatter points + axs.scatter(midx,midz,s=10,c=rho.T, vmin = vmin, vmax = vmax) + + if colorbar: + + if unitType == 'volt': + cbar = plt.colorbar(ph, ax = axs, format="%4.1f",fraction=0.04,orientation="horizontal") + + else: + cbar = plt.colorbar(ph, ax = axs, format="$10^{%.1f}$",fraction=0.04,orientation="horizontal") + + cmin,cmax = cbar.get_clim() + ticks = np.linspace(cmin,cmax,3) + cbar.set_ticks(ticks) + cbar.ax.tick_params(labelsize=10) + + if unitType == 'appConductivity': + cbar.set_label("App.Cond",size=12) + elif unitType == 'appResistivity': + cbar.set_label("App.Res.",size=12) + elif unitType == 'volt': + cbar.set_label("Potential (V)",size=12) + + + if not axlabel: + axs.set_xticklabels([]) + axs.set_yticklabels([]) plt.gca().set_aspect('equal', adjustable='box') @@ -440,7 +454,8 @@ def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx): survey = DC.SurveyDC(SrcList) return survey, Tx, Rx -def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType): + +def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0): """ Write UBC GIF DCIP 2D or 3D observation file @@ -460,6 +475,12 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType): fid = open(fileName,'w') fid.write('! ' + surveyType + ' FORMAT\n') + if iptype!=0: + fid.write('IPTYPE=%i\n'%iptype) + + else: + fid.write('! ' + stype + ' FORMAT\n') + count = 0 for ii in range(DCsurvey.nSrc): @@ -491,18 +512,25 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType): if surveyType == 'SURFACE': - fid.writelines("%e " % ii for ii in mkvc(tx[0,:])) + fid.writelines("%f " % ii for ii in mkvc(tx[0,:])) M = M[:,0] N = N[:,0] if surveyType == 'GENERAL': + # Flip sign for z-elevation to depth + tx[2::2,:] = -tx[2::2,:] + fid.writelines("%e " % ii for ii in mkvc(tx[::2,:])) M = M[:,0::2] N = N[:,0::2] + # Flip sign for z-elevation to depth + M[:,1::2] = -M[:,1::2] + N[:,1::2] = -N[:,1::2] + fid.write('%i\n'% nD) - np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n') + np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%f',delimiter=' ',newline='\n') if dim=='3D': @@ -514,10 +542,11 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType): if surveyType == 'GENERAL': - fid.writelines("%e " % ii for ii in mkvc(tx)) + fid.writelines("%e " % ii for ii in mkvc(tx[0:3,:])) fid.write('%i\n'% nD) np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n') + fid.write('\n') count += nD @@ -620,43 +649,53 @@ def convertObs_DC3D_to_2D(DCsurvey, lineID, flag='local'): return DCsurvey2D -def readUBC_DC3Dobs(fileName): +def readUBC_DC3Dobs(fileName, rtype = 'DC'): """ - Read UBC GIF DCIP 3D observation file and generate survey + Read UBC GIF IP 3D observation file and generate survey :param string fileName:, path to the UBC GIF 3D obs file :rtype: Survey :return: DCIPsurvey """ + zflag = True # Flag for z value provided # Load file - obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!') + if rtype == 'IP': + obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='IPTYPE') + + elif rtype == 'DC': + obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!') + + else: + print "rtype must be 'DC'(default) | 'IP'" # Pre-allocate srcLists = [] Rx = [] d = [] wd = [] - zflag = True # Flag for z value provided + # Countdown for number of obs/tx count = 0 for ii in range(obsfile.shape[0]): + # Skip if blank line if not obsfile[ii]: continue - # First line is transmitter with number of receivers + # First line or end of a transmitter block, read transmitter info if count==0: - - temp = (np.fromstring(obsfile[ii], dtype=float,sep=' ').T) + # Read the line + temp = (np.fromstring(obsfile[ii], dtype=float, sep=' ').T) count = int(temp[-1]) # Check if z value is provided, if False -> nan if len(temp)==5: - tx = np.r_[temp[0:2],np.nan,temp[0:2],np.nan] - zflag = False + tx = np.r_[temp[0:2],np.nan,temp[2:4],np.nan] + + zflag = False # Pass on the flag to the receiver loc else: tx = temp[:-1] @@ -664,8 +703,16 @@ def readUBC_DC3Dobs(fileName): rx = [] continue - temp = np.fromstring(obsfile[ii], dtype=float,sep=' ') + temp = np.fromstring(obsfile[ii], dtype=float,sep=' ') # Get the string + # Filter out negative IP +# if temp[-2] < 0: +# count = count -1 +# print "Negative!" +# +# else: + + # If the Z-location is provided, otherwise put nan if zflag: rx.append(temp[:-2]) @@ -675,7 +722,7 @@ def readUBC_DC3Dobs(fileName): wd.append(temp[-1]) else: - rx.append(np.r_[temp[0:2],np.nan,temp[0:2],np.nan] ) + rx.append(np.r_[temp[0:2],np.nan,temp[2:4],np.nan] ) # Check if there is data with the location if len(temp)==6: d.append(temp[-2]) @@ -683,7 +730,7 @@ def readUBC_DC3Dobs(fileName): count = count -1 - # Reach the end of transmitter block + # Reach the end of transmitter block, append the src, rx and continue if count == 0: rx = np.asarray(rx) Rx = DC.RxDipole(rx[:,:3],rx[:,3:]) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index b694a2f1..93e51f62 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -146,10 +146,15 @@ class BetaSchedule(InversionDirective): class TargetMisfit(InversionDirective): + chifact = 1. + phi_d_star = None + @property def target(self): if getattr(self, '_target', None) is None: - self._target = self.survey.nD*0.5 + if self.phi_d_star is None: + self.phi_d_star = 0.5 * self.survey.nD + self._target = self.chifact * self.phi_d_star # the factor of 0.5 is because we do phid = 0.5*|| dpred - dobs||^2 return self._target @target.setter def target(self, val): @@ -258,6 +263,7 @@ class Update_IRLS(InversionDirective): phi_m_last = None phi_d_last = None + def initialize(self): # Scale the regularization for changes in norm @@ -275,7 +281,7 @@ class Update_IRLS(InversionDirective): self.phi_d_last = self.invProb.phi_d def endIter(self): - # Cool the threshold parameter + # Cool the threshold parameter if required if getattr(self, 'factor', None) is not None: eps = self.reg.eps / self.factor @@ -290,28 +296,44 @@ class Update_IRLS(InversionDirective): # Update the model used for the IRLS weights self.reg.curModel = self.invProb.curModel - # Temporarely set gamma to 1. + # Temporarely set gamma to 1. to get raw phi_m self.reg.gamma = 1. - # Compute change in model objective function and update scaling + # Compute new model objective function value phim_new = self.reg.eval(self.invProb.curModel) + # Update gamma to scale the regularization between IRLS iterations self.reg.gamma = self.phi_m_last / phim_new - self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d + # Set the weighting matrix to None so that it is recomputed next time + # it is called in the inversion + self.reg._W = None class Update_lin_PreCond(InversionDirective): + """ + Create a Jacobi preconditioner for the linear problem + """ + onlyOnStart=False + def initialize(self): + + if getattr(self.opt, 'approxHinv', None) is None: + # Update the pre-conditioner + diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() #* (self.reg.mapping * np.ones(self.reg.curModel.size))**2. + PC = Utils.sdiag((self.prob.mapping.deriv(None).T *diagA)**-1.) + self.opt.approxHinv = PC def endIter(self): # Cool the threshold parameter + if self.onlyOnStart==True: + return if getattr(self.opt, 'approxHinv', None) is not None: # Update the pre-conditioner - diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() * (self.reg.mapping * np.ones(self.reg.curModel.size))**2. - PC = Utils.sdiag(diagA**-1.) + diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() #* (self.reg.mapping * np.ones(self.reg.curModel.size))**2. + PC = Utils.sdiag((self.prob.mapping.deriv(None).T *diagA)**-1.) self.opt.approxHinv = PC - print 'Updated pre-cond' + class Update_Wj(InversionDirective): """ @@ -338,3 +360,19 @@ class Update_Wj(InversionDirective): JtJdiag = JtJdiag / max(JtJdiag) self.reg.wght = JtJdiag + +class Scale_Beta(InversionDirective): + """ + Instead of a linear cooling schedule, beta is allowed to change based + on the ratio between the target misfit and the current data misfit. The + update is done only if the misfit is outside some threshold bounds. + """ + tol = 0.05 + + def endIter(self): + + # Check if misfit is within the tolerance, otherwise adjust beta + val = self.invProb.phi_d / (self.survey.nD*0.5) + + if np.abs(1.-val) > self.tol: + self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d diff --git a/SimPEG/Examples/Inversion_IRLS.py b/SimPEG/Examples/Inversion_IRLS.py index c236860d..06ef4be4 100644 --- a/SimPEG/Examples/Inversion_IRLS.py +++ b/SimPEG/Examples/Inversion_IRLS.py @@ -86,12 +86,12 @@ def run(N=200, plotIt=True): #reg.recModel = mrec reg.wght = np.ones(mesh.nC) reg.mref = np.zeros(mesh.nC) - reg.eps_p = 2e-3 - reg.eps_q = 2e-3 + reg.eps_p = 5e-2 + reg.eps_q = 1e-2 reg.norms = [0., 0., 2., 2.] reg.wght = wr - opt = Optimization.ProjectedGNCG(maxIter=5 ,lower=-2.,upper=2., maxIterCG= 100, tolCG = 1e-3) + opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 20, tolCG = 1e-3) invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta*2.) beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1) #betaest = Directives.BetaEstimate_ByEig() diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index ecdf36ac..f0b6751e 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -330,7 +330,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView): raise NotImplementedError('wrapping in the averaging is not yet implemented') return self._aveF2CCV - def getInterpolationMatCartMesh(self, Mrect, locType='CC'): + def getInterpolationMatCartMesh(self, Mrect, locType='CC', locTypeTo=None): """ Takes a cartesian mesh and returns a projection to translate onto the cartesian grid. """ @@ -338,19 +338,22 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView): assert self.isSymmetric, "Currently we have not taken into account other projections for more complicated CylMeshes" + if locTypeTo is None: + locTypeTo = locType + if locType == 'F': # do this three times for each component - X = self.getInterpolationMatCartMesh(Mrect, locType='Fx') - Y = self.getInterpolationMatCartMesh(Mrect, locType='Fy') - Z = self.getInterpolationMatCartMesh(Mrect, locType='Fz') + X = self.getInterpolationMatCartMesh(Mrect, locType='Fx', locTypeTo=locTypeTo+'x') + Y = self.getInterpolationMatCartMesh(Mrect, locType='Fy', locTypeTo=locTypeTo+'y') + Z = self.getInterpolationMatCartMesh(Mrect, locType='Fz', locTypeTo=locTypeTo+'z') return sp.vstack((X,Y,Z)) if locType == 'E': - X = self.getInterpolationMatCartMesh(Mrect, locType='Ex') - Y = self.getInterpolationMatCartMesh(Mrect, locType='Ey') - Z = spzeros(Mrect.nEz, self.nE) + X = self.getInterpolationMatCartMesh(Mrect, locType='Ex', locTypeTo=locTypeTo+'x') + Y = self.getInterpolationMatCartMesh(Mrect, locType='Ey', locTypeTo=locTypeTo+'y') + Z = spzeros(getattr(Mrect, 'n' + locTypeTo + 'z'), self.nE) return sp.vstack((X,Y,Z)) - grid = getattr(Mrect, 'grid' + locType) + grid = getattr(Mrect, 'grid' + locTypeTo) # This is unit circle stuff, 0 to 2*pi, starting at x-axis, rotating counter clockwise in an x-y slice theta = - np.arctan2(grid[:,0] - self.cartesianOrigin[0], grid[:,1] - self.cartesianOrigin[1]) + np.pi/2 theta[theta < 0] += np.pi*2.0 @@ -366,7 +369,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView): 'Ex': Mrect.tangents[:Mrect.nEx,:], 'Ey': Mrect.tangents[Mrect.nEx:(Mrect.nEx+Mrect.nEy),:], 'Ez': Mrect.tangents[-Mrect.nEz:,:], - }[locType] + }[locTypeTo] if 'F' in locType: normals = np.c_[np.cos(theta), np.sin(theta), np.zeros(theta.size)] proj = ( normals * dotMe ).sum(axis=1) diff --git a/SimPEG/Mesh/MeshIO.py b/SimPEG/Mesh/MeshIO.py index 7501a66f..1c042237 100644 --- a/SimPEG/Mesh/MeshIO.py +++ b/SimPEG/Mesh/MeshIO.py @@ -21,10 +21,9 @@ class TensorMeshIO(object): if '*' in seg: st = seg sp = seg.split('*') - re = np.array(sp[0],dtype=int)*(' ' + sp[1]) + re = int(sp[0])*(' ' + sp[1]) line = line.replace(st,re.strip()) return np.array(line.split(),dtype=float) - # Read the file as line strings, remove lines with comment = ! msh = np.genfromtxt(fileName,delimiter='\n',dtype=np.str,comments='!') diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 195e6203..fc101a61 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -311,6 +311,9 @@ class BaseRegularization(object): tmp = indActive indActive = np.zeros(mesh.nC, dtype=bool) indActive[tmp] = True + if indActive is not None and mapping is None: + mapping = Maps.IdentityMap(nP=indActive.nonzero()[0].size) + self.regmesh = RegularizationMesh(mesh,indActive) self.mapping = mapping or self.mapPair(mesh) self.mapping._assertMatchesPair(self.mapPair) @@ -728,14 +731,14 @@ class Sparse(Simple): @property def W(self): """Full regularization matrix W""" - #if getattr(self, '_W', None) is None: - wlist = (self.Wsmall, self.Wsmooth) - #self._W = sp.vstack(wlist) - return sp.vstack(wlist) + if getattr(self, '_W', None) is None: + wlist = (self.Wsmall, self.Wsmooth) + self._W = sp.vstack(wlist) + return self._W def R(self, f_m , eps, exponent): - eta = (eps**(1-exponent/2.))**0.5 - r = eta / (f_m**2.+ eps**2.)**((1-exponent/2.)/2.) + eta = (eps**(1.-exponent/2.))**0.5 + r = eta / (f_m**2.+ eps**2.)**((1.-exponent/2.)/2.) return r diff --git a/SimPEG/Utils/io_utils.py b/SimPEG/Utils/io_utils.py new file mode 100644 index 00000000..323feb0c --- /dev/null +++ b/SimPEG/Utils/io_utils.py @@ -0,0 +1,137 @@ +from SimPEG import np, Mesh +import time as tm +import vtk, vtk.util.numpy_support as npsup +import re + +def read_GOCAD_ts(tsfile): + """ + + Read GOCAD triangulated surface (*.ts) file + INPUT: + tsfile: Triangulated surface + + OUTPUT: + vrts : Array of vertices in XYZ coordinates [n x 3] + trgl : Array of index for triangles [m x 3]. The order of the vertices + is important and describes the normal + n = cross( (P2 - P1 ) , (P3 - P1) ) + + Author: @fourndo + + + .. note:: + + Remove all attributes from the GoCAD surface before exporting it! + + """ + + + fid = open(tsfile,'r') + line = fid.readline() + + # Skip all the lines until the vertices + while re.match('TFACE',line)==None: + line = fid.readline() + + line = fid.readline() + vrtx = [] + + # Run down all the vertices and save in array + while re.match('VRTX',line): + l_input = re.split('[\s*]',line) + temp = np.array(l_input[2:5]) + vrtx.append(temp.astype(np.float)) + + # Read next line + line = fid.readline() + + vrtx = np.asarray(vrtx) + + # Skip lines to the triangles + while re.match('TRGL',line)==None: + line = fid.readline() + + # Run down the list of triangles + trgl = [] + + # Run down all the vertices and save in array + while re.match('TRGL',line): + l_input = re.split('[\s*]',line) + temp = np.array(l_input[1:4]) + trgl.append(temp.astype(np.int)) + + # Read next line + line = fid.readline() + + trgl = np.asarray(trgl) + + return vrtx, trgl + +def surface2inds(vrtx, trgl, mesh, boundaries=True, internal=True): + """" + Function to read gocad polystructure file and output indexes of mesh with in the structure. + + """ + # Adjust the index + trgl = trgl - 1 + + # Make vtk pts + ptsvtk = vtk.vtkPoints() + ptsvtk.SetData(npsup.numpy_to_vtk(vrtx,deep=1)) + + # Make the polygon connection + polys = vtk.vtkCellArray() + for face in trgl: + poly = vtk.vtkPolygon() + poly.GetPointIds().SetNumberOfIds(len(face)) + for nrv, vert in enumerate(face): + poly.GetPointIds().SetId(nrv,vert) + polys.InsertNextCell(poly) + + # Make the polydata, structure of connections and vrtx + polyData = vtk.vtkPolyData() + polyData.SetPoints(ptsvtk) + polyData.SetPolys(polys) + + # Make implicit func + ImpDistFunc = vtk.vtkImplicitPolyDataDistance() + ImpDistFunc.SetInput(polyData) + + # Convert the mesh + vtkMesh = vtk.vtkRectilinearGrid() + vtkMesh.SetDimensions(mesh.nNx,mesh.nNy,mesh.nNz) + vtkMesh.SetXCoordinates(npsup.numpy_to_vtk(mesh.vectorNx, deep=1)) + vtkMesh.SetYCoordinates(npsup.numpy_to_vtk(mesh.vectorNy, deep=1)) + vtkMesh.SetZCoordinates(npsup.numpy_to_vtk(mesh.vectorNz, deep=1)) + # Add indexes + vtkInd = npsup.numpy_to_vtk(np.arange(mesh.nC), deep=1) + vtkInd.SetName('Index') + vtkMesh.GetCellData().AddArray(vtkInd) + + extractImpDistRectGridFilt = vtk.vtkExtractGeometry() # Object constructor + extractImpDistRectGridFilt.SetImplicitFunction(ImpDistFunc) # + extractImpDistRectGridFilt.SetInputData(vtkMesh) + + if boundaries is True: + extractImpDistRectGridFilt.ExtractBoundaryCellsOn() + + else: + extractImpDistRectGridFilt.ExtractBoundaryCellsOff() + + if internal is True: + extractImpDistRectGridFilt.ExtractInsideOn() + + else: + extractImpDistRectGridFilt.ExtractInsideOff() + + print "Extracting indices from grid..." + # Executing the pipe + extractImpDistRectGridFilt.Update() + + # Get index inside + insideGrid = extractImpDistRectGridFilt.GetOutput() + insideGrid = npsup.vtk_to_numpy(insideGrid.GetCellData().GetArray('Index')) + + + # Return the indexes inside + return insideGrid diff --git a/setup.py b/setup.py index 3383c7f7..5308db47 100644 --- a/setup.py +++ b/setup.py @@ -5,16 +5,17 @@ SimPEG is a python package for simulation and gradient based parameter estimation in the context of geophysical applications. """ -import numpy as np - import os import sys import subprocess from distutils.core import setup +from distutils.command.build_ext import build_ext from setuptools import find_packages from distutils.extension import Extension + + CLASSIFIERS = [ 'Development Status :: 4 - Beta', 'Intended Audience :: Developers', @@ -51,11 +52,16 @@ if args.count("build_ext") > 0 and args.count("--inplace") == 0: try: from Cython.Build import cythonize from Cython.Distutils import build_ext - cythonKwargs = dict(cmdclass={'build_ext': build_ext}) USE_CYTHON = True except Exception, e: USE_CYTHON = False - cythonKwargs = dict() + +class NumpyBuild(build_ext): + def finalize_options(self): + build_ext.finalize_options(self) + __builtins__.__NUMPY_SETUP__ = False + import numpy + self.include_dirs.append(numpy.get_include()) ext = '.pyx' if USE_CYTHON else '.c' @@ -94,8 +100,8 @@ setup( classifiers=CLASSIFIERS, platforms = ["Windows", "Linux", "Solaris", "Mac OS-X", "Unix"], use_2to3 = False, - include_dirs=[np.get_include()], + cmdclass={'build_ext':NumpyBuild}, + setup_requires=['numpy'], ext_modules = extensions, scripts=scripts, - **cythonKwargs ) diff --git a/tests/base/test_regularization.py b/tests/base/test_regularization.py index 5aa21ad5..6dddfb74 100644 --- a/tests/base/test_regularization.py +++ b/tests/base/test_regularization.py @@ -65,10 +65,8 @@ class RegularizationTests(unittest.TestCase): elif mesh.dim == 3: indActive = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5) - mapping = Maps.IdentityMap(nP=indActive.nonzero()[0].size) - for indAct in [indActive, indActive.nonzero()[0]]: # test both bool and integers - reg = r(mesh, mapping=mapping, indActive=indAct) + reg = r(mesh, indActive=indAct) m = np.random.rand(mesh.nC)[indAct] reg.mref = np.ones_like(m)*np.mean(m) diff --git a/tests/mesh/test_cylMesh.py b/tests/mesh/test_cylMesh.py index 5a963398..c7cd44ae 100644 --- a/tests/mesh/test_cylMesh.py +++ b/tests/mesh/test_cylMesh.py @@ -146,6 +146,20 @@ class TestCyl2DMesh(unittest.TestCase): assert np.abs(Pr*(Pc2r*mc) - Pc*mc).max() < 1e-3 + def test_getInterpMatCartMesh_Cells2Nodes(self): + + Mr = Mesh.TensorMesh([100,100,2], x0='CC0') + Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0]) + + mc = np.arange(Mc.nC) + xr = np.linspace(0,0.4,50) + xc = np.linspace(0,0.4,50) + 0.2 + Pr = Mr.getInterpolationMat(np.c_[xr,np.ones(50)*-0.2,np.ones(50)*0.5],'N') + Pc = Mc.getInterpolationMat(np.c_[xc,np.zeros(50),np.ones(50)*0.5],'CC') + Pc2r = Mc.getInterpolationMatCartMesh(Mr, 'CC', locTypeTo='N') + + assert np.abs(Pr*(Pc2r*mc) - Pc*mc).max() < 1e-3 + def test_getInterpMatCartMesh_Faces(self): Mr = Mesh.TensorMesh([100,100,2], x0='CC0') @@ -177,6 +191,37 @@ class TestCyl2DMesh(unittest.TestCase): assert np.abs(mag[dist > 0.1].min() - 1) < TOL + def test_getInterpMatCartMesh_Faces2Edges(self): + + Mr = Mesh.TensorMesh([100,100,2], x0='CC0') + Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0]) + + Pf2e = Mc.getInterpolationMatCartMesh(Mr, 'F', locTypeTo='E') + mf = np.ones(Mc.nF) + + ecart = Pf2e * mf + + excc = Mr.aveEx2CC*Mr.r(ecart, 'E', 'Ex') + eycc = Mr.aveEy2CC*Mr.r(ecart, 'E', 'Ey') + ezcc = Mr.r(ecart, 'E', 'Ez') + + indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5]) + indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5]) + + TOL = 1e-2 + assert np.abs(float(excc[indX]) - 1) < TOL + assert np.abs(float(excc[indY]) - 0) < TOL + assert np.abs(float(eycc[indX]) - 0) < TOL + assert np.abs(float(eycc[indY]) - 1) < TOL + assert np.abs((ezcc - 1).sum()) < TOL + + mag = (excc**2 + eycc**2)**0.5 + dist = ((Mr.gridCC[:,0] + 0.2)**2 + (Mr.gridCC[:,1] + 0.2)**2)**0.5 + + assert np.abs(mag[dist > 0.1].max() - 1) < TOL + assert np.abs(mag[dist > 0.1].min() - 1) < TOL + + def test_getInterpMatCartMesh_Edges(self): Mr = Mesh.TensorMesh([100,100,2], x0='CC0') @@ -185,11 +230,42 @@ class TestCyl2DMesh(unittest.TestCase): Pe = Mc.getInterpolationMatCartMesh(Mr, 'E') me = np.ones(Mc.nE) - erect = Pe * me + ecart = Pe * me - excc = Mr.aveEx2CC*Mr.r(erect, 'E', 'Ex') - eycc = Mr.aveEy2CC*Mr.r(erect, 'E', 'Ey') - ezcc = Mr.r(erect, 'E', 'Ez') + excc = Mr.aveEx2CC*Mr.r(ecart, 'E', 'Ex') + eycc = Mr.aveEy2CC*Mr.r(ecart, 'E', 'Ey') + ezcc = Mr.aveEz2CC*Mr.r(ecart, 'E', 'Ez') + + indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5]) + indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5]) + + TOL = 1e-2 + assert np.abs(float(excc[indX]) - 0) < TOL + assert np.abs(float(excc[indY]) + 1) < TOL + assert np.abs(float(eycc[indX]) - 1) < TOL + assert np.abs(float(eycc[indY]) - 0) < TOL + assert np.abs(ezcc.sum()) < TOL + + mag = (excc**2 + eycc**2)**0.5 + dist = ((Mr.gridCC[:,0] + 0.2)**2 + (Mr.gridCC[:,1] + 0.2)**2)**0.5 + + assert np.abs(mag[dist > 0.1].max() - 1) < TOL + assert np.abs(mag[dist > 0.1].min() - 1) < TOL + + + def test_getInterpMatCartMesh_Edges2Faces(self): + + Mr = Mesh.TensorMesh([100,100,2], x0='CC0') + Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0]) + + Pe2f = Mc.getInterpolationMatCartMesh(Mr, 'E', locTypeTo='F') + me = np.ones(Mc.nE) + + frect = Pe2f * me + + excc = Mr.aveFx2CC*Mr.r(frect, 'F', 'Fx') + eycc = Mr.aveFy2CC*Mr.r(frect, 'F', 'Fy') + ezcc = Mr.r(frect, 'F', 'Fz') indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5]) indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])