From 86b8938d0282eedbdf9e27941a6a0b88a19d4b20 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 15 Oct 2014 11:19:24 -0700 Subject: [PATCH] 3 diagonal estimators implemented --- SimPEG/Utils/matutils.py | 26 ++++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/SimPEG/Utils/matutils.py b/SimPEG/Utils/matutils.py index dbb27bac..57d165ad 100644 --- a/SimPEG/Utils/matutils.py +++ b/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