mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-27 17:17:17 +08:00
Merge branch 'master' of https://github.com/simpeg/simpeg into analytics
This commit is contained in:
seogi_macbook
+1
-1
@@ -1,4 +1,4 @@
|
||||
[bumpversion]
|
||||
current_version = 0.1.10
|
||||
current_version = 0.1.11
|
||||
files = setup.py SimPEG/__init__.py docs/conf.py
|
||||
|
||||
|
||||
+157
-196
@@ -1,12 +1,16 @@
|
||||
from SimPEG import np
|
||||
from SimPEG import np, Utils
|
||||
import BaseDC as DC
|
||||
import BaseDC as IP
|
||||
import warnings
|
||||
|
||||
def getActiveindfromTopo(mesh, topo):
|
||||
# def genActiveindfromTopo(mesh, topo):
|
||||
"""
|
||||
Get active indices from topography
|
||||
"""
|
||||
warnings.warn(
|
||||
"`getActiveindfromTopo` is deprecated and will be removed in future versions. Use `SimPEG.Utils.surface2ind_topo` instead",
|
||||
FutureWarning)
|
||||
from scipy.interpolate import NearestNDInterpolator
|
||||
if mesh.dim==3:
|
||||
nCxy = mesh.nCx*mesh.nCy
|
||||
@@ -28,6 +32,9 @@ def gettopoCC(mesh, airind):
|
||||
"""
|
||||
Get topography from active indices of mesh.
|
||||
"""
|
||||
warnings.warn(
|
||||
"`gettopoCC` is deprecated and will be removed in future versions. Use `SimPEG.Utils.surface2ind_topo` instead",
|
||||
FutureWarning)
|
||||
mesh2D = Mesh.TensorMesh([mesh.hx, mesh.hy], mesh.x0[:2])
|
||||
zc = mesh.gridCC[:,2]
|
||||
AIRIND = airind.reshape((mesh.vnC[0]*mesh.vnC[1],mesh.vnC[2]), order='F')
|
||||
@@ -118,34 +125,27 @@ def readUBC_DC3Dobstopo(filename,mesh,topo,probType="CC"):
|
||||
|
||||
def readUBC_DC2DModel(fileName):
|
||||
"""
|
||||
Read UBC GIF 2DTensor model and generate 2D Tensor model in simpeg
|
||||
Read UBC GIF 2DTensor model and generate 2D Tensor model in simpeg
|
||||
|
||||
Input:
|
||||
:param fileName, path to the UBC GIF 2D model file
|
||||
|
||||
Output:
|
||||
:param SimPEG TensorMesh 2D object
|
||||
:return
|
||||
|
||||
Created on Thu Nov 12 13:14:10 2015
|
||||
|
||||
@author: dominiquef
|
||||
:param string fileName: path to the UBC GIF 2D model file
|
||||
:rtype: TensorMesh
|
||||
:return: SimPEG TensorMesh 2D object
|
||||
|
||||
"""
|
||||
from SimPEG import np, mkvc
|
||||
|
||||
# Open fileand skip header... assume that we know the mesh already
|
||||
obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
|
||||
obsfile = np.genfromtxt(fileName, delimiter=' \n', dtype=np.str, comments='!')
|
||||
|
||||
dim = np.array(obsfile[0].split(),dtype=float)
|
||||
dim = np.array(obsfile[0].split(), dtype=float)
|
||||
|
||||
temp = np.array(obsfile[1].split(),dtype=float)
|
||||
temp = np.array(obsfile[1].split(), dtype=float)
|
||||
|
||||
if len(temp) > 1:
|
||||
model = np.zeros(dim)
|
||||
|
||||
for ii in range(len(obsfile)-1):
|
||||
mm = np.array(obsfile[ii+1].split(),dtype=float)
|
||||
mm = np.array(obsfile[ii+1].split(), dtype=float)
|
||||
model[:,ii] = mm
|
||||
|
||||
model = model[:,::-1]
|
||||
@@ -153,10 +153,10 @@ def readUBC_DC2DModel(fileName):
|
||||
else:
|
||||
|
||||
if len(obsfile[1:])==1:
|
||||
mm = np.array(obsfile[1:].split(),dtype=float)
|
||||
mm = np.array(obsfile[1:].split(), dtype=float)
|
||||
|
||||
else:
|
||||
mm = np.array(obsfile[1:],dtype=float)
|
||||
mm = np.array(obsfile[1:], dtype=float)
|
||||
|
||||
# Permute the second dimension to flip the order
|
||||
model = mm.reshape(dim[1],dim[0])
|
||||
@@ -169,23 +169,19 @@ def readUBC_DC2DModel(fileName):
|
||||
|
||||
return model
|
||||
|
||||
def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cblabel=True, axlabel = True, colorbar = True, contour = 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.
|
||||
Read list of 2D tx-rx location and plot a speudo-section of apparent
|
||||
resistivity.
|
||||
|
||||
Assumes flat topo for now...
|
||||
Assumes flat topo for now...
|
||||
|
||||
Input:
|
||||
:param d2D, z0
|
||||
:switch stype -> Either 'pdp' (pole-dipole) | 'dpdp' (dipole-dipole)
|
||||
:switch dtype=-> Either 'appr' (app. res) | 'appc' (app. con) | 'volt' (potential)
|
||||
Output:
|
||||
:figure scatter plot overlayed on image
|
||||
|
||||
Edited Feb 17th, 2016
|
||||
|
||||
@author: dominiquef
|
||||
:param SurveyDC DCsurvey:
|
||||
:param string surveyType: Either 'pole-dipole' | 'dipole-dipole'
|
||||
:param string unitType: Either 'appResistivity' | 'appConductivity' | 'volt'
|
||||
:rtype: matplotlib.plt
|
||||
:return: figure scatter plot overlayed on image
|
||||
|
||||
"""
|
||||
from SimPEG import np
|
||||
@@ -218,39 +214,39 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cbl
|
||||
Cmid = (Tx[0][0] + Tx[1][0])/2
|
||||
Pmid = (Rx[0][:,0] + Rx[1][:,0])/2
|
||||
|
||||
# Change output for dtype
|
||||
if dtype == 'volt':
|
||||
# Change output for unitType
|
||||
if unitType == 'volt':
|
||||
|
||||
rho = np.hstack([rho,data])
|
||||
|
||||
else:
|
||||
|
||||
# Compute pant leg of apparent rho
|
||||
if stype == 'pdp':
|
||||
if surveyType == 'pole-dipole':
|
||||
|
||||
leg = data * 2*np.pi * MA * ( MA + MN ) / MN
|
||||
|
||||
elif stype == 'dpdp':
|
||||
elif surveyType == 'dipole-dipole':
|
||||
|
||||
leg = data * 2*np.pi / ( 1/MA - 1/MB - 1/NB + 1/NA )
|
||||
|
||||
else:
|
||||
print """dtype must be 'pdp'(pole-dipole) | 'dpdp' (dipole-dipole) """
|
||||
print """unitType must be 'pole-dipole' | 'dipole-dipole' """
|
||||
break
|
||||
|
||||
|
||||
if dtype == 'appc':
|
||||
if unitType == 'appConductivity':
|
||||
|
||||
leg = np.log10(abs(1./leg))
|
||||
rho = np.hstack([rho,leg])
|
||||
|
||||
elif dtype == 'appr':
|
||||
elif unitType == 'appResistivity':
|
||||
|
||||
leg = np.log10(abs(leg))
|
||||
rho = np.hstack([rho,leg])
|
||||
|
||||
else:
|
||||
print """dtype must be 'appr' | 'appc' | 'volt' """
|
||||
print """unitType must be 'appResistivity' | 'appConductivity' | 'volt' """
|
||||
break
|
||||
|
||||
midx = np.hstack([midx, ( Cmid + Pmid )/2 ])
|
||||
@@ -259,7 +255,7 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cbl
|
||||
# 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()
|
||||
@@ -268,36 +264,37 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cbl
|
||||
|
||||
# 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, vmin = vmin, vmax = vmax)
|
||||
plt.gca().tick_params(axis='both', which='major', labelsize=8)
|
||||
|
||||
|
||||
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 dtype == 'volt':
|
||||
|
||||
if unitType == 'volt':
|
||||
cbar = plt.colorbar(ph, ax = axs, format="%4.1f",fraction=0.04,orientation="horizontal")
|
||||
|
||||
else:
|
||||
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 cblabel:
|
||||
if dtype == 'appc':
|
||||
cbar.set_label("App.Cond",size=12)
|
||||
elif dtype == 'appr':
|
||||
cbar.set_label("App.Res.",size=12)
|
||||
elif dtype == 'volt':
|
||||
cbar.set_label("Potential (V)",size=12)
|
||||
|
||||
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:
|
||||
@@ -310,27 +307,24 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cbl
|
||||
|
||||
return ph
|
||||
|
||||
def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
|
||||
def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
|
||||
"""
|
||||
Load in endpoints and survey specifications to generate Tx, Rx location
|
||||
stations.
|
||||
Load in endpoints and survey specifications to generate Tx, Rx location
|
||||
stations.
|
||||
|
||||
Assumes flat topo for now...
|
||||
Assumes flat topo for now...
|
||||
|
||||
Input:
|
||||
:param endl -> input endpoints [x1, y1, z1, x2, y2, z2]
|
||||
:object mesh -> SimPEG mesh object
|
||||
:switch stype -> "dpdp" (dipole-dipole) | "pdp" (pole-dipole) | 'gradient'
|
||||
: param a, n -> pole seperation, number of rx dipoles per tx
|
||||
:param numpy.array endl: input endpoints [[x1, y1] , [x2, y2]]
|
||||
:param Mesh mesh: SimPEG mesh object
|
||||
:param string surveyType: 'dipole-dipole' | 'pole-dipole' | 'gradient'
|
||||
:param float AM_sep: transmitter (A) - receiver (M) seperation
|
||||
:param float b: receiver dipole seperation
|
||||
:param float nrx: pole seperation, number of rx dipoles per tx
|
||||
|
||||
Output:
|
||||
:param Tx, Rx -> List objects for each tx location
|
||||
Lines: P1x, P1y, P1z, P2x, P2y, P2z
|
||||
:rtype: DC.Survey, Src, Rx
|
||||
:returns: DC survey, Source
|
||||
|
||||
Created on Wed December 9th, 2015
|
||||
|
||||
@author: dominiquef
|
||||
!! Require clean up to deal with DCsurvey
|
||||
!! Require clean up to deal with DCsurvey
|
||||
"""
|
||||
|
||||
from SimPEG import np
|
||||
@@ -346,17 +340,17 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
|
||||
dl_x = ( endl[1,0] - endl[0,0] ) / dl_len
|
||||
dl_y = ( endl[1,1] - endl[0,1] ) / dl_len
|
||||
|
||||
nstn = np.floor( dl_len / a )
|
||||
nstn = np.floor( dl_len / AM_sep )
|
||||
|
||||
# Compute discrete pole location along line
|
||||
stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*a
|
||||
stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*a
|
||||
stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*AM_sep
|
||||
stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*AM_sep
|
||||
|
||||
# Create line of P1 locations
|
||||
M = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]]
|
||||
|
||||
# Create line of P2 locations
|
||||
N = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]
|
||||
N = np.c_[stn_x+AM_sep*dl_x, stn_y+AM_sep*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]
|
||||
|
||||
## Build list of Tx-Rx locations depending on survey type
|
||||
# Dipole-dipole: Moving tx with [a] spacing -> [AB a MN1 a MN2 ... a MNn]
|
||||
@@ -366,14 +360,14 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
|
||||
SrcList = []
|
||||
|
||||
|
||||
if stype != 'gradient':
|
||||
if surveyType != 'gradient':
|
||||
|
||||
for ii in range(0, int(nstn)-1):
|
||||
|
||||
|
||||
if stype == 'dpdp':
|
||||
if surveyType == 'dipole-dipole':
|
||||
tx = np.c_[M[ii,:],N[ii,:]]
|
||||
elif stype == 'pdp':
|
||||
elif surveyType == 'pole-dipole':
|
||||
tx = np.c_[M[ii,:],M[ii,:]]
|
||||
|
||||
# Rx.append(np.c_[M[ii+1:indx,:],N[ii+1:indx,:]])
|
||||
@@ -382,33 +376,33 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
|
||||
AB = xy_2_r(tx[0,1],endl[1,0],tx[1,1],endl[1,1])
|
||||
|
||||
# Number of receivers to fit
|
||||
nstn = np.min([np.floor( (AB - b) / a ) , n])
|
||||
nstn = np.min([np.floor( (AB - MN_sep) / AM_sep ) , nrx])
|
||||
|
||||
# Check if there is enough space, else break the loop
|
||||
if nstn <= 0:
|
||||
continue
|
||||
|
||||
# Compute discrete pole location along line
|
||||
stn_x = N[ii,0] + dl_x*b + np.array(range(int(nstn)))*dl_x*a
|
||||
stn_y = N[ii,1] + dl_y*b + np.array(range(int(nstn)))*dl_y*a
|
||||
stn_x = N[ii,0] + dl_x*MN_sep + np.array(range(int(nstn)))*dl_x*AM_sep
|
||||
stn_y = N[ii,1] + dl_y*MN_sep + np.array(range(int(nstn)))*dl_y*AM_sep
|
||||
|
||||
# Create receiver poles
|
||||
# Create line of P1 locations
|
||||
P1 = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]]
|
||||
|
||||
# Create line of P2 locations
|
||||
P2 = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]
|
||||
P2 = np.c_[stn_x+AM_sep*dl_x, stn_y+AM_sep*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]
|
||||
|
||||
Rx.append(np.c_[P1,P2])
|
||||
rxClass = DC.RxDipole(P1, P2)
|
||||
Tx.append(tx)
|
||||
if stype == 'dpdp':
|
||||
if surveyType == 'dipole-dipole':
|
||||
srcClass = DC.SrcDipole([rxClass], M[ii,:],N[ii,:])
|
||||
elif stype == 'pdp':
|
||||
elif surveyType == 'pole-dipole':
|
||||
srcClass = DC.SrcDipole([rxClass], M[ii,:],M[ii,:])
|
||||
SrcList.append(srcClass)
|
||||
|
||||
elif stype == 'gradient':
|
||||
elif surveyType == 'gradient':
|
||||
|
||||
# Gradient survey only requires Tx at end of line and creates a square
|
||||
# grid of receivers at in the middle at a pre-set minimum distance
|
||||
@@ -416,23 +410,23 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
|
||||
Tx.append(np.c_[M[0,:],N[-1,:]])
|
||||
|
||||
# Get the edge limit of survey area
|
||||
min_x = endl[0,0] + dl_x * b
|
||||
min_y = endl[0,1] + dl_y * b
|
||||
min_x = endl[0,0] + dl_x * MN_sep
|
||||
min_y = endl[0,1] + dl_y * MN_sep
|
||||
|
||||
max_x = endl[1,0] - dl_x * b
|
||||
max_y = endl[1,1] - dl_y * b
|
||||
max_x = endl[1,0] - dl_x * MN_sep
|
||||
max_y = endl[1,1] - dl_y * MN_sep
|
||||
|
||||
box_l = np.sqrt( (min_x - max_x)**2 + (min_y - max_y)**2 )
|
||||
box_w = box_l/2.
|
||||
|
||||
nstn = np.floor( box_l / a )
|
||||
nstn = np.floor( box_l / AM_sep )
|
||||
|
||||
# Compute discrete pole location along line
|
||||
stn_x = min_x + np.array(range(int(nstn)))*dl_x*a
|
||||
stn_y = min_y + np.array(range(int(nstn)))*dl_y*a
|
||||
stn_x = min_x + np.array(range(int(nstn)))*dl_x*AM_sep
|
||||
stn_y = min_y + np.array(range(int(nstn)))*dl_y*AM_sep
|
||||
|
||||
# Define number of cross lines
|
||||
nlin = int(np.floor( box_w / a ))
|
||||
nlin = int(np.floor( box_w / AM_sep ))
|
||||
lind = range(-nlin,nlin+1)
|
||||
|
||||
ngrad = nstn * len(lind)
|
||||
@@ -441,12 +435,12 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
|
||||
for ii in range( len(lind) ):
|
||||
|
||||
# Move line in perpendicular direction by dipole spacing
|
||||
lxx = stn_x - lind[ii]*a*dl_y
|
||||
lyy = stn_y + lind[ii]*a*dl_x
|
||||
lxx = stn_x - lind[ii]*AM_sep*dl_y
|
||||
lyy = stn_y + lind[ii]*AM_sep*dl_x
|
||||
|
||||
|
||||
M = np.c_[ lxx, lyy , np.ones(nstn).T*mesh.vectorNz[-1]]
|
||||
N = np.c_[ lxx+a*dl_x, lyy+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]
|
||||
N = np.c_[ lxx+AM_sep*dl_x, lyy+AM_sep*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]
|
||||
|
||||
rx[(ii*nstn):((ii+1)*nstn),:] = np.c_[M,N]
|
||||
|
||||
@@ -455,44 +449,38 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
|
||||
srcClass = DC.SrcDipole([rxClass], M[0,:], N[-1,:])
|
||||
SrcList.append(srcClass)
|
||||
else:
|
||||
print """stype must be either 'pdp', 'dpdp' or 'gradient'. """
|
||||
print """surveyType must be either 'pole-dipole', 'dipole-dipole' or 'gradient'. """
|
||||
|
||||
survey = DC.SurveyDC(SrcList)
|
||||
return survey, Tx, Rx
|
||||
|
||||
def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0):
|
||||
|
||||
def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0):
|
||||
"""
|
||||
Write UBC GIF DCIP 2D or 3D observation file
|
||||
|
||||
Input:
|
||||
:string fileName -> including path where the file is written out
|
||||
:DCsurvey DC survey class object
|
||||
:string dtype -> either '2D' | '3D'
|
||||
:string stype -> either 'SURFACE' | 'GENERAL'
|
||||
|
||||
Output:
|
||||
:param UBC2D-Data file
|
||||
:return
|
||||
|
||||
Last edit: February 16th, 2016
|
||||
|
||||
@author: dominiquef
|
||||
|
||||
:param string fileName: including path where the file is written out
|
||||
:param Survey DCsurvey: DC survey class object
|
||||
:param string dim: either '2D' | '3D'
|
||||
:param string surveyType: either 'SURFACE' | 'GENERAL'
|
||||
:rtype: file
|
||||
:return: UBC2D-Data file
|
||||
"""
|
||||
|
||||
from SimPEG import mkvc
|
||||
|
||||
assert (dtype=='2D') | (dtype=='3D'), "Data must be either '2D' | '3D'"
|
||||
assert (stype=='SURFACE') | (stype=='GENERAL') | (stype=='SIMPLE'), "Data must be either 'SURFACE' | 'GENERAL' | 'SIMPLE'"
|
||||
|
||||
assert (dim=='2D') | (dim=='3D'), "Data must be either '2D' | '3D'"
|
||||
assert (surveyType=='SURFACE') | (surveyType=='GENERAL') | (surveyType=='SIMPLE'), "Data must be either 'SURFACE' | 'GENERAL' | 'SIMPLE'"
|
||||
|
||||
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):
|
||||
@@ -506,33 +494,33 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0):
|
||||
M = rx[0]
|
||||
N = rx[1]
|
||||
|
||||
# Adapt source-receiver location for dtype and stype
|
||||
if dtype=='2D':
|
||||
# Adapt source-receiver location for dim and surveyType
|
||||
if dim=='2D':
|
||||
|
||||
if stype == 'SIMPLE':
|
||||
if surveyType == 'SIMPLE':
|
||||
|
||||
#fid.writelines("%e " % ii for ii in mkvc(tx[0,:]))
|
||||
A = np.repeat(tx[0,0],M.shape[0],axis=0)
|
||||
B = np.repeat(tx[0,1],M.shape[0],axis=0)
|
||||
M = M[:,0]
|
||||
N = N[:,0]
|
||||
|
||||
|
||||
np.savetxt(fid, np.c_[A, B, M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
|
||||
|
||||
|
||||
else:
|
||||
|
||||
if stype == 'SURFACE':
|
||||
if surveyType == 'SURFACE':
|
||||
|
||||
fid.writelines("%f " % ii for ii in mkvc(tx[0,:]))
|
||||
M = M[:,0]
|
||||
N = N[:,0]
|
||||
|
||||
if stype == 'GENERAL':
|
||||
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]
|
||||
@@ -540,31 +528,31 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0):
|
||||
# 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='%f',delimiter=' ',newline='\n')
|
||||
|
||||
if dtype=='3D':
|
||||
if dim=='3D':
|
||||
|
||||
if stype == 'SURFACE':
|
||||
if surveyType == 'SURFACE':
|
||||
|
||||
fid.writelines("%e " % ii for ii in mkvc(tx[0:2,:]))
|
||||
M = M[:,0:2]
|
||||
N = N[:,0:2]
|
||||
|
||||
if stype == 'GENERAL':
|
||||
if surveyType == 'GENERAL':
|
||||
|
||||
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
|
||||
|
||||
fid.close()
|
||||
|
||||
def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'):
|
||||
def convertObs_DC3D_to_2D(DCsurvey, lineID, flag='local'):
|
||||
"""
|
||||
Read DC survey and projects the coordinate system
|
||||
according to the flag = 'Xloc' | 'Yloc' | 'local' (default)
|
||||
@@ -573,15 +561,9 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'):
|
||||
|
||||
The Z value is preserved, but Y coordinates zeroed.
|
||||
|
||||
Input:
|
||||
:param survey3D
|
||||
|
||||
Output:
|
||||
:figure survey2D
|
||||
|
||||
Edited April 6th, 2016
|
||||
|
||||
@author: dominiquef
|
||||
:param DC.Survey survey3D: 3D simpeg DC survey
|
||||
:rtype: DC.Survey
|
||||
:return: survey2D
|
||||
|
||||
"""
|
||||
from SimPEG import np
|
||||
@@ -666,39 +648,34 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'):
|
||||
DCsurvey2D.std = np.asarray(DCsurvey.std)
|
||||
|
||||
return DCsurvey2D
|
||||
|
||||
def readUBC_DC3Dobs(fileName, dtype = 'DC'):
|
||||
|
||||
def readUBC_DC3Dobs(fileName, rtype = 'DC'):
|
||||
"""
|
||||
Read UBC GIF IP 3D observation file and generate survey
|
||||
|
||||
Input:
|
||||
:param fileName, path to the UBC GIF 3D obs file
|
||||
|
||||
Output:
|
||||
:param IPsurvey
|
||||
:return
|
||||
|
||||
@author: dominiquef
|
||||
:param string fileName:, path to the UBC GIF 3D obs file
|
||||
:rtype: Survey
|
||||
:return: DCIPsurvey
|
||||
|
||||
"""
|
||||
zflag = True # Flag for z value provided
|
||||
|
||||
|
||||
# Load file
|
||||
if dtype == 'IP':
|
||||
if rtype == 'IP':
|
||||
obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='IPTYPE')
|
||||
|
||||
elif dtype == 'DC':
|
||||
|
||||
elif rtype == 'DC':
|
||||
obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
|
||||
|
||||
|
||||
else:
|
||||
print "dtype must be 'DC'(default) | 'IP'"
|
||||
|
||||
print "rtype must be 'DC'(default) | 'IP'"
|
||||
|
||||
# Pre-allocate
|
||||
srcLists = []
|
||||
Rx = []
|
||||
d = []
|
||||
wd = []
|
||||
|
||||
|
||||
|
||||
# Countdown for number of obs/tx
|
||||
count = 0
|
||||
@@ -717,7 +694,7 @@ def readUBC_DC3Dobs(fileName, dtype = 'DC'):
|
||||
# Check if z value is provided, if False -> nan
|
||||
if len(temp)==5:
|
||||
tx = np.r_[temp[0:2],np.nan,temp[2:4],np.nan]
|
||||
|
||||
|
||||
zflag = False # Pass on the flag to the receiver loc
|
||||
|
||||
else:
|
||||
@@ -729,12 +706,12 @@ def readUBC_DC3Dobs(fileName, dtype = 'DC'):
|
||||
temp = np.fromstring(obsfile[ii], dtype=float,sep=' ') # Get the string
|
||||
|
||||
# Filter out negative IP
|
||||
# if temp[-2] < 0:
|
||||
# if temp[-2] < 0:
|
||||
# count = count -1
|
||||
# print "Negative!"
|
||||
#
|
||||
#
|
||||
# else:
|
||||
|
||||
|
||||
# If the Z-location is provided, otherwise put nan
|
||||
if zflag:
|
||||
|
||||
@@ -772,17 +749,9 @@ def readUBC_DC2Dobs(fileName):
|
||||
------- NEEDS TO BE UPDATED ------
|
||||
Read UBC GIF 2D observation file and generate arrays for tx-rx location
|
||||
|
||||
Input:
|
||||
:param fileName, path to the UBC GIF 2D model file
|
||||
|
||||
Output:
|
||||
:param rx, tx
|
||||
:return
|
||||
|
||||
Created on Thu Nov 12 13:14:10 2015
|
||||
|
||||
@author: dominiquef
|
||||
|
||||
:param string fileName: path to the UBC GIF 2D model file
|
||||
:rtype: (DC.Src, DC.Rx, ??, ??)
|
||||
:return: source_locs, rx_locs, ??, ??
|
||||
"""
|
||||
|
||||
from SimPEG import np
|
||||
@@ -822,11 +791,9 @@ def readUBC_DC2Dpre(fileName):
|
||||
Read UBC GIF DCIP 2D observation file and generate arrays for tx-rx location
|
||||
|
||||
Input:
|
||||
:param fileName, path to the UBC GIF 3D obs file
|
||||
|
||||
Output:
|
||||
DCsurvey
|
||||
:return
|
||||
:param string fileName: path to the UBC GIF 3D obs file
|
||||
:rtype: DC.Survey
|
||||
:return: DCsurvey
|
||||
|
||||
Created on Mon March 9th, 2016 << Doug's 70th Birthday !! >>
|
||||
|
||||
@@ -888,12 +855,9 @@ def readUBC_DC2DMesh(fileName):
|
||||
"""
|
||||
Read UBC GIF 2DTensor mesh and generate 2D Tensor mesh in simpeg
|
||||
|
||||
Input:
|
||||
:param fileName, path to the UBC GIF mesh file
|
||||
|
||||
Output:
|
||||
:param SimPEG TensorMesh 2D object
|
||||
:return
|
||||
:param string fileName: path to the UBC GIF mesh file
|
||||
:rtype: Mesh.TensorMesh
|
||||
:return: SimPEG TensorMesh 2D object
|
||||
|
||||
Created on Thu Nov 12 13:14:10 2015
|
||||
|
||||
@@ -959,12 +923,9 @@ def xy_2_lineID(DCsurvey):
|
||||
they were collected. May need to generalize for random
|
||||
point locations, but will be more expensive
|
||||
|
||||
Input:
|
||||
:param DCdict Vectors of station location
|
||||
|
||||
Output:
|
||||
:param LineID Vector of integers
|
||||
:return
|
||||
:param numpy.array DCdict: Vectors of station location
|
||||
:rtype: numpy.array
|
||||
:return: LineID Vector of integers
|
||||
|
||||
Created on Thu Feb 11, 2015
|
||||
|
||||
|
||||
+115
-54
@@ -144,6 +144,7 @@ class BetaSchedule(InversionDirective):
|
||||
if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
|
||||
self.invProb.beta /= self.coolingFactor
|
||||
|
||||
|
||||
class TargetMisfit(InversionDirective):
|
||||
|
||||
chifact = 1.
|
||||
@@ -242,12 +243,6 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
|
||||
# Save the file as a npz
|
||||
np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
|
||||
|
||||
|
||||
# class UpdateReferenceModel(Parameter):
|
||||
|
||||
# mref0 = None
|
||||
|
||||
# def nextIter(self):
|
||||
# mref = getattr(self, 'm_prev', None)
|
||||
# if mref is None:
|
||||
# if self.debug: print 'UpdateReferenceModel is using mref0'
|
||||
@@ -258,56 +253,138 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
|
||||
class Update_IRLS(InversionDirective):
|
||||
|
||||
eps_min = None
|
||||
eps_p = None
|
||||
eps_q = None
|
||||
norms = [2.,2.,2.,2.]
|
||||
factor = None
|
||||
gamma = None
|
||||
phi_m_last = None
|
||||
phi_d_last = None
|
||||
f_old = None
|
||||
f_min_change = 1e-2
|
||||
beta_tol = 5e-2
|
||||
|
||||
# Solving parameter for IRLS (mode:2)
|
||||
IRLSiter = 0
|
||||
minGNiter = 5
|
||||
maxIRLSiter = 10
|
||||
iterStart = 0
|
||||
|
||||
# Beta schedule
|
||||
coolingFactor = 2.
|
||||
coolingRate = 1
|
||||
|
||||
mode = 1
|
||||
|
||||
@property
|
||||
def target(self):
|
||||
if getattr(self, '_target', None) is None:
|
||||
self._target = self.survey.nD*0.5
|
||||
return self._target
|
||||
@target.setter
|
||||
def target(self, val):
|
||||
self._target = val
|
||||
|
||||
def initialize(self):
|
||||
|
||||
# Scale the regularization for changes in norm
|
||||
if getattr(self, 'phi_m_last', None) is not None:
|
||||
|
||||
self.reg.curModel = self.invProb.curModel
|
||||
self.reg.gamma = 1.
|
||||
phim_new = self.reg.eval(self.invProb.curModel)
|
||||
self.gamma = self.phi_m_last / phim_new
|
||||
|
||||
self.reg.curModel = self.invProb.curModel
|
||||
self.reg.gamma = self.gamma
|
||||
|
||||
if getattr(self, 'phi_d_last', None) is None:
|
||||
self.phi_d_last = self.invProb.phi_d
|
||||
if self.mode == 1:
|
||||
self.reg.norms = [2., 2., 2., 2.]
|
||||
|
||||
def endIter(self):
|
||||
# Cool the threshold parameter if required
|
||||
if getattr(self, 'factor', None) is not None:
|
||||
eps = self.reg.eps / self.factor
|
||||
|
||||
if getattr(self, 'eps_min', None) is not None:
|
||||
self.reg.eps = np.max([self.eps_min,eps])
|
||||
# After reaching target misfit with l2-norm, switch to IRLS (mode:2)
|
||||
if self.invProb.phi_d < self.target and self.mode == 1:
|
||||
print "Convergence with smooth l2-norm regularization: Start IRLS steps..."
|
||||
|
||||
self.mode = 2
|
||||
print self.eps_p, self.eps_q, self.norms
|
||||
self.reg.eps_p = self.eps_p
|
||||
self.reg.eps_q = self.eps_q
|
||||
self.reg.norms = self.norms
|
||||
self.coolingFactor = 1.
|
||||
self.coolingRate = 1
|
||||
self.iterStart = self.opt.iter
|
||||
self.phi_d_last = self.invProb.phi_d
|
||||
self.phi_m_last = self.invProb.phi_m_last
|
||||
|
||||
self.reg.l2model = self.invProb.curModel
|
||||
self.reg.curModel = self.invProb.curModel
|
||||
|
||||
if getattr(self, 'f_old', None) is None:
|
||||
self.f_old = self.reg.eval(self.invProb.curModel)#self.invProb.evalFunction(self.invProb.curModel, return_g=False, return_H=False)
|
||||
|
||||
# Beta Schedule
|
||||
if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0:
|
||||
if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
|
||||
self.invProb.beta /= self.coolingFactor
|
||||
|
||||
|
||||
# Only update after GN iterations
|
||||
if (self.opt.iter-self.iterStart) % self.minGNiter == 0 and self.mode==2:
|
||||
|
||||
self.IRLSiter += 1
|
||||
|
||||
phim_new = self.reg.eval(self.invProb.curModel)
|
||||
self.f_change = np.abs(self.f_old - phim_new) / self.f_old
|
||||
|
||||
print "Regularization decrease: %6.3e" % (self.f_change)
|
||||
|
||||
# Check for maximum number of IRLS cycles
|
||||
if self.IRLSiter == self.maxIRLSiter:
|
||||
print "Reach maximum number of IRLS cycles: %i" % self.maxIRLSiter
|
||||
self.opt.stopNextIteration = True
|
||||
return
|
||||
|
||||
# Check if the function has changed enough
|
||||
if self.f_change < self.f_min_change and self.IRLSiter > 1:
|
||||
print "Minimum decrease in regularization. End of IRLS"
|
||||
self.opt.stopNextIteration = True
|
||||
return
|
||||
else:
|
||||
self.reg.eps = eps
|
||||
self.f_old = phim_new
|
||||
|
||||
# Get phi_m at the end of current iteration
|
||||
self.phi_m_last = self.invProb.phi_m_last
|
||||
# Cool the threshold parameter if required
|
||||
if getattr(self, 'factor', None) is not None:
|
||||
eps = self.reg.eps / self.factor
|
||||
|
||||
# Update the model used for the IRLS weights
|
||||
self.reg.curModel = self.invProb.curModel
|
||||
if getattr(self, 'eps_min', None) is not None:
|
||||
self.reg.eps = np.max([self.eps_min,eps])
|
||||
else:
|
||||
self.reg.eps = eps
|
||||
|
||||
# Temporarely set gamma to 1. to get raw phi_m
|
||||
self.reg.gamma = 1.
|
||||
# Get phi_m at the end of current iteration
|
||||
self.phi_m_last = self.invProb.phi_m_last
|
||||
|
||||
# Compute new model objective function value
|
||||
phim_new = self.reg.eval(self.invProb.curModel)
|
||||
# Reset the regularization matrices so that it is
|
||||
# recalculated for current model
|
||||
self.reg._Wsmall = None
|
||||
self.reg._Wx = None
|
||||
self.reg._Wy = None
|
||||
self.reg._Wz = None
|
||||
|
||||
# Update gamma to scale the regularization between IRLS iterations
|
||||
self.reg.gamma = self.phi_m_last / phim_new
|
||||
# Update the model used for the IRLS weights
|
||||
self.reg.curModel = self.invProb.curModel
|
||||
|
||||
# Set the weighting matrix to None so that it is recomputed next time
|
||||
# it is called in the inversion
|
||||
self.reg._W = None
|
||||
# Temporarely set gamma to 1. to get raw phi_m
|
||||
self.reg.gamma = 1.
|
||||
|
||||
# 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
|
||||
|
||||
# Reset the regularization matrices again for new gamma
|
||||
self.reg._Wsmall = None
|
||||
self.reg._Wx = None
|
||||
self.reg._Wy = None
|
||||
self.reg._Wz = None
|
||||
|
||||
# Check if misfit is within the tolerance, otherwise scale beta
|
||||
val = self.invProb.phi_d / (self.survey.nD*0.5)
|
||||
|
||||
if np.abs(1.-val) > self.beta_tol:
|
||||
self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d
|
||||
|
||||
class Update_lin_PreCond(InversionDirective):
|
||||
"""
|
||||
@@ -360,19 +437,3 @@ 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
|
||||
|
||||
@@ -2,7 +2,7 @@ from SimPEG import Mesh, Utils, np, sp
|
||||
import SimPEG.DCIP as DC
|
||||
import time
|
||||
|
||||
def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', plotIt=True):
|
||||
def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', unitType='appConductivity', plotIt=True):
|
||||
"""
|
||||
DC Forward Simulation
|
||||
=====================
|
||||
@@ -15,14 +15,14 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p
|
||||
loc = Location of spheres [[x1,y1,z1],[x2,y2,z2]]
|
||||
radi = Radius of spheres [r1,r2]
|
||||
param = Conductivity of background and two spheres [m0,m1,m2]
|
||||
stype = survey type "pdp" (pole dipole) or "dpdp" (dipole dipole)
|
||||
dtype = Data type "appr" (app res) | "appc" (app cond) | "volt" (potential)
|
||||
surveyType = survey type 'pole-dipole' or 'dipole-dipole'
|
||||
unitType = Data type "appResistivity" | "appConductivity" | "volt"
|
||||
Created by @fourndo
|
||||
|
||||
"""
|
||||
|
||||
assert stype in ['pdp', 'dpdp'], "Source type (stype) must be pdp or dpdp (pole dipole or dipole dipole)"
|
||||
assert dtype in ['appr', 'appc', 'volt'], "Data type (dtype) must be appr (app res) or appc (app cond) or volt (potential)"
|
||||
assert surveyType in ['pole-dipole', 'dipole-dipole'], "Source type (surveyType) must be pdp or dpdp (pole dipole or dipole dipole)"
|
||||
assert unitType in ['appResistivity', 'appConductivity', 'volt'], "Unit type (unitType) must be appResistivity or appConductivity or volt (potential)"
|
||||
|
||||
if loc is None:
|
||||
loc = np.c_[[-50.,0.,-50.],[50.,0.,-50.]]
|
||||
@@ -73,8 +73,8 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p
|
||||
locs = np.c_[mesh.gridCC[indx,0],mesh.gridCC[indx,1],np.ones(2).T*mesh.vectorNz[-1]]
|
||||
|
||||
# We will handle the geometry of the survey for you and create all the combination of tx-rx along line
|
||||
# [Tx, Rx] = DC.gen_DCIPsurvey(locs, mesh, stype, param[0], param[1], param[2])
|
||||
survey, Tx, Rx = DC.gen_DCIPsurvey(locs, mesh, stype, param[0], param[1], param[2])
|
||||
# [Tx, Rx] = DC.gen_DCIPsurvey(locs, mesh, surveyType, param[0], param[1], param[2])
|
||||
survey, Tx, Rx = DC.gen_DCIPsurvey(locs, mesh, surveyType, param[0], param[1], param[2])
|
||||
|
||||
# Define some global geometry
|
||||
dl_len = np.sqrt( np.sum((locs[0,:] - locs[1,:])**2) )
|
||||
@@ -118,8 +118,8 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p
|
||||
rxloc_N = np.asarray(Rx[ii][:,3:])
|
||||
|
||||
|
||||
# For usual cases "dpdp" or "gradient"
|
||||
if stype == 'pdp':
|
||||
# For usual cases 'dipole-dipole' or "gradient"
|
||||
if surveyType == 'pole-dipole':
|
||||
# Create an "inifinity" pole
|
||||
tx = np.squeeze(Tx[ii][:,0:1])
|
||||
tinf = tx + np.array([dl_x,dl_y,0])*dl_len*2
|
||||
@@ -157,12 +157,12 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p
|
||||
fig = plt.figure(figsize=(7,7))
|
||||
ax = plt.subplot(2,1,1, aspect='equal')
|
||||
# Plot the location of the spheres for reference
|
||||
circle1=plt.Circle((loc[0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3)
|
||||
circle2=plt.Circle((loc[0,1],loc[2,1]),radi[1],color='k',fill=False, lw=3)
|
||||
circle1=plt.Circle((loc[0,0], loc[2,0]), radi[0], color='w', fill=False, lw=3)
|
||||
circle2=plt.Circle((loc[0,1], loc[2,1]), radi[1], color='k', fill=False, lw=3)
|
||||
ax.add_artist(circle1)
|
||||
ax.add_artist(circle2)
|
||||
|
||||
dat = mesh.plotSlice(np.log10(model), ax =ax, normal = 'Y',
|
||||
dat = mesh.plotSlice(np.log10(model), ax = ax, normal = 'Y',
|
||||
ind = indy,grid=True, clim = np.log10([sig.min(),sig.max()]))
|
||||
|
||||
ax.set_title('3-D model')
|
||||
@@ -188,15 +188,13 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p
|
||||
ax2 = plt.subplot(2,1,2, aspect='equal')
|
||||
|
||||
# Plot the location of the spheres for reference
|
||||
circle1=plt.Circle((loc[0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3)
|
||||
circle2=plt.Circle((loc[0,1],loc[2,1]),radi[1],color='k',fill=False, lw=3)
|
||||
circle1=plt.Circle((loc[0,0], loc[2,0]), radi[0], color='w', fill=False, lw=3)
|
||||
circle2=plt.Circle((loc[0,1], loc[2,1]), radi[1], color='k', fill=False, lw=3)
|
||||
ax2.add_artist(circle1)
|
||||
ax2.add_artist(circle2)
|
||||
|
||||
# Add the speudo section
|
||||
dat = DC.plot_pseudoSection(survey2D,ax2,stype=stype, dtype = dtype)
|
||||
|
||||
# plt.scatter(Tx2d[0][:],Tx[0][2,:],s=40,c='g', marker='v')
|
||||
dat = DC.plot_pseudoSection(survey2D, ax2, surveyType=surveyType, unitType=unitType) # plt.scatter(Tx2d[0][:],Tx[0][2,:],s=40,c='g', marker='v')
|
||||
# plt.scatter(Rx2d[0][:],Rx[0][:,2::3],s=40,c='y')
|
||||
# plt.plot(np.r_[Tx2d[0][0],Rx2d[-1][-1,-1]],np.ones(2)*mesh.vectorNz[-1], color='k')
|
||||
ax2.set_title('Apparent Conductivity data')
|
||||
|
||||
@@ -1,7 +1,7 @@
|
||||
from SimPEG import *
|
||||
|
||||
|
||||
def run(N=200, plotIt=True):
|
||||
def run(N=100, plotIt=True):
|
||||
"""
|
||||
Inversion: Linear Problem
|
||||
=========================
|
||||
@@ -18,6 +18,8 @@ def run(N=200, plotIt=True):
|
||||
mesh = Mesh.TensorMesh([N])
|
||||
|
||||
m0 = np.ones(mesh.nC) * 1e-4
|
||||
mref = np.zeros(mesh.nC)
|
||||
|
||||
nk = 10
|
||||
jk = np.linspace(1.,nk,nk)
|
||||
p = -2.
|
||||
@@ -50,57 +52,47 @@ def run(N=200, plotIt=True):
|
||||
wr = np.sum(prob.G**2.,axis=0)**0.5
|
||||
wr = ( wr/np.max(wr) )
|
||||
|
||||
reg = Regularization.Simple(mesh)
|
||||
reg.wght = wr
|
||||
|
||||
# reg = Regularization.Simple(mesh)
|
||||
# reg.mref = mref
|
||||
# reg.cell_weights = wr
|
||||
#
|
||||
dmis = DataMisfit.l2_DataMisfit(survey)
|
||||
dmis.Wd = 1./wd
|
||||
|
||||
opt = Optimization.ProjectedGNCG(maxIter=30,lower=-2.,upper=2., maxIterCG= 20, tolCG = 1e-4)
|
||||
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
|
||||
invProb.curModel = m0
|
||||
|
||||
beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1)
|
||||
target = Directives.TargetMisfit()
|
||||
|
||||
#
|
||||
# opt = Optimization.ProjectedGNCG(maxIter=20,lower=-2.,upper=2., maxIterCG= 10, tolCG = 1e-4)
|
||||
# invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
|
||||
# invProb.curModel = m0
|
||||
#
|
||||
# beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1)
|
||||
# target = Directives.TargetMisfit()
|
||||
#
|
||||
betaest = Directives.BetaEstimate_ByEig()
|
||||
inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target])
|
||||
|
||||
|
||||
mrec = inv.run(m0)
|
||||
ml2 = mrec
|
||||
print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
|
||||
|
||||
# Switch regularization to sparse
|
||||
phim = invProb.phi_m_last
|
||||
phid = invProb.phi_d
|
||||
# inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target])
|
||||
#
|
||||
#
|
||||
# mrec = inv.run(m0)
|
||||
# ml2 = mrec
|
||||
# print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
|
||||
#
|
||||
# # Switch regularization to sparse
|
||||
# phim = invProb.phi_m_last
|
||||
# phid = invProb.phi_d
|
||||
|
||||
reg = Regularization.Sparse(mesh)
|
||||
reg.mref = mref
|
||||
reg.cell_weights = wr
|
||||
|
||||
#==============================================================================
|
||||
# fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
|
||||
# dmdx = reg.mesh.cellDiffxStencil * mrec
|
||||
# plt.plot(np.sort(dmdx))
|
||||
#==============================================================================
|
||||
|
||||
#reg.recModel = mrec
|
||||
reg.wght = np.ones(mesh.nC)
|
||||
reg.mref = np.zeros(mesh.nC)
|
||||
reg.eps_p = 5e-2
|
||||
reg.eps_q = 1e-2
|
||||
reg.norms = [0., 0., 2., 2.]
|
||||
reg.wght = wr
|
||||
eps_p = 5e-2
|
||||
eps_q = 5e-2
|
||||
norms = [0., 0., 2., 2.]
|
||||
|
||||
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()
|
||||
target = Directives.TargetMisfit()
|
||||
IRLS =Directives.Update_IRLS( phi_m_last = phim, phi_d_last = phid )
|
||||
opt = Optimization.ProjectedGNCG(maxIter=100 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3)
|
||||
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
|
||||
update_Jacobi = Directives.Update_lin_PreCond()
|
||||
IRLS = Directives.Update_IRLS( norms=norms, eps_p=eps_p, eps_q=eps_q)
|
||||
|
||||
inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS])
|
||||
|
||||
m0 = mrec
|
||||
inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest,update_Jacobi])
|
||||
|
||||
# Run inversion
|
||||
mrec = inv.run(m0)
|
||||
@@ -117,7 +109,7 @@ def run(N=200, plotIt=True):
|
||||
axes[0].set_title('Columns of matrix G')
|
||||
|
||||
axes[1].plot(mesh.vectorCCx, mtrue, 'b-')
|
||||
axes[1].plot(mesh.vectorCCx, ml2, 'r-')
|
||||
axes[1].plot(mesh.vectorCCx, reg.l2model, 'r-')
|
||||
#axes[1].legend(('True Model', 'Recovered Model'))
|
||||
axes[1].set_ylim(-1.0,1.25)
|
||||
|
||||
|
||||
+12
-31
@@ -1,22 +1,25 @@
|
||||
from SimPEG import Mesh, Utils, np, SolverLU
|
||||
|
||||
## 2D DC forward modeling example with Tensor and Curvilinear Meshes
|
||||
|
||||
def run(plotIt=True):
|
||||
|
||||
"""
|
||||
Mesh: Basic Forward 2D DC Resistivity
|
||||
=====================================
|
||||
|
||||
2D DC forward modeling example with Tensor and Curvilinear Meshes
|
||||
"""
|
||||
|
||||
# Step1: Generate Tensor and Curvilinear Mesh
|
||||
sz = [40,40]
|
||||
# Tensor Mesh
|
||||
tM = Mesh.TensorMesh(sz)
|
||||
# Curvilinear Mesh
|
||||
rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))
|
||||
|
||||
# Step2: Direct Current (DC) operator
|
||||
def DCfun(mesh, pts):
|
||||
D = mesh.faceDiv
|
||||
G = D.T
|
||||
sigma = 1e-2*np.ones(mesh.nC)
|
||||
Msigi = mesh.getFaceInnerProduct(1./sigma)
|
||||
MsigI = Utils.sdInv(Msigi)
|
||||
A = D*MsigI*G
|
||||
MsigI = mesh.getFaceInnerProduct(sigma, invProp=True, invMat=True)
|
||||
A = -D*MsigI*D.T
|
||||
A[-1,-1] /= mesh.vol[-1] # Remove null space
|
||||
rhs = np.zeros(mesh.nC)
|
||||
txind = Utils.meshutils.closestPoints(mesh, pts)
|
||||
@@ -37,39 +40,17 @@ def run(plotIt=True):
|
||||
if not plotIt: return
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
import matplotlib
|
||||
from matplotlib.mlab import griddata
|
||||
|
||||
#Step4: Making Figure
|
||||
fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
|
||||
label = ["(a)", "(b)"]
|
||||
opts = {}
|
||||
vmin, vmax = phitM.min(), phitM.max()
|
||||
dat = tM.plotImage(phitM, ax=axes[0], clim=(vmin, vmax), grid=True)
|
||||
|
||||
#TODO: At the moment Curvilinear Mesh do not have plotimage
|
||||
|
||||
Xi = tM.gridCC[:,0].reshape(sz[0], sz[1], order='F')
|
||||
Yi = tM.gridCC[:,1].reshape(sz[0], sz[1], order='F')
|
||||
PHIrM = griddata(rM.gridCC[:,0], rM.gridCC[:,1], phirM, Xi, Yi, interp='linear')
|
||||
axes[1].contourf(Xi, Yi, PHIrM, 100, vmin=vmin, vmax=vmax)
|
||||
|
||||
dat = rM.plotImage(phirM, ax=axes[1], clim=(vmin, vmax), grid=True)
|
||||
cb = plt.colorbar(dat[0], ax=axes[0]); cb.set_label("Voltage (V)")
|
||||
cb = plt.colorbar(dat[0], ax=axes[1]); cb.set_label("Voltage (V)")
|
||||
|
||||
tM.plotGrid(ax=axes[0], **opts)
|
||||
axes[0].set_title('TensorMesh')
|
||||
rM.plotGrid(ax=axes[1], **opts)
|
||||
axes[1].set_title('CurvilinearMesh')
|
||||
for i in range(2):
|
||||
axes[i].set_xlim(0.025, 0.975)
|
||||
axes[i].set_ylim(0.025, 0.975)
|
||||
axes[i].text(0., 1.0, label[i], fontsize=20)
|
||||
if i==0:
|
||||
axes[i].set_ylabel("y")
|
||||
else:
|
||||
axes[i].set_ylabel(" ")
|
||||
axes[i].set_xlabel("x")
|
||||
plt.show()
|
||||
|
||||
|
||||
@@ -0,0 +1,41 @@
|
||||
from SimPEG import *
|
||||
from SimPEG.Utils import surface2ind_topo
|
||||
|
||||
|
||||
def run(plotIt=False, nx = 5, ny = 5):
|
||||
"""
|
||||
Here we show how to use :code:`Utils.surface2ind_topo` to identify cells below
|
||||
a topographic surface.
|
||||
|
||||
"""
|
||||
|
||||
mesh = Mesh.TensorMesh([nx,ny], x0='CC') # 2D mesh
|
||||
xtopo = np.linspace(mesh.gridN[:,0].min(), mesh.gridN[:,0].max())
|
||||
topo = 0.4*np.sin(xtopo*5) # define a topographic surface
|
||||
|
||||
Topo = np.hstack([Utils.mkvc(xtopo,2),Utils.mkvc(topo,2)]) #make it an array
|
||||
|
||||
indcc = surface2ind_topo(mesh, Topo,'CC')
|
||||
|
||||
if plotIt:
|
||||
from matplotlib.pylab import plt
|
||||
from scipy.interpolate import interp1d
|
||||
fig, ax = plt.subplots(1,1,figsize=(6,6))
|
||||
mesh.plotGrid(ax=ax, nodes=True, centers=True)
|
||||
ax.plot(xtopo,topo,'k',linewidth=1)
|
||||
# ax.plot(mesh.vectorNx, interp1d(xtopo,topo)(mesh.vectorNx),'--k',linewidth=3)
|
||||
ax.plot(mesh.vectorCCx, interp1d(xtopo,topo)(mesh.vectorCCx),'--k',linewidth=3)
|
||||
|
||||
|
||||
aveN2CC = Utils.sdiag(mesh.aveN2CC.T.sum(1))*mesh.aveN2CC.T
|
||||
a = aveN2CC * indcc
|
||||
a[a > 0] = 1.
|
||||
a[a < 0.25] = np.nan
|
||||
a = a.reshape(mesh.vnN, order='F')
|
||||
masked_array = np.ma.array(a, mask=np.isnan(a))
|
||||
ax.pcolor(mesh.vectorNx,mesh.vectorNy,masked_array.T, cmap = plt.cm.gray,alpha=0.2)
|
||||
plt.show()
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
run(plotIt=True)
|
||||
@@ -8,9 +8,9 @@ import EM_FDEM_Analytic_MagDipoleWholespace
|
||||
import EM_Schenkel_Morrison_Casing
|
||||
import EM_TDEM_1D_Inversion
|
||||
import FLOW_Richards_1D_Celia1990
|
||||
import Forward_BasicDirectCurrent
|
||||
import Inversion_IRLS
|
||||
import Inversion_Linear
|
||||
import Mesh_Basic_ForwardDC
|
||||
import Mesh_Basic_PlotImage
|
||||
import Mesh_Basic_Types
|
||||
import Mesh_Operators_CahnHilliard
|
||||
@@ -20,8 +20,9 @@ import Mesh_QuadTree_HangingNodes
|
||||
import Mesh_Tensor_Creation
|
||||
import MT_1D_ForwardAndInversion
|
||||
import MT_3D_Foward
|
||||
import Utils_surface2ind_topo
|
||||
|
||||
__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward"]
|
||||
__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_ForwardDC", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "Utils_surface2ind_topo"]
|
||||
|
||||
##### AUTOIMPORTS #####
|
||||
|
||||
|
||||
@@ -533,83 +533,6 @@ class ActiveCells(InjectActiveCells):
|
||||
FutureWarning)
|
||||
InjectActiveCells.__init__(self, mesh, indActive, valInactive, nC)
|
||||
|
||||
class InjectActiveCellsTopo(IdentityMap):
|
||||
"""
|
||||
Active model parameters. Extend for cells on topography to air cell (only works for tensor mesh)
|
||||
|
||||
"""
|
||||
|
||||
indActive = None #: Active Cells
|
||||
valInactive = None #: Values of inactive Cells
|
||||
nC = None #: Number of cells in the full model
|
||||
|
||||
def __init__(self, mesh, indActive, nC=None):
|
||||
self.mesh = mesh
|
||||
|
||||
self.nC = nC or mesh.nC
|
||||
|
||||
if indActive.dtype is not bool:
|
||||
z = np.zeros(self.nC,dtype=bool)
|
||||
z[indActive] = True
|
||||
indActive = z
|
||||
self.indActive = indActive
|
||||
|
||||
self.indInactive = np.logical_not(indActive)
|
||||
inds = np.nonzero(self.indActive)[0]
|
||||
self.P = sp.csr_matrix((np.ones(inds.size),(inds, range(inds.size))), shape=(self.nC, self.nP))
|
||||
|
||||
@property
|
||||
def shape(self):
|
||||
return (self.nC, self.nP)
|
||||
|
||||
@property
|
||||
def nP(self):
|
||||
"""Number of parameters in the model."""
|
||||
return self.indActive.sum()
|
||||
|
||||
def _transform(self, m):
|
||||
val_temp = np.zeros(self.mesh.nC)
|
||||
val_temp[self.indActive] = m
|
||||
valInactive = np.zeros(self.mesh.nC)
|
||||
#1D
|
||||
if self.mesh.dim == 1:
|
||||
z_temp = self.mesh.gridCC
|
||||
val_temp[~self.indActive] = val_temp[np.argmax(z_temp[self.indActive])]
|
||||
#2D
|
||||
elif self.mesh.dim == 2:
|
||||
act_temp = self.indActive.reshape((self.mesh.nCx, self.mesh.nCy), order = 'F')
|
||||
val_temp = val_temp.reshape((self.mesh.nCx, self.mesh.nCy), order = 'F')
|
||||
y_temp = self.mesh.gridCC[:,1].reshape((self.mesh.nCx, self.mesh.nCy), order = 'F')
|
||||
for i in range(self.mesh.nCx):
|
||||
act_tempx = act_temp[i,:] == 1
|
||||
val_temp[i,~act_tempx] = val_temp[i,np.argmax(y_temp[i,act_tempx])]
|
||||
valInactive[~self.indActive] = Utils.mkvc(val_temp)[~self.indActive]
|
||||
#3D
|
||||
elif self.mesh.dim == 3:
|
||||
act_temp = self.indActive.reshape((self.mesh.nCx*self.mesh.nCy, self.mesh.nCz), order = 'F')
|
||||
val_temp = val_temp.reshape((self.mesh.nCx*self.mesh.nCy, self.mesh.nCz), order = 'F')
|
||||
z_temp = self.mesh.gridCC[:,2].reshape((self.mesh.nCx*self.mesh.nCy, self.mesh.nCz), order = 'F')
|
||||
for i in range(self.mesh.nCx*self.mesh.nCy):
|
||||
act_tempxy = act_temp[i,:] == 1
|
||||
val_temp[i,~act_tempxy] = val_temp[i,np.argmax(z_temp[i,act_tempxy])]
|
||||
valInactive[~self.indActive] = Utils.mkvc(val_temp)[~self.indActive]
|
||||
|
||||
self.valInactive = valInactive
|
||||
|
||||
return self.P*m + self.valInactive
|
||||
|
||||
def inverse(self, D):
|
||||
return self.P.T*D
|
||||
|
||||
def deriv(self, m):
|
||||
return self.P
|
||||
|
||||
class ActiveCellsTopo(InjectActiveCellsTopo):
|
||||
def __init__(self, mesh, indActive, valInactive, nC=None):
|
||||
warnings.warn(
|
||||
"`ActiveCellsTopo` is deprecated and will be removed in future versions. Use `InjectActiveCellsTopo` instead",
|
||||
FutureWarning)
|
||||
InjectActiveCellsTopo.__init__(self, mesh, indActive, valInactive, nC)
|
||||
|
||||
class Weighting(IdentityMap):
|
||||
"""
|
||||
|
||||
@@ -2,6 +2,7 @@ from SimPEG import Utils, np
|
||||
from BaseMesh import BaseRectangularMesh
|
||||
from DiffOperators import DiffOperators
|
||||
from InnerProducts import InnerProducts
|
||||
from View import CurvView
|
||||
|
||||
# Some helper functions.
|
||||
length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5
|
||||
@@ -10,7 +11,7 @@ normalize2D = lambda x: x/np.kron(np.ones((1, 2)), Utils.mkvc(length2D(x), 2))
|
||||
normalize3D = lambda x: x/np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2))
|
||||
|
||||
|
||||
class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
|
||||
class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvView):
|
||||
"""
|
||||
CurvilinearMesh is a mesh class that deals with curvilinear meshes.
|
||||
|
||||
@@ -330,102 +331,6 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
|
||||
|
||||
|
||||
|
||||
#############################################
|
||||
# Plotting Functions #
|
||||
#############################################
|
||||
|
||||
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.
|
||||
|
||||
|
||||
.. plot::
|
||||
:include-source:
|
||||
|
||||
from SimPEG import Mesh, Utils
|
||||
X, Y = Utils.exampleLrmGrid([3,3],'rotate')
|
||||
M = Mesh.CurvilinearMesh([X, Y])
|
||||
M.plotGrid(showIt=True)
|
||||
|
||||
"""
|
||||
import matplotlib.pyplot as plt
|
||||
import matplotlib
|
||||
from mpl_toolkits.mplot3d import Axes3D
|
||||
mkvc = Utils.mkvc
|
||||
|
||||
axOpts = {'projection':'3d'} if self.dim == 3 else {}
|
||||
if ax is None: ax = plt.subplot(111, **axOpts)
|
||||
|
||||
NN = self.r(self.gridN, 'N', 'N', 'M')
|
||||
if self.dim == 2:
|
||||
|
||||
if lines:
|
||||
X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten()
|
||||
Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten()
|
||||
|
||||
X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten()
|
||||
Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten()
|
||||
|
||||
X = np.r_[X1, X2]
|
||||
Y = np.r_[Y1, Y2]
|
||||
|
||||
ax.plot(X, Y, 'b-')
|
||||
if centers:
|
||||
ax.plot(self.gridCC[:,0],self.gridCC[:,1],'ro')
|
||||
|
||||
# Nx = self.r(self.normals, 'F', 'Fx', 'V')
|
||||
# Ny = self.r(self.normals, 'F', 'Fy', 'V')
|
||||
# Tx = self.r(self.tangents, 'E', 'Ex', 'V')
|
||||
# Ty = self.r(self.tangents, 'E', 'Ey', 'V')
|
||||
|
||||
# ax.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo')
|
||||
|
||||
# nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten()
|
||||
# nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten()
|
||||
# ax.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs')
|
||||
# ax.plot(nX, nY, 'r-')
|
||||
|
||||
# nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten()
|
||||
# nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten()
|
||||
# #ax.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs')
|
||||
# ax.plot(nX, nY, 'g-')
|
||||
|
||||
# tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten()
|
||||
# tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten()
|
||||
# ax.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^')
|
||||
# ax.plot(tX, tY, 'r-')
|
||||
|
||||
# nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten()
|
||||
# nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten()
|
||||
# #ax.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^')
|
||||
# ax.plot(nX, nY, 'g-')
|
||||
|
||||
elif self.dim == 3:
|
||||
X1 = np.c_[mkvc(NN[0][:-1, :, :]), mkvc(NN[0][1:, :, :]), mkvc(NN[0][:-1, :, :])*np.nan].flatten()
|
||||
Y1 = np.c_[mkvc(NN[1][:-1, :, :]), mkvc(NN[1][1:, :, :]), mkvc(NN[1][:-1, :, :])*np.nan].flatten()
|
||||
Z1 = np.c_[mkvc(NN[2][:-1, :, :]), mkvc(NN[2][1:, :, :]), mkvc(NN[2][:-1, :, :])*np.nan].flatten()
|
||||
|
||||
X2 = np.c_[mkvc(NN[0][:, :-1, :]), mkvc(NN[0][:, 1:, :]), mkvc(NN[0][:, :-1, :])*np.nan].flatten()
|
||||
Y2 = np.c_[mkvc(NN[1][:, :-1, :]), mkvc(NN[1][:, 1:, :]), mkvc(NN[1][:, :-1, :])*np.nan].flatten()
|
||||
Z2 = np.c_[mkvc(NN[2][:, :-1, :]), mkvc(NN[2][:, 1:, :]), mkvc(NN[2][:, :-1, :])*np.nan].flatten()
|
||||
|
||||
X3 = np.c_[mkvc(NN[0][:, :, :-1]), mkvc(NN[0][:, :, 1:]), mkvc(NN[0][:, :, :-1])*np.nan].flatten()
|
||||
Y3 = np.c_[mkvc(NN[1][:, :, :-1]), mkvc(NN[1][:, :, 1:]), mkvc(NN[1][:, :, :-1])*np.nan].flatten()
|
||||
Z3 = np.c_[mkvc(NN[2][:, :, :-1]), mkvc(NN[2][:, :, 1:]), mkvc(NN[2][:, :, :-1])*np.nan].flatten()
|
||||
|
||||
X = np.r_[X1, X2, X3]
|
||||
Y = np.r_[Y1, Y2, Y3]
|
||||
Z = np.r_[Z1, Z2, Z3]
|
||||
|
||||
ax.plot(X, Y, 'b', zs=Z)
|
||||
ax.set_zlabel('x3')
|
||||
|
||||
ax.grid(True)
|
||||
ax.set_xlabel('x1')
|
||||
ax.set_ylabel('x2')
|
||||
|
||||
if showIt: plt.show()
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
nc = 5
|
||||
h1 = np.cumsum(np.r_[0, np.ones(nc)/(nc)])
|
||||
|
||||
+78
-40
@@ -552,7 +552,8 @@ class CurvView(object):
|
||||
def __init__(self):
|
||||
pass
|
||||
|
||||
def plotGrid(self, length=0.05, showIt=False):
|
||||
|
||||
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.
|
||||
|
||||
|
||||
@@ -560,60 +561,63 @@ class CurvView(object):
|
||||
:include-source:
|
||||
|
||||
from SimPEG import Mesh, Utils
|
||||
X, Y = Utils.exampleCurvGird([3,3],'rotate')
|
||||
X, Y = Utils.exampleLrmGrid([3,3],'rotate')
|
||||
M = Mesh.CurvilinearMesh([X, Y])
|
||||
M.plotGrid(showIt=True)
|
||||
|
||||
"""
|
||||
import matplotlib.pyplot as plt
|
||||
import matplotlib
|
||||
from mpl_toolkits.mplot3d import Axes3D
|
||||
|
||||
axOpts = {'projection':'3d'} if self.dim == 3 else {}
|
||||
if ax is None: ax = plt.subplot(111, **axOpts)
|
||||
|
||||
NN = self.r(self.gridN, 'N', 'N', 'M')
|
||||
if self.dim == 2:
|
||||
fig = plt.figure(2)
|
||||
fig.clf()
|
||||
ax = plt.subplot(111)
|
||||
X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten()
|
||||
Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten()
|
||||
|
||||
X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten()
|
||||
Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten()
|
||||
if lines:
|
||||
X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten()
|
||||
Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten()
|
||||
|
||||
X = np.r_[X1, X2]
|
||||
Y = np.r_[Y1, Y2]
|
||||
X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten()
|
||||
Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten()
|
||||
|
||||
plt.plot(X, Y)
|
||||
X = np.r_[X1, X2]
|
||||
Y = np.r_[Y1, Y2]
|
||||
|
||||
plt.hold(True)
|
||||
Nx = self.r(self.normals, 'F', 'Fx', 'V')
|
||||
Ny = self.r(self.normals, 'F', 'Fy', 'V')
|
||||
Tx = self.r(self.tangents, 'E', 'Ex', 'V')
|
||||
Ty = self.r(self.tangents, 'E', 'Ey', 'V')
|
||||
ax.plot(X, Y, 'b-')
|
||||
if centers:
|
||||
ax.plot(self.gridCC[:,0],self.gridCC[:,1],'ro')
|
||||
|
||||
plt.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo')
|
||||
# Nx = self.r(self.normals, 'F', 'Fx', 'V')
|
||||
# Ny = self.r(self.normals, 'F', 'Fy', 'V')
|
||||
# Tx = self.r(self.tangents, 'E', 'Ex', 'V')
|
||||
# Ty = self.r(self.tangents, 'E', 'Ey', 'V')
|
||||
|
||||
nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten()
|
||||
nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten()
|
||||
plt.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs')
|
||||
plt.plot(nX, nY, 'r-')
|
||||
# ax.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo')
|
||||
|
||||
nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten()
|
||||
nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten()
|
||||
#plt.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs')
|
||||
plt.plot(nX, nY, 'g-')
|
||||
# nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten()
|
||||
# nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten()
|
||||
# ax.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs')
|
||||
# ax.plot(nX, nY, 'r-')
|
||||
|
||||
tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten()
|
||||
tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten()
|
||||
plt.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^')
|
||||
plt.plot(tX, tY, 'r-')
|
||||
# nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten()
|
||||
# nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten()
|
||||
# #ax.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs')
|
||||
# ax.plot(nX, nY, 'g-')
|
||||
|
||||
nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten()
|
||||
nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten()
|
||||
#plt.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^')
|
||||
plt.plot(nX, nY, 'g-')
|
||||
plt.axis('equal')
|
||||
# tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten()
|
||||
# tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten()
|
||||
# ax.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^')
|
||||
# ax.plot(tX, tY, 'r-')
|
||||
|
||||
# nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten()
|
||||
# nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten()
|
||||
# #ax.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^')
|
||||
# ax.plot(nX, nY, 'g-')
|
||||
|
||||
elif self.dim == 3:
|
||||
fig = plt.figure(3)
|
||||
fig.clf()
|
||||
ax = fig.add_subplot(111, projection='3d')
|
||||
X1 = np.c_[mkvc(NN[0][:-1, :, :]), mkvc(NN[0][1:, :, :]), mkvc(NN[0][:-1, :, :])*np.nan].flatten()
|
||||
Y1 = np.c_[mkvc(NN[1][:-1, :, :]), mkvc(NN[1][1:, :, :]), mkvc(NN[1][:-1, :, :])*np.nan].flatten()
|
||||
Z1 = np.c_[mkvc(NN[2][:-1, :, :]), mkvc(NN[2][1:, :, :]), mkvc(NN[2][:-1, :, :])*np.nan].flatten()
|
||||
@@ -630,16 +634,50 @@ class CurvView(object):
|
||||
Y = np.r_[Y1, Y2, Y3]
|
||||
Z = np.r_[Z1, Z2, Z3]
|
||||
|
||||
plt.plot(X, Y, 'b', zs=Z)
|
||||
ax.plot(X, Y, 'b', zs=Z)
|
||||
ax.set_zlabel('x3')
|
||||
|
||||
ax.grid(True)
|
||||
ax.hold(False)
|
||||
ax.set_xlabel('x1')
|
||||
ax.set_ylabel('x2')
|
||||
|
||||
if showIt: plt.show()
|
||||
|
||||
def plotImage(self, I, ax=None, showIt=False, grid=False, clim=None):
|
||||
if self.dim == 3: raise NotImplementedError('This is not yet done!')
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
import matplotlib
|
||||
from mpl_toolkits.mplot3d import Axes3D
|
||||
import matplotlib.colors as colors
|
||||
import matplotlib.cm as cmx
|
||||
|
||||
if ax is None: ax = plt.subplot(111)
|
||||
jet = cm = plt.get_cmap('jet')
|
||||
cNorm = colors.Normalize(
|
||||
vmin=I.min() if clim is None else clim[0],
|
||||
vmax=I.max() if clim is None else clim[1])
|
||||
|
||||
scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=jet)
|
||||
# ax.set_xlim((self.x0[0], self.h[0].sum()))
|
||||
# ax.set_ylim((self.x0[1], self.h[1].sum()))
|
||||
|
||||
Nx = self.r(self.gridN[:,0],'N','N','M')
|
||||
Ny = self.r(self.gridN[:,1],'N','N','M')
|
||||
cell = self.r(I,'CC','CC','M')
|
||||
|
||||
for ii in range(self.nCx):
|
||||
for jj in range(self.nCy):
|
||||
I = [ii,ii+1,ii+1,ii]
|
||||
J = [jj,jj,jj+1,jj+1]
|
||||
ax.add_patch(plt.Polygon(np.c_[Nx[I,J],Ny[I,J]], facecolor=scalarMap.to_rgba(cell[ii,jj]), edgecolor='k' if grid else 'none'))
|
||||
|
||||
scalarMap._A = [] # http://stackoverflow.com/questions/8342549/matplotlib-add-colorbar-to-a-sequence-of-line-plots
|
||||
ax.set_xlabel('x')
|
||||
ax.set_ylabel('y')
|
||||
if showIt: plt.show()
|
||||
return [scalarMap]
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
from SimPEG import *
|
||||
|
||||
@@ -1008,4 +1008,4 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
|
||||
indx = ((self.xc<=self.lower) & (delx < 0)) | ((self.xc>=self.upper) & (delx > 0))
|
||||
delx[indx] = 0.
|
||||
|
||||
return delx
|
||||
return delx
|
||||
+442
-184
@@ -1,4 +1,6 @@
|
||||
import Utils, Maps, Mesh, numpy as np, scipy.sparse as sp
|
||||
import Utils, Maps, Mesh
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
|
||||
class RegularizationMesh(object):
|
||||
"""
|
||||
@@ -403,7 +405,238 @@ class BaseRegularization(object):
|
||||
|
||||
return mD.T * ( self.W.T * ( self.W * ( mD * v) ) )
|
||||
|
||||
class Tikhonov(BaseRegularization):
|
||||
class Simple(BaseRegularization):
|
||||
"""
|
||||
Simple regularization that does not include length scales in the derivatives.
|
||||
"""
|
||||
|
||||
mrefInSmooth = False #: include mref in the smoothness?
|
||||
alpha_s = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Wsmall'], "Smallness weight")
|
||||
alpha_x = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
|
||||
alpha_y = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
|
||||
alpha_z = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction")
|
||||
cell_weights = 1.
|
||||
|
||||
def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
|
||||
BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
|
||||
|
||||
if isinstance(self.cell_weights,float):
|
||||
self.cell_weights = np.ones(self.regmesh.nC) * self.cell_weights
|
||||
|
||||
@property
|
||||
def Wsmall(self):
|
||||
"""Regularization matrix Wsmall"""
|
||||
if getattr(self,'_Wsmall', None) is None:
|
||||
self._Wsmall = Utils.sdiag((self.alpha_s*self.cell_weights)**0.5)
|
||||
return self._Wsmall
|
||||
|
||||
@property
|
||||
def Wx(self):
|
||||
"""Regularization matrix Wx"""
|
||||
if getattr(self, '_Wx', None) is None:
|
||||
self._Wx = Utils.sdiag((self.alpha_x * (self.regmesh.aveCC2Fx*self.cell_weights))**0.5)*self.regmesh.cellDiffxStencil
|
||||
return self._Wx
|
||||
|
||||
@property
|
||||
def Wy(self):
|
||||
"""Regularization matrix Wy"""
|
||||
if getattr(self, '_Wy', None) is None:
|
||||
self._Wy = Utils.sdiag((self.alpha_y * (self.regmesh.aveCC2Fy*self.cell_weights))**0.5)*self.regmesh.cellDiffyStencil
|
||||
return self._Wy
|
||||
|
||||
@property
|
||||
def Wz(self):
|
||||
"""Regularization matrix Wz"""
|
||||
if getattr(self, '_Wz', None) is None:
|
||||
self._Wz = Utils.sdiag((self.alpha_z * (self.regmesh.aveCC2Fz*self.cell_weights))**0.5)*self.regmesh.cellDiffzStencil
|
||||
return self._Wz
|
||||
|
||||
# @property
|
||||
# def Wsmooth(self):
|
||||
# """Full smoothness regularization matrix W"""
|
||||
# print 'wtf why are we using Wsmooth'
|
||||
# raise NotImplementedError
|
||||
# if getattr(self, '_Wsmooth', None) is None:
|
||||
# wlist = (self.Wx,)
|
||||
# if self.regmesh.dim > 1:
|
||||
# wlist += (self.Wy,)
|
||||
# if self.regmesh.dim > 2:
|
||||
# wlist += (self.Wz,)
|
||||
# self._Wsmooth = sp.vstack(wlist)
|
||||
# return self._Wsmooth
|
||||
#
|
||||
# @property
|
||||
# def W(self):
|
||||
# """Full regularization matrix W"""
|
||||
# print 'wtf why are we using W'
|
||||
# if getattr(self, '_W', None) is None:
|
||||
# wlist = (self.Wsmall, self.Wx)
|
||||
# if self.regmesh.dim > 1:
|
||||
# wlist += (self.Wy,)
|
||||
# if self.regmesh.dim > 2:
|
||||
# wlist += (self.Wz,)
|
||||
# self._W = sp.vstack(wlist)
|
||||
# return self._W
|
||||
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmall(self, m):
|
||||
r = self.Wsmall * ( self.mapping * (m - self.mref) )
|
||||
return 0.5 * r.dot(r)
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmallDeriv(self, m):
|
||||
r = self.Wsmall * ( self.mapping * (m - self.mref) )
|
||||
return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) )
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmall2Deriv(self, m, v = None):
|
||||
rDeriv = self.Wsmall * ( self.mapping.deriv(m - self.mref) )
|
||||
if v is not None:
|
||||
return rDeriv.T * (rDeriv * v)
|
||||
return rDeriv.T * rDeriv
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmoothx(self, m):
|
||||
if self.mrefInSmooth == True:
|
||||
r = self.Wx * ( self.mapping * (m - self.mref) )
|
||||
elif self.mrefInSmooth == False:
|
||||
r = self.Wx * ( self.mapping * (m) )
|
||||
return 0.5 * r.dot(r)
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmoothy(self, m):
|
||||
if self.mrefInSmooth == True:
|
||||
r = self.Wy * ( self.mapping * (m - self.mref) )
|
||||
elif self.mrefInSmooth == False:
|
||||
r = self.Wy * ( self.mapping * (m) )
|
||||
return 0.5 * r.dot(r)
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmoothz(self, m):
|
||||
if self.mrefInSmooth == True:
|
||||
r = self.Wz * ( self.mapping * (m - self.mref) )
|
||||
elif self.mrefInSmooth == False:
|
||||
r = self.Wz * ( self.mapping * (m) )
|
||||
return 0.5 * r.dot(r)
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmooth(self, m):
|
||||
phiSmooth = self._evalSmoothx(m)
|
||||
if self.regmesh.dim > 1:
|
||||
phiSmooth += self._evalSmoothy(m)
|
||||
if self.regmesh.dim > 2:
|
||||
phiSmooth += self._evalSmoothz(m)
|
||||
return phiSmooth
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmoothxDeriv(self, m):
|
||||
if self.mrefInSmooth == True:
|
||||
r = self.Wx * ( self.mapping * ( m - self.mref ) )
|
||||
return r.T * ( self.Wx * self.mapping.deriv(m - self.mref) )
|
||||
elif self.mrefInSmooth == False:
|
||||
r = self.Wx * ( self.mapping * m )
|
||||
return r.T * ( self.Wx * self.mapping.deriv(m) )
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmoothx2Deriv(self, m, v=None):
|
||||
if self.mrefInSmooth == True:
|
||||
rDeriv = self.Wx * ( self.mapping.deriv( m - self.mref ) )
|
||||
elif self.mrefInSmooth == False:
|
||||
rDeriv = self.Wx * ( self.mapping.deriv(m) )
|
||||
|
||||
if v is not None:
|
||||
return rDeriv.T * ( rDeriv * v )
|
||||
return rDeriv.T * rDeriv
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmoothyDeriv(self, m):
|
||||
if self.mrefInSmooth == True:
|
||||
r = self.Wy * ( self.mapping * ( m - self.mref ) )
|
||||
return r.T * ( self.Wy * self.mapping.deriv(m - self.mref) )
|
||||
elif self.mrefInSmooth == False:
|
||||
r = self.Wy * ( self.mapping * m )
|
||||
return r.T * ( self.Wy * self.mapping.deriv(m) )
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmoothy2Deriv(self, m, v=None):
|
||||
if self.mrefInSmooth == True:
|
||||
rDeriv = self.Wy * ( self.mapping.deriv( m - self.mref ) )
|
||||
elif self.mrefInSmooth == False:
|
||||
rDeriv = self.Wy * ( self.mapping.deriv(m) )
|
||||
|
||||
if v is not None:
|
||||
return rDeriv.T * ( rDeriv * v )
|
||||
return rDeriv.T * rDeriv
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmoothzDeriv(self, m):
|
||||
if self.mrefInSmooth == True:
|
||||
r = self.Wz * ( self.mapping * ( m - self.mref ) )
|
||||
return r.T * ( self.Wz * self.mapping.deriv(m - self.mref) )
|
||||
elif self.mrefInSmooth == False:
|
||||
r = self.Wz * ( self.mapping * m )
|
||||
return r.T * ( self.Wz * self.mapping.deriv(m) )
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmoothz2Deriv(self, m, v=None):
|
||||
if self.mrefInSmooth == True:
|
||||
rDeriv = self.Wz * ( self.mapping.deriv( m - self.mref ) )
|
||||
elif self.mrefInSmooth == False:
|
||||
rDeriv = self.Wz * ( self.mapping.deriv(m) )
|
||||
|
||||
if v is not None:
|
||||
return rDeriv.T * ( rDeriv * v )
|
||||
return rDeriv.T * rDeriv
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmoothDeriv(self, m):
|
||||
deriv = self._evalSmoothxDeriv(m)
|
||||
if self.regmesh.dim > 1:
|
||||
deriv += self._evalSmoothyDeriv(m)
|
||||
if self.regmesh.dim > 2:
|
||||
deriv += self._evalSmoothzDeriv(m)
|
||||
return deriv
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmooth2Deriv(self, m, v=None):
|
||||
deriv = self._evalSmoothx2Deriv(m, v)
|
||||
if self.regmesh.dim > 1:
|
||||
deriv += self._evalSmoothy2Deriv(m, v)
|
||||
if self.regmesh.dim > 2:
|
||||
deriv += self._evalSmoothz2Deriv(m, v)
|
||||
return deriv
|
||||
|
||||
|
||||
@Utils.timeIt
|
||||
def eval(self, m):
|
||||
return self._evalSmall(m) + self._evalSmooth(m)
|
||||
|
||||
@Utils.timeIt
|
||||
def evalDeriv(self, m):
|
||||
"""
|
||||
The regularization is:
|
||||
|
||||
.. math::
|
||||
|
||||
R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
|
||||
|
||||
So the derivative is straight forward:
|
||||
|
||||
.. math::
|
||||
|
||||
R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
|
||||
|
||||
"""
|
||||
return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m)
|
||||
|
||||
@Utils.timeIt
|
||||
def eval2Deriv(self, m, v=None):
|
||||
return self._evalSmall2Deriv(m, v) + self._evalSmooth2Deriv(m, v)
|
||||
|
||||
|
||||
|
||||
class Tikhonov(Simple):
|
||||
"""
|
||||
L2 Tikhonov regularization with both smallness and smoothness (first order
|
||||
derivative) contributions.
|
||||
@@ -493,56 +726,131 @@ class Tikhonov(BaseRegularization):
|
||||
self._Wzz = Utils.sdiag((self.regmesh.vol*self.alpha_zz)**0.5)*self.regmesh.faceDiffz*self.regmesh.cellDiffz
|
||||
return self._Wzz
|
||||
|
||||
|
||||
@property
|
||||
def Wsmooth(self):
|
||||
def Wsmooth2(self):
|
||||
"""Full smoothness regularization matrix W"""
|
||||
if getattr(self, '_Wsmooth', None) is None:
|
||||
wlist = (self.Wx, self.Wxx)
|
||||
wlist = (self.Wxx)
|
||||
if self.regmesh.dim > 1:
|
||||
wlist += (self.Wy, self.Wyy)
|
||||
wlist += (self.Wyy)
|
||||
if self.regmesh.dim > 2:
|
||||
wlist += (self.Wz, self.Wzz)
|
||||
wlist += (self.Wzz)
|
||||
self._Wsmooth = sp.vstack(wlist)
|
||||
return self._Wsmooth
|
||||
|
||||
@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 self._W
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmall(self, m):
|
||||
r = self.Wsmall * ( self.mapping * (m - self.mref) )
|
||||
return 0.5 * r.dot(r)
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmooth(self, m):
|
||||
def _evalSmoothxx(self, m):
|
||||
if self.mrefInSmooth == True:
|
||||
r = self.Wsmooth * ( self.mapping * (m - self.mref) )
|
||||
r = self.Wxx * ( self.mapping * (m - self.mref) )
|
||||
elif self.mrefInSmooth == False:
|
||||
r = self.Wsmooth * ( self.mapping * (m) )
|
||||
r = self.Wxx * ( self.mapping * (m) )
|
||||
return 0.5 * r.dot(r)
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmoothyy(self, m):
|
||||
if self.mrefInSmooth == True:
|
||||
r = self.Wyy * ( self.mapping * (m - self.mref) )
|
||||
elif self.mrefInSmooth == False:
|
||||
r = self.Wyy * ( self.mapping * (m) )
|
||||
return 0.5 * r.dot(r)
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmoothzz(self, m):
|
||||
if self.mrefInSmooth == True:
|
||||
r = self.Wzz * ( self.mapping * (m - self.mref) )
|
||||
elif self.mrefInSmooth == False:
|
||||
r = self.Wzz * ( self.mapping * (m) )
|
||||
return 0.5 * r.dot(r)
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmooth2(self, m):
|
||||
phiSmooth2 = self._evalSmoothxx(m)
|
||||
if self.regmesh.dim > 1:
|
||||
phiSmooth2 += self._evalSmoothyy(m)
|
||||
if self.regmesh.dim > 2:
|
||||
phiSmooth2 += self._evalSmoothzz(m)
|
||||
return phiSmooth2
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmoothxxDeriv(self, m):
|
||||
if self.mrefInSmooth == True:
|
||||
r = self.Wxx * ( self.mapping * ( m - self.mref ) )
|
||||
return r.T * ( self.Wxx * self.mapping.deriv(m - self.mref) )
|
||||
elif self.mrefInSmooth == False:
|
||||
r = self.Wxx * ( self.mapping * m )
|
||||
return r.T * ( self.Wxx * self.mapping.deriv(m) )
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmoothyyDeriv(self, m):
|
||||
if self.mrefInSmooth == True:
|
||||
r = self.Wyy * ( self.mapping * ( m - self.mref ) )
|
||||
return r.T * ( self.Wyy * self.mapping.deriv(m - self.mref) )
|
||||
elif self.mrefInSmooth == False:
|
||||
r = self.Wyy * ( self.mapping * m )
|
||||
return r.T * ( self.Wyy * self.mapping.deriv(m) )
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmoothzzDeriv(self, m):
|
||||
if self.mrefInSmooth == True:
|
||||
r = self.Wzz * ( self.mapping * ( m - self.mref ) )
|
||||
return r.T * ( self.Wzz * self.mapping.deriv(m - self.mref) )
|
||||
elif self.mrefInSmooth == False:
|
||||
r = self.Wzz * ( self.mapping * m )
|
||||
return r.T * ( self.Wzz * self.mapping.deriv(m) )
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmoothxx2Deriv(self, m, v=None):
|
||||
if self.mrefInSmooth == True:
|
||||
rDeriv = self.Wxx * ( self.mapping.deriv( m - self.mref ) )
|
||||
elif self.mrefInSmooth == False:
|
||||
rDeriv = self.Wxx * self.mapping.deriv(m)
|
||||
if v is not None:
|
||||
return rDeriv.T * (rDeriv * v)
|
||||
return rDeriv.T * rDeriv
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmoothyy2Deriv(self, m, v=None):
|
||||
if self.mrefInSmooth == True:
|
||||
rDeriv = self.Wyy * ( self.mapping.deriv( m - self.mref ) )
|
||||
elif self.mrefInSmooth == False:
|
||||
rDeriv = self.Wyy * self.mapping.deriv(m)
|
||||
if v is not None:
|
||||
return rDeriv.T * (rDeriv * v)
|
||||
return rDeriv.T * rDeriv
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmoothzz2Deriv(self, m, v=None):
|
||||
if self.mrefInSmooth == True:
|
||||
rDeriv = self.Wzz * ( self.mapping.deriv( m - self.mref ) )
|
||||
elif self.mrefInSmooth == False:
|
||||
rDeriv = self.Wzz * self.mapping.deriv(m)
|
||||
if v is not None:
|
||||
return rDeriv.T * (rDeriv * v)
|
||||
return rDeriv.T * rDeriv
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmoothDeriv2(self, m):
|
||||
deriv = self._evalSmoothxxDeriv(m)
|
||||
if self.regmesh.dim > 1:
|
||||
deriv += self._evalSmoothyyDeriv(m)
|
||||
if self.regmesh.dim > 2:
|
||||
deriv += self._evalSmoothzzDeriv(m)
|
||||
return deriv
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmooth2Deriv2(self, m, v=None):
|
||||
deriv = self._evalSmoothxx2Deriv(m, v)
|
||||
if self.regmesh.dim > 1:
|
||||
deriv += self._evalSmoothyy2Deriv(m, v)
|
||||
if self.regmesh.dim > 2:
|
||||
deriv += self._evalSmoothzz2Deriv(m, v)
|
||||
return deriv
|
||||
|
||||
|
||||
@Utils.timeIt
|
||||
def eval(self, m):
|
||||
return self._evalSmall(m) + self._evalSmooth(m)
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmallDeriv(self,m):
|
||||
r = self.Wsmall * ( self.mapping * (m - self.mref) )
|
||||
return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) )
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmoothDeriv(self,m):
|
||||
if self.mrefInSmooth == True:
|
||||
r = self.Wsmooth * ( self.mapping * ( m - self.mref ) )
|
||||
return r.T * ( self.Wsmooth * self.mapping.deriv(m - self.mref) )
|
||||
elif self.mrefInSmooth == False:
|
||||
r = self.Wsmooth * ( self.mapping * m )
|
||||
return r.T * ( self.Wsmooth * self.mapping.deriv(m) )
|
||||
return self._evalSmall(m) + self._evalSmooth(m) + self._evalSmooth2(m)
|
||||
|
||||
@Utils.timeIt
|
||||
def evalDeriv(self, m):
|
||||
@@ -560,184 +868,134 @@ class Tikhonov(BaseRegularization):
|
||||
R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
|
||||
|
||||
"""
|
||||
return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m)
|
||||
return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + self._evalSmoothDeriv2(m)
|
||||
|
||||
def eval2Deriv(self, m, v=None):
|
||||
"""
|
||||
The regularization is:
|
||||
|
||||
.. math::
|
||||
|
||||
R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
|
||||
|
||||
So the derivative is straight forward:
|
||||
|
||||
.. math::
|
||||
|
||||
R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
|
||||
|
||||
"""
|
||||
return self._evalSmall2Deriv(m, v) + self._evalSmooth2Deriv(m, v) + self._evalSmooth2Deriv2(m, v)
|
||||
|
||||
|
||||
class Simple(Tikhonov):
|
||||
|
||||
class Sparse(Simple):
|
||||
"""
|
||||
Simple regularization that does not include length scales in the derivatives.
|
||||
The regularization is:
|
||||
|
||||
.. math::
|
||||
|
||||
R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top R^\\top R W(m-m_\\text{ref})}
|
||||
|
||||
where the IRLS weight
|
||||
|
||||
.. math::
|
||||
|
||||
R = \eta TO FINISH LATER!!!
|
||||
|
||||
So the derivative is straight forward:
|
||||
|
||||
.. math::
|
||||
|
||||
R(m) = \mathbf{W^\\top R^\\top R W (m-m_\\text{ref})}
|
||||
|
||||
The IRLS weights are recomputed after each beta solves.
|
||||
It is strongly recommended to do a few Gauss-Newton iterations
|
||||
before updating.
|
||||
"""
|
||||
|
||||
mrefInSmooth = False #: SMOOTH and SMOOTH_MOD_DIF options
|
||||
alpha_s = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Wsmall'], "Smallness weight")
|
||||
alpha_x = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
|
||||
alpha_y = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
|
||||
alpha_z = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction")
|
||||
wght = 1.
|
||||
|
||||
# set default values
|
||||
eps_p = 1e-1 # Threshold value for the model norm
|
||||
eps_q = 1e-1 # Threshold value for the model gradient norm
|
||||
curModel = None # Requires model to compute the weights
|
||||
l2model = None
|
||||
gamma = 1. # Model norm scaling to smooth out convergence
|
||||
norms = [0., 2., 2., 2.] # Values for norm on (m, dmdx, dmdy, dmdz)
|
||||
cell_weights = 1. # Consider overwriting with sensitivity weights
|
||||
|
||||
def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
|
||||
BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
|
||||
Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
|
||||
|
||||
if isinstance(self.wght,float):
|
||||
self.wght = np.ones(self.regmesh.nC) * self.wght
|
||||
if isinstance(self.cell_weights,float):
|
||||
self.cell_weights = np.ones(self.regmesh.nC) * self.cell_weights
|
||||
|
||||
@property
|
||||
def Wsmall(self):
|
||||
"""Regularization matrix Wsmall"""
|
||||
if getattr(self,'_Wsmall', None) is None:
|
||||
self._Wsmall = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.wght)**0.5)
|
||||
if getattr(self, 'curModel', None) is None:
|
||||
self.Rs = Utils.speye(self.regmesh.nC)
|
||||
|
||||
else:
|
||||
f_m = self.mapping * (self.curModel - self.reg.mref)
|
||||
self.rs = self.R(f_m , self.eps_p, self.norms[0])
|
||||
self.Rs = Utils.sdiag( self.rs )
|
||||
|
||||
self._Wsmall = Utils.sdiag((self.alpha_s*self.gamma*self.cell_weights)**0.5)*self.Rs
|
||||
|
||||
return self._Wsmall
|
||||
|
||||
@property
|
||||
def Wx(self):
|
||||
"""Regularization matrix Wx"""
|
||||
if getattr(self, '_Wx', None) is None:
|
||||
self._Wx = Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.regmesh.cellDiffxStencil
|
||||
if getattr(self,'_Wx', None) is None:
|
||||
if getattr(self, 'curModel', None) is None:
|
||||
self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
|
||||
|
||||
else:
|
||||
f_m = self.regmesh.cellDiffxStencil * (self.mapping * self.curModel)
|
||||
self.rx = self.R( f_m , self.eps_q, self.norms[1])
|
||||
self.Rx = Utils.sdiag( self.rx )
|
||||
|
||||
self._Wx = Utils.sdiag(( self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.cell_weights))**0.5)*self.Rx*self.regmesh.cellDiffxStencil
|
||||
|
||||
return self._Wx
|
||||
|
||||
@property
|
||||
def Wy(self):
|
||||
"""Regularization matrix Wy"""
|
||||
if getattr(self, '_Wy', None) is None:
|
||||
self._Wy = Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol * self.alpha_y*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.regmesh.cellDiffyStencil
|
||||
if getattr(self,'_Wy', None) is None:
|
||||
if getattr(self, 'curModel', None) is None:
|
||||
self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0])
|
||||
|
||||
else:
|
||||
f_m = self.regmesh.cellDiffyStencil * (self.mapping * self.curModel)
|
||||
self.ry = self.R( f_m , self.eps_q, self.norms[2])
|
||||
self.Ry = Utils.sdiag( self.ry )
|
||||
|
||||
self._Wy = Utils.sdiag((self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.cell_weights))**0.5)*self.Ry*self.regmesh.cellDiffyStencil
|
||||
|
||||
return self._Wy
|
||||
|
||||
@property
|
||||
def Wz(self):
|
||||
"""Regularization matrix Wz"""
|
||||
if getattr(self, '_Wz', None) is None:
|
||||
self._Wz = Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.regmesh.cellDiffzStencil
|
||||
if getattr(self,'_Wz', None) is None:
|
||||
if getattr(self, 'curModel', None) is None:
|
||||
self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0])
|
||||
|
||||
else:
|
||||
f_m = self.regmesh.cellDiffzStencil * (self.mapping * self.curModel)
|
||||
self.rz = self.R( f_m , self.eps_q, self.norms[3])
|
||||
self.Rz = Utils.sdiag( self.rz )
|
||||
|
||||
self._Wz = Utils.sdiag((self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.cell_weights))**0.5)*self.Rz*self.regmesh.cellDiffzStencil
|
||||
|
||||
return self._Wz
|
||||
|
||||
@property
|
||||
def Wsmooth(self):
|
||||
"""Full smoothness regularization matrix W"""
|
||||
if getattr(self, '_Wsmooth', None) is None:
|
||||
wlist = (self.Wx,)
|
||||
if self.regmesh.dim > 1:
|
||||
wlist += (self.Wy,)
|
||||
if self.regmesh.dim > 2:
|
||||
wlist += (self.Wz,)
|
||||
self._Wsmooth = sp.vstack(wlist)
|
||||
return self._Wsmooth
|
||||
|
||||
@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 self._W
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmall(self, m):
|
||||
r = self.Wsmall * ( self.mapping * (m - self.mref) )
|
||||
return 0.5 * r.dot(r)
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmooth(self, m):
|
||||
if self.mrefInSmooth == True:
|
||||
r = self.Wsmooth * ( self.mapping * (m - self.mref) )
|
||||
elif self.mrefInSmooth == False:
|
||||
r = self.Wsmooth * ( self.mapping * m)
|
||||
return 0.5 * r.dot(r)
|
||||
|
||||
|
||||
class Sparse(Simple):
|
||||
|
||||
# set default values
|
||||
eps_p = 1e-1
|
||||
eps_q = 1e-1
|
||||
curModel = None # use a model to compute the weights
|
||||
gamma = 1.
|
||||
norms = [0., 2., 2., 2.]
|
||||
wght = 1.
|
||||
|
||||
def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
|
||||
Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
|
||||
|
||||
if isinstance(self.wght,float):
|
||||
self.wght = np.ones(self.regmesh.nC) * self.wght
|
||||
|
||||
@property
|
||||
def Wsmall(self):
|
||||
"""Regularization matrix Wsmall"""
|
||||
if getattr(self, 'curModel', None) is None:
|
||||
self.Rs = Utils.speye(self.regmesh.nC)
|
||||
|
||||
else:
|
||||
f_m = self.curModel - self.reg.mref
|
||||
self.rs = self.R(f_m , self.eps_p, self.norms[0])
|
||||
#print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs))
|
||||
self.Rs = Utils.sdiag( self.rs )
|
||||
|
||||
return Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma*self.wght)**0.5)*self.Rs
|
||||
|
||||
|
||||
@property
|
||||
def Wx(self):
|
||||
"""Regularization matrix Wx"""
|
||||
|
||||
if getattr(self, 'curModel', None) is None:
|
||||
self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
|
||||
|
||||
else:
|
||||
f_m = self.regmesh.cellDiffxStencil * self.curModel
|
||||
self.rx = self.R( f_m , self.eps_q, self.norms[1])
|
||||
self.Rx = Utils.sdiag( self.rx )
|
||||
|
||||
return Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.Rx*self.regmesh.cellDiffxStencil
|
||||
|
||||
@property
|
||||
def Wy(self):
|
||||
"""Regularization matrix Wy"""
|
||||
|
||||
if getattr(self, 'curModel', None) is None:
|
||||
self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0])
|
||||
|
||||
else:
|
||||
f_m = self.regmesh.cellDiffyStencil * self.curModel
|
||||
self.ry = self.R( f_m , self.eps_q, self.norms[2])
|
||||
self.Ry = Utils.sdiag( self.ry )
|
||||
|
||||
return Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.Ry*self.regmesh.cellDiffyStencil
|
||||
|
||||
@property
|
||||
def Wz(self):
|
||||
"""Regularization matrix Wz"""
|
||||
|
||||
if getattr(self, 'curModel', None) is None:
|
||||
self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0])
|
||||
|
||||
else:
|
||||
f_m = self.regmesh.cellDiffzStencil * self.curModel
|
||||
self.rz = self.R( f_m , self.eps_q, self.norms[3])
|
||||
self.Rz = Utils.sdiag( self.rz )
|
||||
|
||||
return Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.Rz*self.regmesh.cellDiffzStencil
|
||||
|
||||
@property
|
||||
def Wsmooth(self):
|
||||
"""Full smoothness regularization matrix W"""
|
||||
#if getattr(self, '_Wsmooth', None) is None:
|
||||
wlist = (self.Wx,)
|
||||
if self.regmesh.dim > 1:
|
||||
wlist += (self.Wy,)
|
||||
if self.regmesh.dim > 2:
|
||||
wlist += (self.Wz,)
|
||||
#self._Wsmooth = sp.vstack(wlist)
|
||||
return sp.vstack(wlist)
|
||||
|
||||
@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 self._W
|
||||
|
||||
def R(self, f_m , eps, exponent):
|
||||
|
||||
# Eta scaling is important for mix-norms...do not mess with it
|
||||
eta = (eps**(1.-exponent/2.))**0.5
|
||||
r = eta / (f_m**2.+ eps**2.)**((1.-exponent/2.)/2.)
|
||||
|
||||
|
||||
@@ -7,3 +7,4 @@ from CounterUtils import *
|
||||
import ModelBuilder
|
||||
import SolverUtils
|
||||
from coordutils import *
|
||||
from modelutils import *
|
||||
|
||||
@@ -0,0 +1,63 @@
|
||||
from matutils import mkvc, ndgrid
|
||||
import numpy as np
|
||||
|
||||
def surface2ind_topo(mesh, topo, gridLoc='CC'):
|
||||
# def genActiveindfromTopo(mesh, topo):
|
||||
"""
|
||||
Get active indices from topography
|
||||
"""
|
||||
|
||||
|
||||
if mesh.dim == 3:
|
||||
from scipy.interpolate import NearestNDInterpolator
|
||||
Ftopo = NearestNDInterpolator(topo[:,:2], topo[:,2])
|
||||
|
||||
if gridLoc == 'CC':
|
||||
XY = ndgrid(mesh.vectorCCx, mesh.vectorCCy)
|
||||
Zcc = mesh.gridCC[:,2].reshape((np.prod(mesh.vnC[:2]), mesh.nCz), order='F')
|
||||
|
||||
gridTopo = Ftopo(XY)
|
||||
actind = [gridTopo[ixy] <= Zcc[ixy,:] for ixy in range(np.prod(mesh.vnC[0]))]
|
||||
actind = np.hstack(actind)
|
||||
|
||||
elif gridLoc == 'N':
|
||||
|
||||
XY = ndgrid(mesh.vectorNx, mesh.vectorNy)
|
||||
gridTopo = Ftopo(XY).reshape(mesh.vnN[:2], order='F')
|
||||
|
||||
if mesh._meshType not in ['TENSOR', 'CYL', 'BASETENSOR']:
|
||||
raise NotImplementedError('Nodal surface2ind_topo not implemented for %s mesh'%mesh._meshType)
|
||||
|
||||
Nz = mesh.vectorNz[1:] # TODO: this will only work for tensor meshes
|
||||
actind = np.array([False]*mesh.nC).reshape(mesh.vnC, order='F')
|
||||
|
||||
for ii in range(mesh.nCx):
|
||||
for jj in range(mesh.nCy):
|
||||
actind[ii,jj,:] = [np.all(gridTopo[ii:ii+2, jj:jj+2] >= Nz[kk]) for kk in range(len(Nz)) ]
|
||||
|
||||
elif mesh.dim == 2:
|
||||
from scipy.interpolate import interp1d
|
||||
Ftopo = interp1d(topo[:,0], topo[:,1])
|
||||
|
||||
if gridLoc == 'CC':
|
||||
gridTopo = Ftopo(mesh.gridCC[:,0])
|
||||
actind = mesh.gridCC[:,1] <= gridTopo
|
||||
|
||||
elif gridLoc == 'N':
|
||||
|
||||
gridTopo = Ftopo(mesh.vectorNx)
|
||||
if mesh._meshType not in ['TENSOR', 'CYL', 'BASETENSOR']:
|
||||
raise NotImplementedError('Nodal surface2ind_topo not implemented for %s mesh'%mesh._meshType)
|
||||
|
||||
Ny = mesh.vectorNy[1:] # TODO: this will only work for tensor meshes
|
||||
actind = np.array([False]*mesh.nC).reshape(mesh.vnC, order='F')
|
||||
|
||||
for ii in range(mesh.nCx):
|
||||
actind[ii,:] = [np.all(gridTopo[ii:ii+2] > Ny[kk]) for kk in range(len(Ny)) ]
|
||||
|
||||
else:
|
||||
raise NotImplementedError('surface2ind_topo not implemented for 1D mesh')
|
||||
|
||||
return mkvc(actind)
|
||||
|
||||
|
||||
+1
-1
@@ -15,7 +15,7 @@ import Directives
|
||||
import Inversion
|
||||
import Tests
|
||||
|
||||
__version__ = '0.1.10'
|
||||
__version__ = '0.1.11'
|
||||
__author__ = 'Rowan Cockett'
|
||||
__license__ = 'MIT'
|
||||
__copyright__ = 'Copyright 2014 Rowan Cockett'
|
||||
|
||||
+2
-2
@@ -51,9 +51,9 @@ copyright = u'2013, SimPEG Developers'
|
||||
# built documents.
|
||||
#
|
||||
# The short X.Y version.
|
||||
version = '0.1.10'
|
||||
version = '0.1.11'
|
||||
# The full version, including alpha/beta/rc tags.
|
||||
release = '0.1.10'
|
||||
release = '0.1.11'
|
||||
|
||||
# The language for content autogenerated by Sphinx. Refer to documentation
|
||||
# for a list of supported languages.
|
||||
|
||||
@@ -20,8 +20,8 @@ INPUT:
|
||||
loc = Location of spheres [[x1,y1,z1],[x2,y2,z2]]
|
||||
radi = Radius of spheres [r1,r2]
|
||||
param = Conductivity of background and two spheres [m0,m1,m2]
|
||||
stype = survey type "pdp" (pole dipole) or "dpdp" (dipole dipole)
|
||||
dtype = Data type "appr" (app res) | "appc" (app cond) | "volt" (potential)
|
||||
surveyType = survey type 'pole-dipole' or 'dipole-dipole'
|
||||
unitType = Data type "appResistivity" | "appConductivity" | "volt"
|
||||
Created by @fourndo
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,25 @@
|
||||
.. _examples_Mesh_Basic_ForwardDC:
|
||||
|
||||
.. --------------------------------- ..
|
||||
.. ..
|
||||
.. THIS FILE IS AUTO GENEREATED ..
|
||||
.. ..
|
||||
.. SimPEG/Examples/__init__.py ..
|
||||
.. ..
|
||||
.. --------------------------------- ..
|
||||
|
||||
|
||||
Mesh: Basic Forward 2D DC Resistivity
|
||||
=====================================
|
||||
|
||||
2D DC forward modeling example with Tensor and Curvilinear Meshes
|
||||
|
||||
|
||||
.. plot::
|
||||
|
||||
from SimPEG import Examples
|
||||
Examples.Mesh_Basic_ForwardDC.run()
|
||||
|
||||
.. literalinclude:: ../../SimPEG/Examples/Mesh_Basic_ForwardDC.py
|
||||
:language: python
|
||||
:linenos:
|
||||
+8
-5
@@ -1,4 +1,4 @@
|
||||
.. _examples_Forward_BasicDirectCurrent:
|
||||
.. _examples_Utils_surface2ind_topo:
|
||||
|
||||
.. --------------------------------- ..
|
||||
.. ..
|
||||
@@ -8,14 +8,17 @@
|
||||
.. ..
|
||||
.. --------------------------------- ..
|
||||
|
||||
Forward BasicDirectCurrent
|
||||
==========================
|
||||
|
||||
Here we show how to use :code:`Utils.surface2ind_topo` to identify cells below
|
||||
a topographic surface.
|
||||
|
||||
|
||||
|
||||
.. plot::
|
||||
|
||||
from SimPEG import Examples
|
||||
Examples.Forward_BasicDirectCurrent.run()
|
||||
Examples.Utils_surface2ind_topo.run()
|
||||
|
||||
.. literalinclude:: ../../SimPEG/Examples/Forward_BasicDirectCurrent.py
|
||||
.. literalinclude:: ../../SimPEG/Examples/Utils_surface2ind_topo.py
|
||||
:language: python
|
||||
:linenos:
|
||||
Reference in New Issue
Block a user