From 8ac2cd1c4b646be5b31e267b4793ad8c47164a7e Mon Sep 17 00:00:00 2001
From: Arthus <arthus@gfz-potsdam.de>
Date: Tue, 22 Jun 2021 14:10:13 +0200
Subject: [PATCH] Add a routine to extract block diagonals and fix temporal
errors.
---
corbam/inversion.py | 31 ++++++++++++++++++++-----------
corbam/utils.py | 15 ++++++++++++++-
tests/test_inversion.py | 4 ++--
3 files changed, 36 insertions(+), 14 deletions(-)
diff --git a/corbam/inversion.py b/corbam/inversion.py
index b893b4c..03999da 100644
--- a/corbam/inversion.py
+++ b/corbam/inversion.py
@@ -45,7 +45,8 @@ from scipy import linalg
import pyfield
-from corbam.utils import scaling, read_data, dif2nez, nez2dif
+from corbam.utils import scaling, read_data, dif2nez, nez2dif, \
+ extractBlockDiag, blockDiag
__all__ = ['grad_d', 'grad_i', 'grad_f', 'Inversion']
@@ -218,7 +219,7 @@ class Inversion:
self.idx_NEZ_F))
# Temporal errors
- errs_T_full = data[['dt', 'dt', 'dt']].loc[self.idx_full]**2
+ errs_T_full = data[['dt']].loc[self.idx_full]**2
errs_T_full = np.asfortranarray(errs_T_full).ravel()
self.errs_T_full = errs_T_full
@@ -300,15 +301,14 @@ class Inversion:
self.cov_NEZ_full += lapl_appr
# Noisy input correction term for complete records
- corr_full = 2*np.diagonal((alph/theta_Dip)**2
- * cov_NEZ_obs_Dip[self.idx_NEZ_full[:, None],
- self.idx_NEZ_full[None, :]]
- + (lamb/theta_WoDip)**2
- * cov_NEZ_obs_WoDip[self.idx_NEZ_full[:,
- None],
- self.idx_NEZ_full[None,
- :]])
- self.corr_full = np.diag(corr_full * self.errs_T_full)
+ corr_full = \
+ 2*extractBlockDiag((alph/theta_Dip)**2
+ * cov_NEZ_obs_Dip[self.idx_NEZ_full[:, None],
+ self.idx_NEZ_full[None, :]]
+ + (lamb/theta_WoDip)**2
+ * cov_NEZ_obs_WoDip[self.idx_NEZ_full[:, None],
+ self.idx_NEZ_full[None, :]])
+ self.corr_full = blockDiag(corr_full * self.errs_T_full[:, None, None])
self.cov_NEZ_full += noisy * self.corr_full
# Invert covariance using Cholesky
L = linalg.cholesky(self.cov_NEZ_full, lower=True)
@@ -383,6 +383,9 @@ class Inversion:
+ (lamb/theta_WoDip)**2 \
* cov_NEZ_obs_WoDip[self.idx_NEZ_icmp[:, None],
self.idx_NEZ_icmp[None, :]]
+
+ corr_icmp2 = extractBlockDiag(corr_icmp)
+
corr_icmp = corr_icmp.reshape(self.n_icmp, 3,
self.n_icmp, 3)
@@ -392,6 +395,12 @@ class Inversion:
self.grad_DIF)
corr_icmp = 2*np.diagonal(corr_icmp)
+ corr_icmp2 = np.einsum('ij,ijl,il->i',
+ self.grad_DIF,
+ corr_icmp2,
+ self.grad_DIF)
+ print(np.allclose(corr_icmp, 2*corr_icmp2))
+
self.corr_icmp = np.diag(corr_icmp * self.errs_T_icmp)
self.cov_DIF_icmp += noisy * self.corr_icmp
diff --git a/corbam/utils.py b/corbam/utils.py
index ed442d8..7025133 100644
--- a/corbam/utils.py
+++ b/corbam/utils.py
@@ -35,7 +35,20 @@ import pyfield
path = Path(__file__).parent
__all__ = ['read_data', 'load', 'scaling', 'scale_coeffs', 'scale_vars',
- 'scale_covs', 'nez2dif', 'dif2nez']
+ 'scale_covs', 'nez2dif', 'dif2nez', 'extractBlockDiag',
+ 'blockDiag']
+
+
+def extractBlockDiag(A, dim=3):
+ n = int(A.shape[0] / dim)
+ return np.copy(A.reshape(n, dim, n, dim)[np.arange(n), :, np.arange(n), :])
+
+
+def blockDiag(B):
+ n, dim, _ = B.shape
+ ret = np.zeros((n, dim, n, dim))
+ ret[np.arange(n), :, np.arange(n), :] = B
+ return ret.reshape(n*dim, n*dim)
def read_data(fname):
diff --git a/tests/test_inversion.py b/tests/test_inversion.py
index a4acdd8..7eef99c 100644
--- a/tests/test_inversion.py
+++ b/tests/test_inversion.py
@@ -24,10 +24,10 @@ from corbam.inversion import Inversion
from corbam.utils import load
-pars = load('../tests/test_parfile.py')
+pars = load('./tests/test_parfile.py')
inv = Inversion(pars.data)
-with np.load('../dat/benchmark.npz') as fh:
+with np.load('./dat/benchmark.npz') as fh:
times = fh['times']
loc = fh['loc']
x = fh['x']
--
GitLab