diff --git a/geomagio/adjusted/AdjustedMatrix.py b/geomagio/adjusted/AdjustedMatrix.py index 38714db9c0b1e0a22deac425b321ba44314edd8c..4d163723659976615b707e91e765ecb7c5296671 100644 --- a/geomagio/adjusted/AdjustedMatrix.py +++ b/geomagio/adjusted/AdjustedMatrix.py @@ -3,9 +3,9 @@ from obspy import Stream, UTCDateTime from pydantic import BaseModel from typing import Any, List, Optional +from ..residual.Reading import Reading, get_absolutes_xyz, get_ordinates from .. import ChannelConverter from .. import pydantic_utcdatetime -from ..residual.Reading import Reading, get_absolutes_xyz, get_ordinates from .Metric import Metric @@ -51,11 +51,12 @@ class AdjustedMatrix(BaseModel): return adjusted def get_metrics(self, readings: List[Reading]) -> List[Metric]: - """Computes mean absolute error and standard deviation for X, Y, Z, and dF between expected and predicted values. + """Computes mean absolute error and standard deviation between expected and predicted values + Metrics are computed for X, Y, Z, and dF values Attributes ---------- - readings: list of Readings + readings: list of valid readings matrix: composed matrix Outputs diff --git a/geomagio/adjusted/Affine.py b/geomagio/adjusted/Affine.py index 055461d5728ab285abce503df445a30edd123948..1df35fca3d114160ca2c08e1e102d9385900c041 100644 --- a/geomagio/adjusted/Affine.py +++ b/geomagio/adjusted/Affine.py @@ -4,12 +4,12 @@ from obspy import UTCDateTime from pydantic import BaseModel from typing import List, Optional, Tuple -from .. import pydantic_utcdatetime from ..residual.Reading import ( Reading, get_absolutes_xyz, get_ordinates, ) +from .. import pydantic_utcdatetime from .AdjustedMatrix import AdjustedMatrix from .transform import RotationTranslationXY, TranslateOrigins, Transform @@ -22,9 +22,8 @@ class Affine(BaseModel): observatory: 3-letter observatory code starttime: beginning time for matrix creation endtime: end time for matrix creation - acausal: when True, utilizes readings from after set endtime - update_interval: window of time a matrix is representative of - transforms: list of methods for matrix calculations + update_interval: window of time(in seconds) a matrix is representative of + transforms: methods for matrix calculations """ observatory: str = None @@ -46,11 +45,12 @@ class Affine(BaseModel): Attributes ---------- - readings: list of readings containing absolutes + readings: readings containing absolutes + epochs: optional time markers for unreliable observations Outputs ------- - Ms: list of AdjustedMatrix objects created from calculations + Ms: AdjustedMatrix objects created from calculations """ # default set to create one matrix between starttime and endtime update_interval = self.update_interval or (self.endtime - self.starttime) @@ -141,15 +141,13 @@ def get_epochs( Attributes ---------- - epoch_start: float value signifying start of last valid interval - epoch_end: float value signifying end of last valid interval - epochs: list of floats signifying bad data times + epochs: bad data times time: current time epoch is being evaluated at Outputs ------- - epoch_start: float value signifying start of current valid interval - epoch_end: float value signifying end of current valid interval + epoch_start: start of current valid interval + epoch_end: end of current valid interval """ epoch_start = None epoch_end = None diff --git a/geomagio/adjusted/transform/Rescale3D.py b/geomagio/adjusted/transform/Rescale3D.py index 40154b5ed9b7d3ed7f8af705158e492b2fe62558..aa1c58b2f9bfeae68e716c6eb5e5c0f4dfa9a819 100644 --- a/geomagio/adjusted/transform/Rescale3D.py +++ b/geomagio/adjusted/transform/Rescale3D.py @@ -1,5 +1,6 @@ import numpy as np from typing import List, Optional, Tuple + from .LeastSq import LeastSq diff --git a/geomagio/adjusted/transform/Transform.py b/geomagio/adjusted/transform/Transform.py index 2263290ac43fbcd5cd381c7148344280d2597738..b3b7a391ce306bf9f2d6b8530e01535d85e93531 100644 --- a/geomagio/adjusted/transform/Transform.py +++ b/geomagio/adjusted/transform/Transform.py @@ -40,8 +40,12 @@ class Transform(BaseModel): Calculate time-dependent weights according to exponential decay. Inputs: - times: array of times, or any time-like index whose relative values represent spacing between events + ------- + readings: list of valid readings + time: time weights are calculated for + Output: + ------- weights: array of vector distances/metrics """ @@ -56,7 +60,7 @@ class Transform(BaseModel): # if memory is actually infinite, return equal weights if np.isinf(self.memory): - return filter_iqrs(multiseries=baselines, weights=np.ones(times.shape)) + weights = np.ones(times.shape) # initialize weights weights = np.zeros(times.shape) diff --git a/test/adjusted_test/adjusted_test.py b/test/adjusted_test/adjusted_test.py index 6e4ef63685c77635dbbe5f0dc2394c721f06e542..04dc921e5b1be69eca14d040af12691b8521fde4 100644 --- a/test/adjusted_test/adjusted_test.py +++ b/test/adjusted_test/adjusted_test.py @@ -10,8 +10,8 @@ from geomagio.adjusted.transform import ( LeastSq, QRFactorization, Rescale3D, - SVD, RotationTranslationXY, + SVD, ShearYZ, TranslateOrigins, ZRotationHscale, @@ -24,152 +24,6 @@ from geomagio.residual import WebAbsolutesFactory from test.residual_test.residual_test import get_spreadsheet_directory_readings -def get_sythetic_variables(): - with open("etc/adjusted/synthetic.json") as file: - data = json.load(file) - 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"]]) - weights = np.arange(0, len(ordinates[0])) - return ordinates, absolutes, weights - - -def get_expected_synthetic_result(key): - with open("etc/adjusted/synthetic.json") as file: - expected = json.load(file) - return expected["results"][key] - - -def test_LeastSq_synthetic(): - ordinates, absolutes, weights = get_sythetic_variables() - assert_array_almost_equal( - LeastSq().calculate( - ordinates=ordinates, - absolutes=absolutes, - weights=weights, - ), - get_expected_synthetic_result("LeastSq"), - decimal=3, - ) - - -def test_QRFactorization_synthetic(): - ordinates, absolutes, weights = get_sythetic_variables() - assert_array_almost_equal( - QRFactorization().calculate( - ordinates=ordinates, - absolutes=absolutes, - weights=weights, - ), - get_expected_synthetic_result("QRFactorization"), - decimal=3, - ) - - -def test_Rescale3D_synthetic(): - ordinates, absolutes, weights = get_sythetic_variables() - assert_array_almost_equal( - Rescale3D().calculate( - ordinates=ordinates, - absolutes=absolutes, - weights=weights, - ), - get_expected_synthetic_result("Rescale3D"), - decimal=3, - ) - - -def test_SVD_synthetic(): - ordinates, absolutes, weights = get_sythetic_variables() - assert_array_almost_equal( - SVD().calculate( - ordinates=ordinates, - absolutes=absolutes, - weights=weights, - ), - get_expected_synthetic_result("SVD"), - decimal=3, - ) - - -def test_RotationTranslationXY_synthetic(): - ordinates, absolutes, weights = get_sythetic_variables() - assert_array_almost_equal( - RotationTranslationXY().calculate( - ordinates=ordinates, - absolutes=absolutes, - weights=weights, - ), - get_expected_synthetic_result("RotationTranslationXY"), - decimal=3, - ) - - -def test_ShearYZ_synthetic(): - ordinates, absolutes, weights = get_sythetic_variables() - assert_array_almost_equal( - ShearYZ().calculate( - ordinates=ordinates, - absolutes=absolutes, - weights=weights, - ), - get_expected_synthetic_result("ShearYZ"), - decimal=3, - ) - - -def test_TranslateOrigins_synthetic(): - ordinates, absolutes, weights = get_sythetic_variables() - assert_array_almost_equal( - TranslateOrigins().calculate( - ordinates=ordinates, - absolutes=absolutes, - weights=weights, - ), - get_expected_synthetic_result("TranslateOrigins"), - decimal=3, - ) - - -def test_ZRotationHscale_synthetic(): - ordinates, absolutes, weights = get_sythetic_variables() - assert_array_almost_equal( - ZRotationHscale().calculate( - ordinates=ordinates, - absolutes=absolutes, - weights=weights, - ), - get_expected_synthetic_result("ZRotationHscale"), - decimal=3, - ) - - -def test_ZRotationHscaleZbaseline_synthetic(): - ordinates, absolutes, weights = get_sythetic_variables() - assert_array_almost_equal( - ZRotationHscaleZbaseline().calculate( - ordinates=ordinates, - absolutes=absolutes, - weights=weights, - ), - get_expected_synthetic_result("ZRotationHscaleZbaseline"), - decimal=3, - ) - - -def test_ZRotationShear_synthetic(): - ordinates, absolutes, weights = get_sythetic_variables() - assert_array_almost_equal( - ZRotationShear().calculate( - ordinates=ordinates, - absolutes=absolutes, - weights=weights, - ), - get_expected_synthetic_result("ZRotationShear"), - decimal=3, - ) - - def format_result(result) -> dict: Ms = [] for M in result: @@ -186,6 +40,12 @@ def get_excpected_matrices(observatory, key): return expected[key] +def get_expected_synthetic_result(key): + with open("etc/adjusted/synthetic.json") as file: + expected = json.load(file) + return expected["results"][key] + + def get_readings_BOU201911202001(): readings = WebAbsolutesFactory().get_readings( observatory="BOU", @@ -195,7 +55,44 @@ def get_readings_BOU201911202001(): return readings -def test_BOU201911202001_short_causal(): +def get_sythetic_variables(): + with open("etc/adjusted/synthetic.json") as file: + data = json.load(file) + 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"]]) + weights = np.arange(0, len(ordinates[0])) + return ordinates, absolutes, weights + + +def test_BOU201911202001_infinite_one_interval(): + readings = get_readings_BOU201911202001() + result = Affine( + observatory="BOU", + starttime=UTCDateTime("2019-11-01T00:00:00Z"), + endtime=UTCDateTime("2020-01-31T23:59:00Z"), + transforms=[ + RotationTranslationXY(memory=np.inf, acausal=True), + TranslateOrigins(memory=np.inf, acausal=True), + ], + update_interval=None, + ).calculate( + readings=readings, + ) + + matrices = format_result([adjusted_matrix.matrix for adjusted_matrix in result]) + expected_matrices = get_excpected_matrices("BOU", "inf_one_interval") + for i in range(len(matrices)): + assert_array_almost_equal( + matrices[i], + expected_matrices[i], + decimal=3, + err_msg=f"Matrix {i} not equal", + ) + assert_equal(len(matrices), 1) + + +def test_BOU201911202001_infinite_weekly(): readings = get_readings_BOU201911202001() starttime = UTCDateTime("2019-11-01T00:00:00Z") @@ -209,13 +106,15 @@ def test_BOU201911202001_short_causal(): endtime=endtime, update_interval=update_interval, transforms=[ - RotationTranslationXY(memory=(86400 * 100), acausal=False), - TranslateOrigins(memory=(86400 * 10), acausal=False), + RotationTranslationXY(memory=np.inf, acausal=True), + TranslateOrigins(memory=np.inf, 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_causal") + expected_matrices = get_excpected_matrices("BOU", "inf_weekly") for i in range(len(matrices)): assert_array_almost_equal( matrices[i], @@ -226,6 +125,22 @@ def test_BOU201911202001_short_causal(): assert_equal(len(matrices), ((endtime - starttime) // update_interval) + 1) +def test_BOU201911202001_invalid_readings(): + readings = [] + starttime = UTCDateTime("2019-11-01T00:00:00Z") + result = Affine( + observatory="BOU", + starttime=starttime, + endtime=UTCDateTime("2020-01-31T23:59:00Z"), + transforms=[ + RotationTranslationXY(memory=np.inf, acausal=True), + TranslateOrigins(memory=np.inf, acausal=True), + ], + update_interval=None, + ).calculate(readings=readings,)[0] + assert result == AdjustedMatrix(time=starttime) + + def test_BOU201911202001_short_acausal(): readings = get_readings_BOU201911202001() @@ -259,7 +174,7 @@ def test_BOU201911202001_short_acausal(): assert_equal(len(matrices), ((endtime - starttime) // update_interval) + 1) -def test_BOU201911202001_infinite_weekly(): +def test_BOU201911202001_short_causal(): readings = get_readings_BOU201911202001() starttime = UTCDateTime("2019-11-01T00:00:00Z") @@ -273,15 +188,13 @@ def test_BOU201911202001_infinite_weekly(): endtime=endtime, update_interval=update_interval, transforms=[ - RotationTranslationXY(memory=np.inf, acausal=True), - TranslateOrigins(memory=np.inf, acausal=True), + RotationTranslationXY(memory=(86400 * 100), acausal=False), + TranslateOrigins(memory=(86400 * 10), acausal=False), ], - ).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") + expected_matrices = get_excpected_matrices("BOU", "short_causal") for i in range(len(matrices)): assert_array_almost_equal( matrices[i], @@ -292,23 +205,32 @@ def test_BOU201911202001_infinite_weekly(): assert_equal(len(matrices), ((endtime - starttime) // update_interval) + 1) -def test_BOU201911202001_infinite_one_interval(): - readings = get_readings_BOU201911202001() +def test_CMO2015_infinite_one_interval(): + readings = get_spreadsheet_directory_readings( + observatory="CMO", + starttime=UTCDateTime("2015-01-01T00:00:00Z"), + endtime=UTCDateTime("2015-12-31T23:59:00Z"), + path="etc/residual/Caldata/", + ) + + assert len(readings) == 146 + result = Affine( - observatory="BOU", - starttime=UTCDateTime("2019-11-01T00:00:00Z"), - endtime=UTCDateTime("2020-01-31T23:59:00Z"), + observatory="CMO", + starttime=UTCDateTime("2015-02-01T00:00:00Z"), + endtime=UTCDateTime("2015-11-27T23:59:00Z"), transforms=[ RotationTranslationXY(memory=np.inf, acausal=True), TranslateOrigins(memory=np.inf, acausal=True), ], + acausal=True, update_interval=None, ).calculate( readings=readings, ) matrices = format_result([adjusted_matrix.matrix for adjusted_matrix in result]) - expected_matrices = get_excpected_matrices("BOU", "inf_one_interval") + expected_matrices = get_excpected_matrices("CMO", "inf_one_interval") for i in range(len(matrices)): assert_array_almost_equal( matrices[i], @@ -316,26 +238,11 @@ def test_BOU201911202001_infinite_one_interval(): decimal=3, err_msg=f"Matrix {i} not equal", ) - assert_equal(len(matrices), 1) - -def test_BOU201911202001_invalid_readings(): - readings = [] - starttime = UTCDateTime("2019-11-01T00:00:00Z") - result = Affine( - observatory="BOU", - starttime=starttime, - endtime=UTCDateTime("2020-01-31T23:59:00Z"), - transforms=[ - RotationTranslationXY(memory=np.inf, acausal=True), - TranslateOrigins(memory=np.inf, acausal=True), - ], - update_interval=None, - ).calculate(readings=readings,)[0] - assert result == AdjustedMatrix(time=starttime) + assert_equal(len(matrices), 1) -def test_CMO2015_causal(): +def test_CMO2015_infinite_weekly(): readings = get_spreadsheet_directory_readings( observatory="CMO", starttime=UTCDateTime("2015-01-01T00:00:00Z"), @@ -353,17 +260,18 @@ def test_CMO2015_causal(): observatory="CMO", starttime=starttime, endtime=endtime, - update_interval=update_interval, transforms=[ - RotationTranslationXY(memory=(86400 * 100), acausal=False), - TranslateOrigins(memory=(86400 * 10), acausal=False), + RotationTranslationXY(memory=np.inf, acausal=True), + TranslateOrigins(memory=np.inf, acausal=True), ], + update_interval=update_interval, + acausal=True, ).calculate( readings=readings, ) matrices = format_result([adjusted_matrix.matrix for adjusted_matrix in result]) - expected_matrices = get_excpected_matrices("CMO", "short_causal") + expected_matrices = get_excpected_matrices("CMO", "inf_weekly") for i in range(len(matrices)): assert_array_almost_equal( matrices[i], @@ -374,7 +282,7 @@ def test_CMO2015_causal(): assert_equal(len(matrices), ((endtime - starttime) // update_interval) + 1) -def test_CMO2015_acausal(): +def test_CMO2015_short_acausal(): readings = get_spreadsheet_directory_readings( observatory="CMO", starttime=UTCDateTime("2015-01-01T00:00:00Z"), @@ -413,7 +321,7 @@ def test_CMO2015_acausal(): assert_equal(len(matrices), ((endtime - starttime) // update_interval) + 1) -def test_CMO2015_infinite_weekly(): +def test_CMO2015_short_causal(): readings = get_spreadsheet_directory_readings( observatory="CMO", starttime=UTCDateTime("2015-01-01T00:00:00Z"), @@ -431,54 +339,17 @@ def test_CMO2015_infinite_weekly(): observatory="CMO", starttime=starttime, endtime=endtime, - transforms=[ - RotationTranslationXY(memory=np.inf, acausal=True), - TranslateOrigins(memory=np.inf, acausal=True), - ], update_interval=update_interval, - acausal=True, - ).calculate( - readings=readings, - ) - - matrices = format_result([adjusted_matrix.matrix for adjusted_matrix in result]) - expected_matrices = get_excpected_matrices("CMO", "inf_weekly") - for i in range(len(matrices)): - assert_array_almost_equal( - matrices[i], - expected_matrices[i], - decimal=3, - err_msg=f"Matrix {i} not equal", - ) - assert_equal(len(matrices), ((endtime - starttime) // update_interval) + 1) - - -def test_CMO2015_infinite_one_interval(): - readings = get_spreadsheet_directory_readings( - observatory="CMO", - starttime=UTCDateTime("2015-01-01T00:00:00Z"), - endtime=UTCDateTime("2015-12-31T23:59:00Z"), - path="etc/residual/Caldata/", - ) - - assert len(readings) == 146 - - result = Affine( - observatory="CMO", - starttime=UTCDateTime("2015-02-01T00:00:00Z"), - endtime=UTCDateTime("2015-11-27T23:59:00Z"), transforms=[ - RotationTranslationXY(memory=np.inf, acausal=True), - TranslateOrigins(memory=np.inf, acausal=True), + RotationTranslationXY(memory=(86400 * 100), acausal=False), + TranslateOrigins(memory=(86400 * 10), acausal=False), ], - acausal=True, - update_interval=None, ).calculate( readings=readings, ) matrices = format_result([adjusted_matrix.matrix for adjusted_matrix in result]) - expected_matrices = get_excpected_matrices("CMO", "inf_one_interval") + expected_matrices = get_excpected_matrices("CMO", "short_causal") for i in range(len(matrices)): assert_array_almost_equal( matrices[i], @@ -486,8 +357,7 @@ def test_CMO2015_infinite_one_interval(): decimal=3, err_msg=f"Matrix {i} not equal", ) - - assert_equal(len(matrices), 1) + assert_equal(len(matrices), ((endtime - starttime) // update_interval) + 1) def test_get_epochs(): @@ -516,3 +386,133 @@ def test_get_epochs(): epochs=epochs, time=UTCDateTime("2020-01-07T00:00:00Z"), ) + + +def test_LeastSq_synthetic(): + ordinates, absolutes, weights = get_sythetic_variables() + assert_array_almost_equal( + LeastSq().calculate( + ordinates=ordinates, + absolutes=absolutes, + weights=weights, + ), + get_expected_synthetic_result("LeastSq"), + decimal=3, + ) + + +def test_QRFactorization_synthetic(): + ordinates, absolutes, weights = get_sythetic_variables() + assert_array_almost_equal( + QRFactorization().calculate( + ordinates=ordinates, + absolutes=absolutes, + weights=weights, + ), + get_expected_synthetic_result("QRFactorization"), + decimal=3, + ) + + +def test_Rescale3D_synthetic(): + ordinates, absolutes, weights = get_sythetic_variables() + assert_array_almost_equal( + Rescale3D().calculate( + ordinates=ordinates, + absolutes=absolutes, + weights=weights, + ), + get_expected_synthetic_result("Rescale3D"), + decimal=3, + ) + + +def test_RotationTranslationXY_synthetic(): + ordinates, absolutes, weights = get_sythetic_variables() + assert_array_almost_equal( + RotationTranslationXY().calculate( + ordinates=ordinates, + absolutes=absolutes, + weights=weights, + ), + get_expected_synthetic_result("RotationTranslationXY"), + decimal=3, + ) + + +def test_ShearYZ_synthetic(): + ordinates, absolutes, weights = get_sythetic_variables() + assert_array_almost_equal( + ShearYZ().calculate( + ordinates=ordinates, + absolutes=absolutes, + weights=weights, + ), + get_expected_synthetic_result("ShearYZ"), + decimal=3, + ) + + +def test_SVD_synthetic(): + ordinates, absolutes, weights = get_sythetic_variables() + assert_array_almost_equal( + SVD().calculate( + ordinates=ordinates, + absolutes=absolutes, + weights=weights, + ), + get_expected_synthetic_result("SVD"), + decimal=3, + ) + + +def test_TranslateOrigins_synthetic(): + ordinates, absolutes, weights = get_sythetic_variables() + assert_array_almost_equal( + TranslateOrigins().calculate( + ordinates=ordinates, + absolutes=absolutes, + weights=weights, + ), + get_expected_synthetic_result("TranslateOrigins"), + decimal=3, + ) + + +def test_ZRotationHscale_synthetic(): + ordinates, absolutes, weights = get_sythetic_variables() + assert_array_almost_equal( + ZRotationHscale().calculate( + ordinates=ordinates, + absolutes=absolutes, + weights=weights, + ), + get_expected_synthetic_result("ZRotationHscale"), + decimal=3, + ) + + +def test_ZRotationHscaleZbaseline_synthetic(): + ordinates, absolutes, weights = get_sythetic_variables() + assert_array_almost_equal( + ZRotationHscaleZbaseline().calculate( + ordinates=ordinates, + absolutes=absolutes, + weights=weights, + ), + get_expected_synthetic_result("ZRotationHscaleZbaseline"), + decimal=3, + ) + + +def test_ZRotationShear_synthetic(): + ordinates, absolutes, weights = get_sythetic_variables() + assert_array_almost_equal( + ZRotationShear().calculate( + ordinates=ordinates, + absolutes=absolutes, + weights=weights, + ), + get_expected_synthetic_result("ZRotationShear"), + decimal=3, + )