From c59c31602b07f8d1eec738068548ed11ba718600 Mon Sep 17 00:00:00 2001
From: pcain-usgs <pcain@usgs.gov>
Date: Mon, 25 Jan 2021 13:13:43 -0700
Subject: [PATCH] Capture log number metric
---
geomagio/adjusted/AdjustedMatrix.py | 3 ++-
geomagio/adjusted/Affine.py | 8 +++++---
geomagio/adjusted/Transform.py | 27 ++++++++++++++++-----------
test/adjusted_test/adjusted_test.py | 19 ++++++++++++++++++-
4 files changed, 41 insertions(+), 16 deletions(-)
diff --git a/geomagio/adjusted/AdjustedMatrix.py b/geomagio/adjusted/AdjustedMatrix.py
index 980d44d1d..5913d2a5a 100644
--- a/geomagio/adjusted/AdjustedMatrix.py
+++ b/geomagio/adjusted/AdjustedMatrix.py
@@ -1,6 +1,6 @@
from obspy import UTCDateTime
from pydantic import BaseModel
-from typing import Optional, Any
+from typing import Optional, Any, List
from .. import pydantic_utcdatetime
@@ -20,5 +20,6 @@ class AdjustedMatrix(BaseModel):
matrix: Any
pier_correction: float
+ metrics: List[float] = [0.0, 0.0]
starttime: Optional[UTCDateTime] = None
endtime: Optional[UTCDateTime] = None
diff --git a/geomagio/adjusted/Affine.py b/geomagio/adjusted/Affine.py
index 183b8a8e4..ffd395afd 100644
--- a/geomagio/adjusted/Affine.py
+++ b/geomagio/adjusted/Affine.py
@@ -245,7 +245,7 @@ class Affine(BaseModel):
times = self.get_times(readings)
Ms = []
weights = []
-
+ metrics = []
inputs = ordinates
for transform in self.transforms:
@@ -268,7 +268,7 @@ class Affine(BaseModel):
inputs = np.dot(
M, np.vstack([inputs[0], inputs[1], inputs[2], np.ones_like(inputs[0])])
)[0:3]
-
+ metrics.append(transform.metric)
Ms.append(M)
# compose affine transform matrices using reverse ordered matrices
@@ -277,7 +277,9 @@ class Affine(BaseModel):
[reading.pier_correction for reading in readings], weights=weights
)
- return AdjustedMatrix(matrix=M_composed, pier_correction=pier_correction)
+ return AdjustedMatrix(
+ matrix=M_composed, pier_correction=pier_correction, metrics=metrics
+ )
def get_weights(
self,
diff --git a/geomagio/adjusted/Transform.py b/geomagio/adjusted/Transform.py
index 9937cb583..dd98f1b67 100644
--- a/geomagio/adjusted/Transform.py
+++ b/geomagio/adjusted/Transform.py
@@ -13,8 +13,13 @@ class Transform(object):
Defaults to infinite(equal weighting)
"""
- def __init__(self, memory: float = np.inf):
+ def __init__(self, memory: float = np.inf, metric: float = 0.0):
self.memory = memory
+ self.metric = metric
+
+ def update_metric(self, sigma):
+ # indicates the maximum possible number of lost significant values
+ self.metric = np.log(abs(sigma[0] / sigma[-1]))
def get_weights(self, times: UTCDateTime, time: int = None) -> List[float]:
"""
@@ -176,7 +181,7 @@ class NoConstraints(LeastSq):
# regression matrix M that minimizes L2 norm
M_r, res, rank, sigma = spl.lstsq(ord_stacked.T, abs_stacked.T)
-
+ self.update_metric(sigma)
if rank < 3:
print("Poorly conditioned or singular matrix, returning NaNs")
return np.nan * np.ones((4, 4))
@@ -222,7 +227,7 @@ class ZRotationShear(LeastSq):
# regression matrix M that minimizes L2 norm
M_r, res, rank, sigma = spl.lstsq(ord_stacked.T, abs_stacked.T)
-
+ self.update_metric(sigma)
if rank < 3:
print("Poorly conditioned or singular matrix, returning NaNs")
return np.nan * np.ones((4, 4))
@@ -269,7 +274,7 @@ class ZRotationHscale(LeastSq):
# regression matrix M that minimizes L2 norm
M_r, res, rank, sigma = spl.lstsq(ord_stacked.T, abs_stacked.T)
-
+ self.update_metric(sigma)
if rank < 3:
print("Poorly conditioned or singular matrix, returning NaNs")
return np.nan * np.ones((4, 4))
@@ -316,7 +321,7 @@ class ZRotationHscaleZbaseline(LeastSq):
# regression matrix M that minimizes L2 norm
M_r, res, rank, sigma = spl.lstsq(ord_stacked.T, abs_stacked.T)
-
+ self.update_metric(sigma)
if rank < 3:
print("Poorly conditioned or singular matrix, returning NaNs")
return np.nan * np.ones((4, 4))
@@ -360,7 +365,7 @@ class RotationTranslation3D(SVD):
# Singular value decomposition, then rotation matrix from L&R eigenvectors
# (the determinant guarantees a rotation, and not a reflection)
U, S, Vh = np.linalg.svd(H)
-
+ self.update_metric(S)
if np.sum(S) < 3:
print("Poorly conditioned or singular matrix, returning NaNs")
return np.nan * np.ones((4, 4))
@@ -410,7 +415,7 @@ class Rescale3D(LeastSq):
# regression matrix M that minimizes L2 norm
M_r, res, rank, sigma = spl.lstsq(ord_stacked.T, abs_stacked.T)
-
+ self.update_metric(sigma)
if rank < 3:
print("Poorly conditioned or singular matrix, returning NaNs")
return np.nan * np.ones((4, 4))
@@ -465,7 +470,7 @@ class TranslateOrigins(LeastSq):
# regression matrix M that minimizes L2 norm
M_r, res, rank, sigma = spl.lstsq(ord_stacked.T, abs_stacked.T)
-
+ self.update_metric(sigma)
if rank < 3:
print("Poorly conditioned or singular matrix, returning NaNs")
return np.nan * np.ones((4, 4))
@@ -505,7 +510,7 @@ class ShearYZ(LeastSq):
abs_stacked = self.get_weighted_values(values=abs_stacked, weights=weights)
# regression matrix M that minimizes L2 norm
M_r, res, rank, sigma = spl.lstsq(ord_stacked.T, abs_stacked.T)
-
+ self.update_metric(sigma)
if rank < 3:
print("Poorly conditioned or singular matrix, returning NaNs")
return np.nan * np.ones((4, 4))
@@ -554,7 +559,7 @@ class RotationTranslationXY(SVD):
# Singular value decomposition, then rotation matrix from L&R eigenvectors
# (the determinant guarantees a rotation, and not a reflection)
U, S, Vh = np.linalg.svd(H)
-
+ self.update_metric(S)
if np.sum(S) < 2:
print("Poorly conditioned or singular matrix, returning NaNs")
return np.nan * np.ones((4, 4))
@@ -637,7 +642,7 @@ class QRFactorization(SVD):
# regression matrix M that minimizes L2 norm
M_r, res, rank, sigma = spl.lstsq(ord_stacked.T, abs_stacked.T)
-
+ self.update_metric(sigma)
if rank < 2:
print("Poorly conditioned or singular matrix, returning NaNs")
return np.nan * np.ones((4, 4))
diff --git a/test/adjusted_test/adjusted_test.py b/test/adjusted_test/adjusted_test.py
index d37c3aa04..eaebd4471 100644
--- a/test/adjusted_test/adjusted_test.py
+++ b/test/adjusted_test/adjusted_test.py
@@ -1,6 +1,6 @@
import json
import numpy as np
-from numpy.testing import assert_equal, assert_array_almost_equal
+from numpy.testing import assert_equal, assert_array_almost_equal, assert_array_less
from obspy.core import UTCDateTime
from geomagio.adjusted.SpreadsheetSummaryFactory import SpreadsheetSummaryFactory
@@ -235,6 +235,7 @@ def test_BOU201911202001_short_causal():
).calculate(readings=readings)
matrices = format_result([adjusted_matrix.matrix for adjusted_matrix in result])
+ metrics = [adjusted_matrix.metrics for adjusted_matrix in result]
expected_matrices = get_excpected_matrices("BOU", "short_causal")
for i in range(len(matrices)):
assert_array_almost_equal(
@@ -243,6 +244,7 @@ def test_BOU201911202001_short_causal():
decimal=3,
err_msg=f"Matrix {i} not equal",
)
+ assert_array_less(metrics[i], 5.0)
assert_equal(len(matrices), ((endtime - starttime) // update_interval) + 1)
@@ -265,6 +267,7 @@ def test_BOU201911202001_short_acausal():
)
matrices = format_result([adjusted_matrix.matrix for adjusted_matrix in result])
+ metrics = [adjusted_matrix.metrics for adjusted_matrix in result]
expected_matrices = get_excpected_matrices("BOU", "short_acausal")
for i in range(len(matrices)):
assert_array_almost_equal(
@@ -273,6 +276,7 @@ def test_BOU201911202001_short_acausal():
decimal=3,
err_msg=f"Matrix {i} not equal",
)
+ assert_array_less(metrics[i], 5.0)
assert_equal(len(matrices), ((endtime - starttime) // update_interval) + 1)
@@ -299,6 +303,7 @@ def test_BOU201911202001_infinite_weekly():
)
matrices = format_result([adjusted_matrix.matrix for adjusted_matrix in result])
+ metrics = [adjusted_matrix.metrics for adjusted_matrix in result]
expected_matrices = get_excpected_matrices("BOU", "inf_weekly")
for i in range(len(matrices)):
assert_array_almost_equal(
@@ -307,6 +312,7 @@ def test_BOU201911202001_infinite_weekly():
decimal=3,
err_msg=f"Matrix {i} not equal",
)
+ assert_array_less(metrics[i], 5.0)
assert_equal(len(matrices), ((endtime - starttime) // update_interval) + 1)
@@ -327,6 +333,7 @@ def test_BOU201911202001_infinite_one_interval():
)
matrices = format_result([adjusted_matrix.matrix for adjusted_matrix in result])
+ metrics = [adjusted_matrix.metrics for adjusted_matrix in result]
expected_matrices = get_excpected_matrices("BOU", "inf_one_interval")
for i in range(len(matrices)):
assert_array_almost_equal(
@@ -335,6 +342,7 @@ def test_BOU201911202001_infinite_one_interval():
decimal=3,
err_msg=f"Matrix {i} not equal",
)
+ assert_array_less(metrics[i], 5.0)
assert_equal(len(matrices), 1)
@@ -362,6 +370,7 @@ def test_CMO2015_causal():
)
matrices = format_result([adjusted_matrix.matrix for adjusted_matrix in result])
+ metrics = [adjusted_matrix.metrics for adjusted_matrix in result]
expected_matrices = get_excpected_matrices("CMO", "short_causal")
for i in range(len(matrices)):
assert_array_almost_equal(
@@ -370,6 +379,7 @@ def test_CMO2015_causal():
decimal=3,
err_msg=f"Matrix {i} not equal",
)
+ assert_array_less(metrics[i], 5.0)
assert_equal(len(matrices), ((endtime - starttime) // update_interval) + 1)
@@ -398,6 +408,7 @@ def test_CMO2015_acausal():
)
matrices = format_result([adjusted_matrix.matrix for adjusted_matrix in result])
+ metrics = [adjusted_matrix.metrics for adjusted_matrix in result]
expected_matrices = get_excpected_matrices("CMO", "short_acausal")
for i in range(len(matrices)):
assert_array_almost_equal(
@@ -406,6 +417,7 @@ def test_CMO2015_acausal():
decimal=3,
err_msg=f"Matrix {i} not equal",
)
+ assert_array_less(metrics[i], 5.0)
assert_equal(len(matrices), ((endtime - starttime) // update_interval) + 1)
@@ -438,6 +450,7 @@ def test_CMO2015_infinite_weekly():
)
matrices = format_result([adjusted_matrix.matrix for adjusted_matrix in result])
+ metrics = [adjusted_matrix.metrics for adjusted_matrix in result]
expected_matrices = get_excpected_matrices("CMO", "inf_weekly")
for i in range(len(matrices)):
assert_array_almost_equal(
@@ -446,6 +459,7 @@ def test_CMO2015_infinite_weekly():
decimal=3,
err_msg=f"Matrix {i} not equal",
)
+ assert_array_less(metrics[i], 5.0)
assert_equal(len(matrices), ((endtime - starttime) // update_interval) + 1)
@@ -474,6 +488,7 @@ def test_CMO2015_infinite_one_interval():
)
matrices = format_result([adjusted_matrix.matrix for adjusted_matrix in result])
+ metrics = [adjusted_matrix.metrics for adjusted_matrix in result]
expected_matrices = get_excpected_matrices("CMO", "inf_one_interval")
for i in range(len(matrices)):
assert_array_almost_equal(
@@ -482,4 +497,6 @@ def test_CMO2015_infinite_one_interval():
decimal=3,
err_msg=f"Matrix {i} not equal",
)
+ assert_array_less(metrics[i], 5.0)
+
assert_equal(len(matrices), 1)
--
GitLab