Create function to generate linear forward operator (F)

Clean project
This commit is contained in:
D Fournier
2015-12-21 08:42:50 -08:00
parent e9f02e1d77
commit d6a2490a7b
16 changed files with 1408 additions and 4746 deletions
+95
View File
@@ -0,0 +1,95 @@
def fwr_MAG_data(mesh,B,M,rxLoc,model,flag):
"""
Forward model magnetic data using integral equation
INPUT:
xn, yn, zn = Mesh nodes location
B = Inducing field parameter [Binc, Bdecl, B0]
M = Magnetization matrix [Minc, Mdecl]
rxLox = Observation location informat [obsx, obsy, obsz]
model = Model associated with mesh
OUTPUT:
dobs =Observation array in format [obsx, obsy, obsz, data]
Created on Oct 7, 2015
@author: dominiquef
"""
#%%
from SimPEG import np, Utils, sp, mkvc
from get_T_mat import get_T_mat
xn = mesh.vectorNx;
yn = mesh.vectorNy;
zn = mesh.vectorNz;
mcell = (len(xn)-1) * (len(yn)-1) * (len(zn)-1)
ndata = rxLoc.shape[0]
# Convert declination from north to cartesian
Md = (450.-float(M[1]))%360.
# Create magnetization matrix
mx = np.cos(np.deg2rad(M[0])) * np.cos(np.deg2rad(Md))
my = np.cos(np.deg2rad(M[0])) * np.sin(np.deg2rad(Md))
mz = np.sin(np.deg2rad(M[0]))
Mx = Utils.sdiag(np.ones([mcell])*mx*B[2])
My = Utils.sdiag(np.ones([mcell])*my*B[2])
Mz = Utils.sdiag(np.ones([mcell])*mz*B[2])
#matplotlib.pyplot.spy(scipy.sparse.csr_matrix(Mx))
#plt.show()
Mxyz = sp.vstack((Mx,My,Mz));
#%% Create TMI projector
# Convert Bdecination from north to cartesian
D = (450.-float(B[1]))%360.
Ptmi = mkvc(np.r_[np.cos(np.deg2rad(B[0]))*np.cos(np.deg2rad(D)),np.cos(np.deg2rad(B[0]))*np.sin(np.deg2rad(D)),np.sin(np.deg2rad(B[0]))],2).T;
if flag=='tmi':
d = np.zeros(ndata)
elif flag=='xyz':
d = np.zeros((3,ndata))
# Loop through all observations and create forward operator (ndata-by-mcell)
print "Begin forward modeling " +str(int(ndata)) + " data points..."
# Add counter to dsiplay progress. Good for large problems
progress = -1;
for ii in range(ndata):
tx, ty, tz = get_T_mat(xn,yn,zn,rxLoc[ii,:])
Gxyz = np.vstack((tx,ty,tz))*Mxyz
if flag=='xyz':
d[:,ii] = Gxyz.dot(model)
elif flag=='tmi':
d[ii] = Ptmi.dot(Gxyz.dot(model))
#%%
# Forward operator
d_iter = np.floor(float(ii)/float(ndata)*10.);
if d_iter > progress:
arg = "Done " + str(d_iter*10) + " %"
print arg
progress = d_iter;
print "Done 100% ...forward modeling completed!!\n"
return d
+132
View File
@@ -0,0 +1,132 @@
def fwr_MAG_F(mesh,B,M,rxLoc,flag):
"""
Forward model magnetic data using integral equation
INPUT:
mesh = Mesh in SimPEG format
B = Inducing field parameter [Binc, Bdecl, B0]
M = Magnetization information
[OPTIONS]
1- [Minc, Mdecl] : Assumes uniform magnetization orientation
2- [mx1,mx2,..., my1,...,mz1] : cell-based defined magnetization direction
3- diag(M): Block diagonal matrix with [Mx, My, Mz] along the diagonal
rxLox = Observation location informat [obsx, obsy, obsz]
flag = 'tmi' | 'xyz' | 'full'
[OPTIONS]
1- tmi : Magnetization direction used and data are projected onto the
inducing field direction F.shape([ndata, nc])
2- xyz : Magnetization direction used and data are given in 3-components
F.shape([3*ndata, nc])
3- full: Full tensor matrix stored with shape([3*ndata, 3*nc])
OUTPUT:
F = Linear forward modeling operation
Created on Dec, 20th 2015
@author: dominiquef
"""
#%%
from SimPEG import np, Utils, sp, mkvc
from get_T_mat import get_T_mat
xn = mesh.vectorNx;
yn = mesh.vectorNy;
zn = mesh.vectorNz;
mcell = (len(xn)-1) * (len(yn)-1) * (len(zn)-1)
ndata = rxLoc.shape[0]
#%% Create TMI projector
# Convert Bdecination from north to cartesian
D = (450.-float(B[1]))%360.
Ptmi = mkvc(np.r_[np.cos(np.deg2rad(B[0]))*np.cos(np.deg2rad(D)),
np.cos(np.deg2rad(B[0]))*np.sin(np.deg2rad(D)),
np.sin(np.deg2rad(B[0]))],2).T;
# Pre-allocate space
if flag=='tmi' | flag == 'xyz':
# If assumes uniform magnetization direction
if len(M) == 2:
# Convert declination from north to cartesian
Md = (450.-float(M[1]))%360.
# Create magnetization matrix
mx = np.cos(np.deg2rad(M[0])) * np.cos(np.deg2rad(Md))
my = np.cos(np.deg2rad(M[0])) * np.sin(np.deg2rad(Md))
mz = np.sin(np.deg2rad(M[0]))
Mx = Utils.sdiag(np.ones([mcell])*mx*B[2])
My = Utils.sdiag(np.ones([mcell])*my*B[2])
Mz = Utils.sdiag(np.ones([mcell])*mz*B[2])
Mxyz = sp.vstack((Mx,My,Mz));
# Otherwise if given a vector 3*ncells
elif len(M) == mesh.nC * 3:
Mxyz = sp.spdiags(M,0,mesh.nC * 3,mesh.nC * 3)
if flag == 'tmi':
F = np.zeros(ndata, mesh.nC)
elif flag == 'xyz':
F = np.zeros(int(3*ndata), mesh.nC)
elif flag == 'full':
F = np.zeros(int(3*ndata), int(3*mesh.nC))
else:
print """Flag must be either 'tmi' | 'xyz' | 'full', please revised"""
return
# Loop through all observations and create forward operator (ndata-by-mcell)
print "Begin calculation of forward operator: " + flag
# Add counter to dsiplay progress. Good for large problems
progress = -1;
for ii in range(ndata):
tx, ty, tz = get_T_mat(xn,yn,zn,rxLoc[ii,:])
if flag=='tmi':
F[ii,:] = Ptmi.dot(np.vstack((tx,ty,tz)))*Mxyz
elif flag == 'xyz':
F[ii,:] = tx*Mxyz
F[ii+ndata,:] = ty*Mxyz
F[ii+2*ndata,:] = tz*Mxyz
elif flag == 'full':
F[ii,:] = tx
F[ii+ndata,:] = ty
F[ii+2*ndata,:] = tz
# Display progress
counter = np.floor(float(ii)/float(ndata)*10.);
if counter > progress:
arg = "Done " + str(counter*10) + " %"
print arg
progress = counter;
print "Done 100% ...forward modeling completed!!\n"
return F
+133
View File
@@ -0,0 +1,133 @@
'''
Created on Sep 27, 2015
@author: dominiquef
'''
def get_T_mat(xn,yn,zn,rxLoc):
"""
Load in the nodes of a tensor mesh and computes the magnetic tensor
for a given observation location [obsx, obsy, obsz]
OUTPUT:
Tx = [Txx Txy Txz]
Ty = [Tyx Tyy Tyz]
Tz = [Tzx Tzy Tzz]
where each elements have dimension 1-by-mcell.
Only the upper half 5 elements have to be computed since symetric.
Currently done as for-loops but will eventually be changed to vector
indexing, once the topography has been figured out.
"""
from SimPEG import np, mkvc
ncx = len(xn)-1
ncy = len(yn)-1
ncz = len(zn)-1
mcell = ncx*ncy*ncz
# Pre-allocate space for 1D array
Tx = np.zeros((1,3*mcell))
Ty = np.zeros((1,3*mcell))
Tz = np.zeros((1,3*mcell))
yn2,xn2,zn2 = np.meshgrid(yn[1:], xn[1:], zn[1:])
yn1,xn1,zn1 = np.meshgrid(yn[0:ncy], xn[0:ncx], zn[0:ncz])
yn2 = mkvc(yn2)
yn1 = mkvc(yn1)
zn2 = mkvc(zn2)
zn1 = mkvc(zn1)
xn2 = mkvc(xn2)
xn1 = mkvc(xn1)
#%%
#==============================================================================
dz2 = rxLoc[2] - zn1;
dz1 = rxLoc[2] - zn2;
dy2 = yn2 - rxLoc[1];
dy1 = yn1 - rxLoc[1];
dx2 = xn2 - rxLoc[0];
dx1 = xn1 - rxLoc[0];
R1 = ( dy2**2 + dx2**2 );
R2 = ( dy2**2 + dx1**2 );
R3 = ( dy1**2 + dx2**2 );
R4 = ( dy1**2 + dx1**2 );
arg1 = np.sqrt( dz2**2 + R2 );
arg2 = np.sqrt( dz2**2 + R1 );
arg3 = np.sqrt( dz1**2 + R1 );
arg4 = np.sqrt( dz1**2 + R2 );
arg5 = np.sqrt( dz2**2 + R3 );
arg6 = np.sqrt( dz2**2 + R4 );
arg7 = np.sqrt( dz1**2 + R4 );
arg8 = np.sqrt( dz1**2 + R3 );
Tx[0,0:mcell] = np.arctan2( dy1 * dz2 , ( dx2 * arg5 ) ) +\
- np.arctan2( dy2 * dz2 , ( dx2 * arg2 ) ) +\
np.arctan2( dy2 * dz1 , ( dx2 * arg3 ) ) +\
- np.arctan2( dy1 * dz1 , ( dx2 * arg8 ) ) +\
np.arctan2( dy2 * dz2 , ( dx1 * arg1 ) ) +\
- np.arctan2( dy1 * dz2 , ( dx1 * arg6 ) ) +\
np.arctan2( dy1 * dz1 , ( dx1 * arg7 ) ) +\
- np.arctan2( dy2 * dz1 , ( dx1 * arg4 ) );
Ty[0,0:mcell] = np.log( ( dz2 + arg2 ) / (dz1 + arg3 ) ) +\
-np.log( ( dz2 + arg1 ) / (dz1 + arg4 ) ) +\
np.log( ( dz2 + arg6 ) / (dz1 + arg7 ) ) +\
-np.log( ( dz2 + arg5 ) / (dz1 + arg8 ) );
Ty[0,mcell:2*mcell] = np.arctan2( dx1 * dz2 , ( dy2 * arg1 ) ) +\
- np.arctan2( dx2 * dz2 , ( dy2 * arg2 ) ) +\
np.arctan2( dx2 * dz1 , ( dy2 * arg3 ) ) +\
- np.arctan2( dx1 * dz1 , ( dy2 * arg4 ) ) +\
np.arctan2( dx2 * dz2 , ( dy1 * arg5 ) ) +\
- np.arctan2( dx1 * dz2 , ( dy1 * arg6 ) ) +\
np.arctan2( dx1 * dz1 , ( dy1 * arg7 ) ) +\
- np.arctan2( dx2 * dz1 , ( dy1 * arg8 ) );
R1 = (dy2**2 + dz1**2);
R2 = (dy2**2 + dz2**2);
R3 = (dy1**2 + dz1**2);
R4 = (dy1**2 + dz2**2);
Ty[0,2*mcell:] = np.log( ( dx1 + np.sqrt( dx1**2 + R1 ) ) / (dx2 + np.sqrt( dx2**2 + R1 ) ) ) +\
-np.log( ( dx1 + np.sqrt( dx1**2 + R2 ) ) / (dx2 + np.sqrt( dx2**2 + R2 ) ) ) +\
np.log( ( dx1 + np.sqrt( dx1**2 + R4 ) ) / (dx2 + np.sqrt( dx2**2 + R4 ) ) ) +\
-np.log( ( dx1 + np.sqrt( dx1**2 + R3 ) ) / (dx2 + np.sqrt( dx2**2 + R3 ) ) );
R1 = (dx2**2 + dz1**2);
R2 = (dx2**2 + dz2**2);
R3 = (dx1**2 + dz1**2);
R4 = (dx1**2 + dz2**2);
Tx[0,2*mcell:] = np.log( ( dy1 + np.sqrt( dy1**2 + R1 ) ) / (dy2 + np.sqrt( dy2**2 + R1 ) ) ) +\
-np.log( ( dy1 + np.sqrt( dy1**2 + R2 ) ) / (dy2 + np.sqrt( dy2**2 + R2 ) ) ) +\
np.log( ( dy1 + np.sqrt( dy1**2 + R4 ) ) / (dy2 + np.sqrt( dy2**2 + R4 ) ) ) +\
-np.log( ( dy1 + np.sqrt( dy1**2 + R3 ) ) / (dy2 + np.sqrt( dy2**2 + R3 ) ) );
Tz[0,2*mcell:] = -( Ty[0,mcell:2*mcell] + Tx[0,0:mcell] );
Tz[0,mcell:2*mcell] = Ty[0,2*mcell:];
Tx[0,mcell:2*mcell] = Ty[0,0:mcell];
Tz[0,0:mcell] = Tx[0,2*mcell:];
Tx = Tx/(4*np.pi);
Ty = Ty/(4*np.pi);
Tz = Tz/(4*np.pi);
return Tx,Ty,Tz
+117
View File
@@ -0,0 +1,117 @@
'''
Created on Jul 17, 2013
@author: dominiquef
'''
def get_UBC_mesh(meshfile):
""" Read UBC mesh file and extract parameters
Works for the condenced version (20 * 3) --> [20 20 20] """
fid = open(meshfile,'r')
from numpy import zeros
# Go through the log file and extract data and the last achieved misfit
for ii in range (1, 6):
line = fid.readline()
line = line.split(' ')
# First line: number of cells in i, j, k
if ii == 1:
numcell=[]
for jj in range(len(line)):
t = int(line[jj])
numcell.append(t)
nX = numcell[0]
nY = numcell[1]
nZ = numcell[2]
# Second line: origin coordinate (X,Y,Z)
elif ii==2:
origin = []
for jj in range(len(line)):
t = float(line[jj])
origin.append(t)
# Other lines for the xn, yn, zn (nodes location)
elif ii==3:
xn=zeros((nX+1,1), dtype=float)
xn[0] = origin[0]
count_entry = 0;
count = 0;
while (count<nX):
if line[count_entry].find('*') != -1:
ndx = line[count_entry].split('*')
for kk in range(int(ndx[0])):
xn[count+1] = xn[count] + (ndx[1])
count = count+1
count_entry=count_entry+1
else:
t = float(line[count_entry])
xn[count+1]= xn[count] +t
count = count+1;
count_entry=count_entry+1
elif ii==4:
yn=zeros((nY+1,1), dtype=float)
yn[0] = origin[0]
count_entry = 0;
count = 0;
while (count<nY):
if line[count_entry].find('*') != -1:
ndx = line[count_entry].split('*')
for kk in range(int(ndx[0])):
yn[count+1] = yn[count] + (ndx[1])
count = count+1
count_entry=count_entry+1
else:
t = float(line[count_entry])
yn[count+1]= yn[count] +t
count = count+1;
count_entry=count_entry+1
elif ii==5:
zn=zeros((nZ+1,1), dtype=float)
zn[0] = origin[0]
count_entry = 0;
count = 0;
while (count<nZ):
if line[count_entry].find('*') != -1:
ndx = line[count_entry].split('*')
for kk in range(int(ndx[0])):
zn[count+1] = zn[count] + (ndx[1])
count = count+1
count_entry=count_entry+1
else:
t = float(line[count_entry])
zn[count+1]= zn[count] +t
count = count+1;
count_entry=count_entry+1
fid.close();
return xn,yn,zn
+58
View File
@@ -0,0 +1,58 @@
'''
Created on Jul 17, 2013
@author: dominiquef
'''
def read_MAG_obs(obs_file):
"""Read input files for the lp_norm script"""
from numpy import zeros
fid = open(obs_file,'r')
# First line has the declination, inclination and amplitude of B0
line = fid.readline()
line = line.split()
Incl = float(line[0])
Decl = float(line[1])
B0 = float(line[2])
# Second line has the magnetization orientation and a flag
line = fid.readline()
line = line.split()
Minc = float(line[0])
Mdec = float(line[1])
FLAG = float(line[2])
# Third line has the number of rows
line = fid.readline()
line = line.split()
ndat = int(line[0])
# Pre-allocate space for obsx, obsy, obsz, data, uncert
obsx = zeros((ndat,1), dtype=float)
obsy = zeros((ndat,1), dtype=float)
obsz = zeros((ndat,1), dtype=float)
data = zeros((ndat,1), dtype=float)
unct = zeros((ndat,1), dtype=float)
for ii in range(ndat):
line = fid.readline()
line = line.split()
obsx[ii] = line[0]
obsy[ii] = line[1]
obsz[ii] = line[2]
if len(line)>3:
data[ii] = line[3]
if len(line)>4:
unct[ii] = line[4]
return Decl, Incl, B0, Mdec, Minc, obsx, obsy, obsz, data, unct
+52
View File
@@ -0,0 +1,52 @@
def read_MAGfwr_inp(input_file):
"""Read input files for forward modeling MAG data with integral form
INPUT:
input_file: File name containing the forward parameter
OUTPUT:
mshfile
obsfile
modfile
magfile
topofile
# All files should be in the working directory, otherwise the path must
# be specified.
Created on Jul 17, 2013
@author: dominiquef
"""
fid = open(input_file,'r')
line = fid.readline()
l_input = line.split('!')
mshfile = l_input[0].rstrip()
line = fid.readline()
l_input = line.split('!')
obsfile = l_input[0].rstrip()
line = fid.readline()
l_input = line.split('!')
modfile = l_input[0].rstrip()
line = fid.readline()
l_input = line.split('!')
if l_input=='null':
magfile = []
else:
magfile = l_input[0].rstrip()
line = fid.readline()
l_input = line.split('!')
if l_input=='null':
topofile = []
else:
topofile = l_input[0].rstrip()
return mshfile, obsfile, modfile, magfile, topofile