3 diagonal estimators implemented

SimPEG/Utils/matutils.py

@@ -341,26 +341,36 @@ def invPropertyTensor(M, tensor, returnMatrix=False):
     return T
 
 
-def diagEst(matFun, n, k=None, type='Probing1s'):
+def diagEst(matFun, n, k=None, type='Probing'):
     """ Based on Saad http://www-users.cs.umn.edu/~saad/PDF/umsi-2005-082.pdf, and http://www.cita.utoronto.ca/~niels/diagonal.pdf"""
 
     if k is None:
         k = np.floor(n/10.)
 
-    if type =='Probing1s':
-        getv = lambda n: np.random.random_integers(-1,high=1,size=n)
+    if type =='Ones':
+        def getv(n,i=None):
+            v = np.random.randn(n)
+            v[v<0]  = -1.
+            v[v>=0] = 1.
+            return v
 
-    elif type == 'ProbingRandn':
-        getv = lambda n: np.random.randn(n)
+    elif type == 'Random':
+        def getv(n,i=None):
+            return np.random.randn(n)
+
+    else: #if type == 'Probing':
+        def getv(n,i):
+            v = np.zeros(n)
+            v[i:n:k] = 1.
+            return v
 
     #d = np.zeros(n)
     Mv = np.zeros(n)
     vv = np.zeros(n)
 
     for i in range(0,k):
-        print k
-        vk = getv(n)
-        Mv += (matFun(vk))*vk
+        vk = getv(n,i)
+        Mv += matFun(vk)*vk
         vv += vk*vk
 
     d   = Mv/vv