From 35efd9039e2d33487816cbee1a7f2d59da787f1a Mon Sep 17 00:00:00 2001
From: pcain-usgs <pcain@usgs.gov>
Date: Fri, 22 Jan 2021 17:03:19 -0700
Subject: [PATCH] Add weights to synthetic testing
---
etc/adjusted/synthetic.json | 136 ++++++++++++++--------------
geomagio/adjusted/Transform.py | 78 +++++++---------
test/adjusted_test/adjusted_test.py | 93 ++++++-------------
3 files changed, 128 insertions(+), 179 deletions(-)
diff --git a/etc/adjusted/synthetic.json b/etc/adjusted/synthetic.json
index f93e5cd40..c3dd5594f 100644
--- a/etc/adjusted/synthetic.json
+++ b/etc/adjusted/synthetic.json
@@ -136,22 +136,22 @@
"results": {
"NoConstraints": [
[
- 0.5364919027495764,
- -0.20707814501736455,
- 0.3372529089344527,
- -1.7929218549826358
+ 0.5364919027495736,
+ -0.20707814501736121,
+ 0.33725290893445425,
+ -1.7929218549826373
],
[
- 0.532504593054401,
- 0.8470924780256476,
- -0.3269654921326787,
- -4.004987214867758
+ 0.5325045930544003,
+ 0.847092478025647,
+ -0.32696549213267934,
+ -4.004987214867749
],
[
- -0.2588476812717039,
- 0.4215661361680678,
- 0.6706148784369707,
- -2.9215661361680687
+ -0.25884768127170443,
+ 0.4215661361680674,
+ 0.6706148784369698,
+ -2.921566136168065
],
[
0,
@@ -162,22 +162,22 @@
],
"ZRotationShear": [
[
- 0.623705221351958,
- -0.16841150209489725,
+ 0.5587910280835219,
+ -0.1593988978820583,
0,
- -1.2028974670750636
+ -0.9594702768899068
],
[
- 0.4479515917846801,
- 0.8096053056474853,
+ 0.5108856705906101,
+ 0.8008676180137191,
0,
- -4.577013754155448
+ -4.813015546664324
],
[
0,
0,
- 0.6098875103018411,
- -1.809849332552786
+ 0.7036074465845426,
+ -2.273788534881388
],
[
0,
@@ -188,22 +188,22 @@
],
"ZRotationHscale": [
[
- 0.6561504712144426,
- -0.31217620673287216,
+ 0.6395521012137755,
+ -0.3193786953272828,
0,
- -0.7215676859558748
+ -0.6145604575448345
],
[
- 0.31217620673287216,
- 0.6561504712144426,
+ 0.3193786953272828,
+ 0.6395521012137755,
0,
- -3.563702076556801
+ -3.550295624100921
],
[
0,
0,
- 0.6098875103018406,
- -1.809849332552785
+ 0.7036074465845424,
+ -2.2737885348813887
],
[
0,
@@ -214,14 +214,14 @@
],
"ZRotationHscaleZbaseline": [
[
- 0.016821600812797096,
- -0.01778941712811468,
+ 0.053926953353156136,
+ 0.03582167951143748,
0,
0
],
[
- 0.01778941712811468,
- 0.016821600812797096,
+ -0.03582167951143748,
+ 0.053926953353156136,
0,
0
],
@@ -229,7 +229,7 @@
0,
0,
1,
- -2.9674214586052265
+ -3.2332029498017487
],
[
0,
@@ -240,22 +240,22 @@
],
"RotationTranslation3D": [
[
- 0.8054252079416996,
- -0.31075264909884437,
- 0.5047009267775519,
- -2.662452958809875
+ 0.8100811050998874,
+ -0.3030939677547838,
+ 0.5018990435045756,
+ -2.8734904692082317
],
[
- 0.5054022574334629,
- 0.804909769980755,
- -0.3109482599576447,
- -3.803929616959254
+ 0.4979757275740692,
+ 0.8075341522117618,
+ -0.3160834822617412,
+ -3.7722406528522203
],
[
- -0.309610711367079,
- 0.5055225546576319,
- 0.8053497092190405,
- -3.5091860177384517
+ -0.3094976218119015,
+ 0.5059867979723319,
+ 0.8051015975456041,
+ -3.511977037527143
],
[
0,
@@ -266,21 +266,21 @@
],
"Rescale3D": [
[
- 0.02758794737157335,
+ 0.10754195315849033,
0,
0,
0
],
[
0,
- 0.010807638086966475,
+ 0.012787625579547376,
0,
0
],
[
0,
0,
- 0.03534928189623147,
+ 0.032798154583522594,
0
],
[
@@ -295,19 +295,19 @@
1,
0,
0,
- -2.9812470914161855
+ -3.0518410718922895
],
[
0,
1,
0,
- -4.02728108273549
+ -3.8667699341149695
],
[
0,
0,
1,
- -2.9674214586052265
+ -3.2332029498017487
],
[
0,
@@ -324,14 +324,14 @@
0
],
[
- -1.3484311868544543,
+ -0.6692050417082805,
1,
0,
0
],
[
- 0.5519469765286653,
- -0.10894882663577743,
+ 0.35315324212931615,
+ -0.06394207042397665,
1,
0
],
@@ -344,22 +344,22 @@
],
"RotationTranslationXY": [
[
- 0.9030081198758902,
- -0.4296234810135611,
+ 0.8946494196067308,
+ -0.44676886193795917,
0,
- -0.9742971123340419
+ -0.9794537295011263
],
[
- 0.4296234810135612,
- 0.9030081198758901,
+ 0.4467688619379593,
+ 0.894649419606731,
0,
- -4.897795293523813
+ -4.958977628364134
],
[
0,
0,
1,
- -2.9674214586052274
+ -3.233202949801749
],
[
0,
@@ -370,22 +370,22 @@
],
"QRFactorization": [
[
- 0.8122230186734241,
- -0.22848881373914198,
+ 0.7380347209075413,
+ -0.2626732024800775,
0,
- -1.51463869998581
+ -1.1628964087154696
],
[
- 0.5833470390231107,
- 1.0670858952888813,
+ 0.6747627366229761,
+ 1.1147956698369132,
0,
- -6.006314806000902
+ -6.5703958608516055
],
[
0,
0,
1,
- -2.9674214586052274
+ -3.233202949801749
],
[
0,
diff --git a/geomagio/adjusted/Transform.py b/geomagio/adjusted/Transform.py
index ad79e5d24..0e2dcd58f 100644
--- a/geomagio/adjusted/Transform.py
+++ b/geomagio/adjusted/Transform.py
@@ -98,11 +98,8 @@ class LeastSq(Transform):
if weights is None:
return values
weights = np.sqrt(weights)
- return (
- values[0] * weights,
- values[1] * weights,
- values[2] * weights,
- )
+ weights = np.vstack((weights, weights, weights)).T.ravel()
+ return values * weights
class SVD(Transform):
@@ -146,8 +143,6 @@ class NoConstraints(LeastSq):
weights: List[float],
) -> np.array:
"""Calculates affine with no constraints using least squares."""
- ordinates = self.get_weighted_values(ordinates, weights)
- absolutes = self.get_weighted_values(absolutes, weights)
# LHS, or dependent variables
#
# [A[0,0], A[1,0], A[2,0], A[0,1], A[1,1], A[2,1], ...]
@@ -176,6 +171,9 @@ class NoConstraints(LeastSq):
ord_stacked[10, 2::3] = ordinates[2]
ord_stacked[11, 2::3] = 1.0
+ ord_stacked = self.get_weighted_values(ord_stacked, weights)
+ abs_stacked = self.get_weighted_values(abs_stacked, weights)
+
# regression matrix M that minimizes L2 norm
M_r, res, rank, sigma = spl.lstsq(ord_stacked.T, abs_stacked.T)
@@ -200,8 +198,6 @@ class ZRotationShear(LeastSq):
weights: List[float],
) -> np.array:
"""Calculates affine using least squares, constrained to rotate about the Z axis."""
- ordinates = self.get_weighted_values(ordinates, weights)
- absolutes = self.get_weighted_values(absolutes, weights)
# LHS, or dependent variables
#
abs_stacked = self.get_stacked_absolutes(absolutes)
@@ -221,6 +217,9 @@ class ZRotationShear(LeastSq):
ord_stacked[6, 2::3] = ordinates[2]
ord_stacked[7, 2::3] = 1.0
+ ord_stacked = self.get_weighted_values(ord_stacked, weights)
+ abs_stacked = self.get_weighted_values(abs_stacked, weights)
+
# regression matrix M that minimizes L2 norm
M_r, res, rank, sigma = spl.lstsq(ord_stacked.T, abs_stacked.T)
@@ -246,8 +245,6 @@ class ZRotationHscale(LeastSq):
absolutes: Tuple[List[float], List[float], List[float]],
weights: List[float],
) -> np.array:
- ordinates = self.get_weighted_values(ordinates, weights)
- absolutes = self.get_weighted_values(absolutes, weights)
# LHS, or dependent variables
#
abs_stacked = self.get_stacked_absolutes(absolutes)
@@ -267,6 +264,9 @@ class ZRotationHscale(LeastSq):
ord_stacked[4, 2::3] = ordinates[2]
ord_stacked[5, 2::3] = 1.0
+ ord_stacked = self.get_weighted_values(ord_stacked, weights)
+ abs_stacked = self.get_weighted_values(abs_stacked, weights)
+
# regression matrix M that minimizes L2 norm
M_r, res, rank, sigma = spl.lstsq(ord_stacked.T, abs_stacked.T)
@@ -344,15 +344,14 @@ class RotationTranslation3D(SVD):
# generate weighted "covariance" matrix
H = np.dot(
self.get_stacked_values(
- values=ordinates,
- weighted_values=weighted_ordinates,
- ndims=3,
+ values=ordinates, weighted_values=weighted_ordinates, ndims=3,
+ ),
+ np.dot(
+ np.diag(weights),
+ self.get_stacked_values(
+ values=absolutes, weighted_values=weighted_absolutes, ndims=3,
+ ).T,
),
- self.get_stacked_values(
- values=absolutes,
- weighted_values=weighted_absolutes,
- ndims=3,
- ).T,
)
# Singular value decomposition, then rotation matrix from L&R eigenvectors
# (the determinant guarantees a rotation, and not a reflection)
@@ -388,8 +387,6 @@ class Rescale3D(LeastSq):
absolutes: Tuple[List[float], List[float], List[float]],
weights: List[float],
) -> np.array:
- absolutes = self.get_weighted_values(values=absolutes, weights=weights)
- ordinates = self.get_weighted_values(values=ordinates, weights=weights)
abs_stacked = self.get_stacked_absolutes(absolutes=absolutes)
# RHS, or independent variables
# (reduces degrees of freedom by 13:
@@ -404,6 +401,9 @@ class Rescale3D(LeastSq):
ord_stacked[1, 1::3] = ordinates[1]
ord_stacked[2, 2::3] = ordinates[2]
+ abs_stacked = self.get_weighted_values(values=abs_stacked, weights=weights)
+ ord_stacked = self.get_weighted_values(values=ord_stacked, weights=weights)
+
# regression matrix M that minimizes L2 norm
M_r, res, rank, sigma = spl.lstsq(ord_stacked.T, abs_stacked.T)
@@ -483,8 +483,6 @@ class ShearYZ(LeastSq):
absolutes: Tuple[List[float], List[float], List[float]],
weights: List[float],
) -> np.array:
- ordinates = self.get_weighted_values(values=ordinates, weights=weights)
- # return generate_affine_7(ord_hez, abs_xyz, weights)
# RHS, or independent variables
# (reduces degrees of freedom by 13:
# - 2 for making x independent of y,z;
@@ -498,9 +496,10 @@ class ShearYZ(LeastSq):
ord_stacked[2, 0::3] = ordinates[0]
ord_stacked[2, 1::3] = ordinates[1]
ord_stacked[2, 2::3] = 1.0
-
- # regression matrix M that minimizes L2 norm
abs_stacked = self.get_stacked_absolutes(absolutes=absolutes)
+ ord_stacked = self.get_weighted_values(values=ord_stacked, weights=weights)
+ 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)
if rank < 3:
@@ -534,16 +533,12 @@ class RotationTranslationXY(SVD):
# generate weighted "covariance" matrix
H = np.dot(
self.get_stacked_values(
- values=ordinates,
- weighted_values=weighted_ordinates,
- ndims=2,
+ values=ordinates, weighted_values=weighted_ordinates, ndims=2,
),
np.dot(
np.diag(weights),
self.get_stacked_values(
- values=absolutes,
- weighted_values=weighted_absolutes,
- ndims=2,
+ values=absolutes, weighted_values=weighted_absolutes, ndims=2,
).T,
),
)
@@ -588,12 +583,7 @@ class QRFactorization(SVD):
if weights is None:
return values
weights = np.sqrt(weights)
- return np.array(
- [
- values[0] * weights,
- values[1] * weights,
- ]
- )
+ return np.array([values[0] * weights, values[1] * weights,])
def calculate(
self,
@@ -611,25 +601,19 @@ class QRFactorization(SVD):
# return generate_affine_9(ord_hez, abs_xyz, weights)
# LHS, or dependent variables
abs_stacked = self.get_stacked_values(
- values=absolutes,
- weighted_values=weighted_absolutes,
- ndims=2,
+ values=absolutes, weighted_values=weighted_absolutes, ndims=2,
)
# RHS, or independent variables
ord_stacked = self.get_stacked_values(
- values=ordinates,
- weighted_values=weighted_ordinates,
- ndims=2,
+ values=ordinates, weighted_values=weighted_ordinates, ndims=2,
)
abs_stacked = self.get_weighted_values_lstsq(
- values=abs_stacked,
- weights=weights,
+ values=abs_stacked, weights=weights,
)
ord_stacked = self.get_weighted_values_lstsq(
- values=ord_stacked,
- weights=weights,
+ values=ord_stacked, weights=weights,
)
# regression matrix M that minimizes L2 norm
diff --git a/test/adjusted_test/adjusted_test.py b/test/adjusted_test/adjusted_test.py
index 503806fb4..0e76e9725 100644
--- a/test/adjusted_test/adjusted_test.py
+++ b/test/adjusted_test/adjusted_test.py
@@ -54,7 +54,8 @@ def get_sythetic_variables():
variables = data["variables"]
ordinates = np.array([variables["h_ord"], variables["e_ord"], variables["z_ord"]])
absolutes = np.array([variables["x_abs"], variables["y_abs"], variables["z_abs"]])
- return ordinates, absolutes
+ weights = np.arange(0, len(ordinates[0]))
+ return ordinates, absolutes, weights
def get_expected_synthetic_result(key):
@@ -64,12 +65,10 @@ def get_expected_synthetic_result(key):
def test_NoConstraints_synthetic():
- ordinates, absolutes = get_sythetic_variables()
+ ordinates, absolutes, weights = get_sythetic_variables()
assert_array_almost_equal(
NoConstraints().calculate(
- ordinates=ordinates,
- absolutes=absolutes,
- weights=None,
+ ordinates=ordinates, absolutes=absolutes, weights=weights,
),
get_expected_synthetic_result("NoConstraints"),
decimal=3,
@@ -77,12 +76,10 @@ def test_NoConstraints_synthetic():
def test_ZRotationShear_synthetic():
- ordinates, absolutes = get_sythetic_variables()
+ ordinates, absolutes, weights = get_sythetic_variables()
assert_array_almost_equal(
ZRotationShear().calculate(
- ordinates=ordinates,
- absolutes=absolutes,
- weights=None,
+ ordinates=ordinates, absolutes=absolutes, weights=weights,
),
get_expected_synthetic_result("ZRotationShear"),
decimal=3,
@@ -90,12 +87,10 @@ def test_ZRotationShear_synthetic():
def test_ZRotationHscale_synthetic():
- ordinates, absolutes = get_sythetic_variables()
+ ordinates, absolutes, weights = get_sythetic_variables()
assert_array_almost_equal(
ZRotationHscale().calculate(
- ordinates=ordinates,
- absolutes=absolutes,
- weights=None,
+ ordinates=ordinates, absolutes=absolutes, weights=weights,
),
get_expected_synthetic_result("ZRotationHscale"),
decimal=3,
@@ -103,12 +98,10 @@ def test_ZRotationHscale_synthetic():
def test_ZRotationHscaleZbaseline_synthetic():
- ordinates, absolutes = get_sythetic_variables()
+ ordinates, absolutes, weights = get_sythetic_variables()
assert_array_almost_equal(
ZRotationHscaleZbaseline().calculate(
- ordinates=ordinates,
- absolutes=absolutes,
- weights=None,
+ ordinates=ordinates, absolutes=absolutes, weights=weights,
),
get_expected_synthetic_result("ZRotationHscaleZbaseline"),
decimal=3,
@@ -116,12 +109,10 @@ def test_ZRotationHscaleZbaseline_synthetic():
def test_RotationTranslation3D_synthetic():
- ordinates, absolutes = get_sythetic_variables()
+ ordinates, absolutes, weights = get_sythetic_variables()
assert_array_almost_equal(
RotationTranslation3D().calculate(
- ordinates=ordinates,
- absolutes=absolutes,
- weights=None,
+ ordinates=ordinates, absolutes=absolutes, weights=weights,
),
get_expected_synthetic_result("RotationTranslation3D"),
decimal=3,
@@ -129,12 +120,10 @@ def test_RotationTranslation3D_synthetic():
def test_Rescale3D_synthetic():
- ordinates, absolutes = get_sythetic_variables()
+ ordinates, absolutes, weights = get_sythetic_variables()
assert_array_almost_equal(
Rescale3D().calculate(
- ordinates=ordinates,
- absolutes=absolutes,
- weights=None,
+ ordinates=ordinates, absolutes=absolutes, weights=weights,
),
get_expected_synthetic_result("Rescale3D"),
decimal=3,
@@ -142,12 +131,10 @@ def test_Rescale3D_synthetic():
def test_TranslateOrigins_synthetic():
- ordinates, absolutes = get_sythetic_variables()
+ ordinates, absolutes, weights = get_sythetic_variables()
assert_array_almost_equal(
TranslateOrigins().calculate(
- ordinates=ordinates,
- absolutes=absolutes,
- weights=None,
+ ordinates=ordinates, absolutes=absolutes, weights=weights,
),
get_expected_synthetic_result("TranslateOrigins"),
decimal=3,
@@ -155,25 +142,19 @@ def test_TranslateOrigins_synthetic():
def test_ShearYZ_synthetic():
- ordinates, absolutes = get_sythetic_variables()
+ ordinates, absolutes, weights = get_sythetic_variables()
assert_array_almost_equal(
- ShearYZ().calculate(
- ordinates=ordinates,
- absolutes=absolutes,
- weights=None,
- ),
+ ShearYZ().calculate(ordinates=ordinates, absolutes=absolutes, weights=weights,),
get_expected_synthetic_result("ShearYZ"),
decimal=3,
)
def test_RotationTranslationXY_synthetic():
- ordinates, absolutes = get_sythetic_variables()
+ ordinates, absolutes, weights = get_sythetic_variables()
assert_array_almost_equal(
RotationTranslationXY().calculate(
- ordinates=ordinates,
- absolutes=absolutes,
- weights=None,
+ ordinates=ordinates, absolutes=absolutes, weights=weights,
),
get_expected_synthetic_result("RotationTranslationXY"),
decimal=3,
@@ -181,12 +162,10 @@ def test_RotationTranslationXY_synthetic():
def test_QRFactorization_synthetic():
- ordinates, absolutes = get_sythetic_variables()
+ ordinates, absolutes, weights = get_sythetic_variables()
assert_array_almost_equal(
QRFactorization().calculate(
- ordinates=ordinates,
- absolutes=absolutes,
- weights=None,
+ ordinates=ordinates, absolutes=absolutes, weights=weights,
),
get_expected_synthetic_result("QRFactorization"),
decimal=3,
@@ -259,9 +238,7 @@ def test_BOU201911202001_short_acausal():
endtime=endtime,
update_interval=update_interval,
acausal=True,
- ).calculate(
- readings=readings,
- )
+ ).calculate(readings=readings,)
matrices = format_result([adjusted_matrix.matrix for adjusted_matrix in result])
expected_matrices = get_excpected_matrices("BOU", "short_acausal")
@@ -293,9 +270,7 @@ def test_BOU201911202001_infinite_weekly():
RotationTranslationXY(memory=np.inf),
TranslateOrigins(memory=np.inf),
],
- ).calculate(
- readings=readings,
- )
+ ).calculate(readings=readings,)
matrices = format_result([adjusted_matrix.matrix for adjusted_matrix in result])
expected_matrices = get_excpected_matrices("BOU", "inf_weekly")
@@ -321,9 +296,7 @@ def test_BOU201911202001_infinite_one_interval():
TranslateOrigins(memory=np.inf),
],
update_interval=None,
- ).calculate(
- readings=readings,
- )
+ ).calculate(readings=readings,)
matrices = format_result([adjusted_matrix.matrix for adjusted_matrix in result])
expected_matrices = get_excpected_matrices("BOU", "inf_one_interval")
@@ -356,9 +329,7 @@ def test_CMO2015_causal():
starttime=starttime,
endtime=endtime,
update_interval=update_interval,
- ).calculate(
- readings=readings,
- )
+ ).calculate(readings=readings,)
matrices = format_result([adjusted_matrix.matrix for adjusted_matrix in result])
expected_matrices = get_excpected_matrices("CMO", "short_causal")
@@ -392,9 +363,7 @@ def test_CMO2015_acausal():
endtime=endtime,
update_interval=update_interval,
acausal=True,
- ).calculate(
- readings=readings,
- )
+ ).calculate(readings=readings,)
matrices = format_result([adjusted_matrix.matrix for adjusted_matrix in result])
expected_matrices = get_excpected_matrices("CMO", "short_acausal")
@@ -432,9 +401,7 @@ def test_CMO2015_infinite_weekly():
],
update_interval=update_interval,
acausal=True,
- ).calculate(
- readings=readings,
- )
+ ).calculate(readings=readings,)
matrices = format_result([adjusted_matrix.matrix for adjusted_matrix in result])
expected_matrices = get_excpected_matrices("CMO", "inf_weekly")
@@ -468,9 +435,7 @@ def test_CMO2015_infinite_one_interval():
],
acausal=True,
update_interval=None,
- ).calculate(
- readings=readings,
- )
+ ).calculate(readings=readings,)
matrices = format_result([adjusted_matrix.matrix for adjusted_matrix in result])
expected_matrices = get_excpected_matrices("CMO", "inf_one_interval")
--
GitLab