From c3476321a87f8da1289fa4520d2aef77fdc859c0 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 6 Aug 2013 18:07:21 -0700 Subject: [PATCH] Removed EMforward code for merging into master. (This is still being developed in branch: eldadsWork) --- SimPEG/EMforward.py | 79 ------------------------------------- SimPEG/interpmat.py | 94 --------------------------------------------- SimPEG/ops.py | 62 ------------------------------ 3 files changed, 235 deletions(-) delete mode 100644 SimPEG/EMforward.py delete mode 100644 SimPEG/interpmat.py delete mode 100644 SimPEG/ops.py diff --git a/SimPEG/EMforward.py b/SimPEG/EMforward.py deleted file mode 100644 index 9c332c9b..00000000 --- a/SimPEG/EMforward.py +++ /dev/null @@ -1,79 +0,0 @@ -import numpy as np -from utils import mkvc -import scipy.sparse.linalg.dsolve as dsl -from InnerProducts import getFaceInnerProduct, getEdgeInnerProduct - -def getMisfit(m,mesh,forward,param): - - mu0 = 4*np.pi*1e-7 - omega = forward['omega'] #[param['indomega']] - rhs = forward['rhs'] #[:,param['indrhs']] - mis = 0 - dmis = m*0 - - # Maxwell's system for E - for i in range(len(omega)): - for j in range(rhs.shape[1]): - Curl = mesh.edgeCurl - #Grad = mesh.nodalGrad - sigma = np.exp(m) - Me,PP = getEdgeInnerProduct(mesh,sigma) - Mf = 1/mu0 * getFaceInnerProduct(mesh) # assume mu = mu0 - - A = Curl.T * Mf * Curl - 1j * omega[i] * Me - b = mkvc(np.array(rhs[:,j])) - e = dsl.spsolve(A,b) - e = mkvc(e,2) - #print np.linalg.norm(A*e-b)/np.linalg.norm(b) - P = forward['projection'] - d = P*e - r = mkvc(d - param.dobs[i,j,:],2) - - mis = mis + 0.5*(r.T*r) - # get derivatives - lam = dsl.spsolve(A.T,P.T*r) - lam = mkvc(lam,2) - Gij = - 1j * omega[i] * PP.T*sp.diag((PP*e)*mesh.vol) - dmis = dmis - Gij.T*lam - - - return mis, dmis, d - - - -if __name__ == '__main__': - from TensorMesh import TensorMesh - from interpmat import interpmat - from scipy import sparse as sp - - h = [np.ones(7),np.ones(8),np.ones(9)] - mesh = TensorMesh(h) - xs = np.array([3.1,4.3,5.4,6.5]) - ys = np.array([3.2,4.1,5.4,6.2]) - zs = np.array([4.3,4.2,4.1,4.1]); - - xyz = mesh.gridEx - x = xyz[:,0]; y = xyz[:,1]; z = xyz[:,2] - x = list(set(x)); y = list(set(y)); z = list(set(z)) - Px = interpmat(x,y,z,xs,ys,zs) - xyz = mesh.gridEy - x = xyz[:,0]; y = xyz[:,1]; z = xyz[:,2] - x = list(set(x)); y = list(set(y)); z = list(set(z)) - Py = interpmat(x,y,z,xs,ys,zs) - xyz = mesh.gridEz - x = xyz[:,0]; y = xyz[:,1]; z = xyz[:,2] - x = list(set(x)); y = list(set(y)); z = list(set(z)) - Pz = interpmat(x,y,z,xs,ys,zs) - P = sp.hstack((Px,Py,Pz)) - - ne = np.sum(mesh.nE) - Q = np.matrix(np.random.randn(ne,5)) - omega = [1,2,3] - forward = {'omega':omega, 'rhs':Q,'projection':P} - dobs = np.ones([np.size(xs),5,np.size(omega)]) - param = {'dobs':dobs} - - - - m = np.ones(mesh.nC) - getMisfit(m,mesh,forward,param) diff --git a/SimPEG/interpmat.py b/SimPEG/interpmat.py deleted file mode 100644 index bfb80291..00000000 --- a/SimPEG/interpmat.py +++ /dev/null @@ -1,94 +0,0 @@ -from scipy import sparse as sp -import numpy as np - - -def interpmat(x,y,z,xr,yr,zr): -# -# This function does local linear interpolation -# computed for each receiver point in turn -# -# [Q] = linint(x,y,z,xr,yr,zr) -# Interpolation matrix -# - - nx = np.size(x) - ny = np.size(y) - nz = np.size(z) - - nps = np.size(xr) - - #Q = spalloc(np,nx*ny*nz,8*np); - Q = sp.lil_matrix((nps,nx*ny*nz)) - ind_x = np.array([0,0]) - ind_y = np.array([0,0]) - ind_z = np.array([0,0]) - dx, dy, dz = np.zeros(2), np.zeros(2), np.zeros(2) - for i in range(0, nps): - im = np.argmin(abs(xr[i]-x)) - print i,im - if xr[i] - x[im] >= 0: # Point on the left - ind_x[0] = im; ind_x[1] = im+1 - else: # Point on the right - ind_x[0] = im-1; ind_x[1] = im - - - dx[0] = xr[i] - x[ind_x[0]] - dx[1] = x[ind_x[1]] - xr[i] - - im = np.argmin(abs(yr[i] - y)) - if yr[i] - y[im] >= 0: # Point on the left - ind_y[0] = im; ind_y[1] = im+1 - else: # Point on the right - ind_y[0] = im-1; ind_y[1] = im - - - dy[0] = yr[i] - y[ind_y[0]] - dy[1] = y[ind_y[1]] - yr[i]; - - im = np.argmin(abs(zr[i] - z)); - if zr[i] -z[im] >= 0: # Point on the left - ind_z[0] = im; ind_z[1] = im+1 - else: # Point on the right - ind_z[0] = im-1; ind_z[1] = im; - - dz[0] = zr[i] - z[ind_z[0]]; dz[1] = z[ind_z[1]] - zr[i] - - Dx = x[ind_x[1]] - x[ind_x[0]] - Dy = y[ind_y[1]] - y[ind_y[0]] - Dz = z[ind_z[1]] - z[ind_z[0]] - #dv = Dx*Dy*Dz - - # Get the row in the matrix - v = np.zeros([nx, ny,nz]) - - v[ ind_x[0], ind_y[0], ind_z[0]] = (1-dx[0]/Dx)*(1-dy[0]/Dy)*(1-dz[0]/Dz) - v[ ind_x[0], ind_y[1], ind_z[0]] = (1-dx[0]/Dx)*(1-dy[1]/Dy)*(1-dz[0]/Dz) - v[ ind_x[1], ind_y[0], ind_z[0]] = (1-dx[1]/Dx)*(1-dy[0]/Dy)*(1-dz[0]/Dz) - v[ ind_x[1], ind_y[1], ind_z[0]] = (1-dx[1]/Dx)*(1-dy[1]/Dy)*(1-dz[0]/Dz) - v[ ind_x[0], ind_y[0], ind_z[1]] = (1-dx[0]/Dx)*(1-dy[0]/Dy)*(1-dz[1]/Dz) - v[ ind_x[0], ind_y[1], ind_z[1]] = (1-dx[0]/Dx)*(1-dy[1]/Dy)*(1-dz[1]/Dz) - v[ ind_x[1], ind_y[0], ind_z[1]] = (1-dx[1]/Dx)*(1-dy[0]/Dy)*(1-dz[1]/Dz) - v[ ind_x[1], ind_y[1], ind_z[1]] = (1-dx[1]/Dx)*(1-dy[1]/Dy)*(1-dz[1]/Dz) - - - print(np.shape(v.flatten('F'))) - print(np.shape(Q)) - - Q[i,:] = v.flatten('F') - - - return Q.tocsr() - - -if __name__ == '__main__': - - x = np.array([1.1, 2.1, 3.6, 4.9]) - y = np.array([1.2, 2.2, 3.3, 4.9, 5.6]) - z = np.array([0.8, 1.7, 4.9, 6.5]) - - xr = np.array([2.5,3.2]) - yr = np.array([2.4,3.6]) - zr = np.array([2.5,3.9]) - - A = interpmat(x,y,z,xr,yr,zr) - \ No newline at end of file diff --git a/SimPEG/ops.py b/SimPEG/ops.py deleted file mode 100644 index 4997cde9..00000000 --- a/SimPEG/ops.py +++ /dev/null @@ -1,62 +0,0 @@ -import numpy as np -from utils import mkvc -import scipy.sparse.linalg as spla -import scipy.sparse as sp - - -def matmul(A,B): - - # first check shape - if np.shape(A)[1] != np.shape(B)[0]: - print 'error in sizes' - return - - # Check types - sA = sp.issparse(A) - sB = sp.issparse(B) - - if ((sA == False) & (sB == True)): # doesno't work unless we trick it - return (B.T.dot(A.T)).T - else: - return A.dot(B) - - -def dot(A,B): - A = mkvc(A,1) - B = mkvc(B,1) - return np.dot(A,B) - -def inner(A,B): - A = mkvc(A,1) - B = mkvc(B,1) - return np.dot(A,B) - - -if __name__ == '__main__': - import numpy as np - from utils import mkvc - import scipy.sparse as sp - - # generate sparse and dense matrices - A = sp.rand(100, 200, density=0.05, format='csr', dtype=None) - B = sp.rand(200, 150, density=0.05, format='csr', dtype=None) - C = np.random.rand(200,150) - D = np.random.rand(150,100) - b = mkvc(np.arange(200),1) - c = np.reshape(b,(1,200)) - matmul(A,B) - matmul(A,C) - matmul(C,D) - matmul(D,A) - matmul(A,b) - dot(c,b) - dot(C,C) - print np.shape(c), np.shape(b)[0] - print matmul(c,b),dot(c,b) - - - - - - -