1. Added output for the edgeInnerProduct so we can use it in the derivative phase

2. The getMisfit may be working but needs a little bit more work

We need to go through our code to make sure that all vectors that can be multiplied by a matrix are set to
vec = mkvc(vec,2)

SimPEG/EMforward.py

@@ -1,13 +1,15 @@
 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):
+def getMisfit(m,mesh,forward,param):
 
     mu0   = 4*np.pi*1e-7
     omega = forward['omega'] #[param['indomega']]
     rhs   = forward['rhs']   #[:,param['indrhs']]
-    misfit = 0
+    mis   = 0
+    dmis  = m*0
 
     # Maxwell's system for E
     for i in range(len(omega)):
@@ -15,8 +17,8 @@ def getMisfit(m,mesh,forward):
             Curl  = mesh.edgeCurl
             #Grad  = mesh.nodalGrad
             sigma = np.exp(m)
-            Me,PP = mesh.getEdgeMass(sigma)
-            Mf    = 1/mu0 * mesh.getFaceMass()   # assume mu = mu0
+            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]))
@@ -30,20 +32,45 @@ def getMisfit(m,mesh,forward):
             mis  = mis + 0.5*(r.T*r)
             # get derivatives
             lam  = dsl.spsolve(A.T,P.T*r)
-            Gij  = PP.T*diag((PP*e)*mesh.vol) 
+            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)
-    ne   = np.sum(mesh.nE)
-    Q    = np.matrix(np.random.randn(ne,5))
-    P    = np.matrix(Q.T)
-    forward = {'omega':[1,2,3], 'rhs':Q,'projection':P}
+    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]
+    Px  = interpmat(x,y,z,xs,ys,zs)
+    xyz = mesh.gridEy
+    x   = xyz[:,0]; y = xyz[:,1]; z = xyz[:,2]
+    Py  = interpmat(x,y,z,xs,ys,zs)
+    xyz = mesh.gridEz
+    x   = xyz[:,0]; y = xyz[:,1]; z = xyz[:,2]
+    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),np.shape(Q,2),np.size(omega))
+    param   = {'dobs':dobs}
+
+
 
     m = np.ones(mesh.nC)
-    getMisfit(m,mesh,forward)
+    getMisfit(m,mesh,forward,param)

SimPEG/InnerProducts.py

@@ -169,13 +169,15 @@ def getEdgeInnerProduct(mesh, sigma):
     # Cell volume
     v = np.sqrt(mesh.vol)
     v3 = np.r_[v, v, v]
-    V = sdiag(v3)*Sigma*sdiag(v3)  # to keep symmetry
+    #V = sdiag(v3)*Sigma*sdiag(v3)  # to keep symmetry
 
-    A = P000.T*V*P000 + P001.T*V*P001 + P010.T*V*P010 + P011.T*V*P011 + P100.T*V*P100 + P101.T*V*P101 + P110.T*V*P110 + P111.T*V*P111
+    #A = P000.T*V*P000 + P001.T*V*P001 + P010.T*V*P010 + P011.T*V*P011 + P100.T*V*P100 + P101.T*V*P101 + P110.T*V*P110 + P111.T*V*P111
 
-    A = 0.125*A
-
-    return A
+    #A = 0.125*A
+    P = sp.vstack((sdiag(v3)*P000,sdiag(v3)*P001,sdiag(v3)*P010,sdiag(v3)*P011,
+                sdiag(v3)*P100,sdiag(v3)*P101,sdiag(v3)*P110,sdiag(v3)*P111))
+    A = 0.125* (P.T*Sigma*P)
+    return A, P