diff --git a/geomagio/adjusted/AdjustedMatrix.py b/geomagio/adjusted/AdjustedMatrix.py index 5913d2a5af9a013584efd6ffa6264287f69b45b4..c8b8dcc649c14a3008389c0438b49edf5400c0eb 100644 --- a/geomagio/adjusted/AdjustedMatrix.py +++ b/geomagio/adjusted/AdjustedMatrix.py @@ -1,8 +1,13 @@ from obspy import UTCDateTime from pydantic import BaseModel -from typing import Optional, Any, List +from typing import ( + Any, + List, + Optional, +) from .. import pydantic_utcdatetime +from .Metric import Metric class AdjustedMatrix(BaseModel): @@ -20,6 +25,6 @@ class AdjustedMatrix(BaseModel): matrix: Any pier_correction: float - metrics: List[float] = [0.0, 0.0] + metrics: Optional[List[Metric]] = None starttime: Optional[UTCDateTime] = None endtime: Optional[UTCDateTime] = None diff --git a/geomagio/adjusted/Affine.py b/geomagio/adjusted/Affine.py index ffd395afdc4b2eb1a62207f89de65fbbb2880d15..7c53cf5e0100a6e8c62dd06a8e65b9cbad9ccd88 100644 --- a/geomagio/adjusted/Affine.py +++ b/geomagio/adjusted/Affine.py @@ -5,7 +5,9 @@ from pydantic import BaseModel from typing import List, Optional, Tuple from .AdjustedMatrix import AdjustedMatrix +from .. import ChannelConverter from .. import pydantic_utcdatetime +from .Metric import Metric from ..residual import Reading from .Transform import Transform, TranslateOrigins, RotationTranslationXY @@ -273,13 +275,44 @@ class Affine(BaseModel): # compose affine transform matrices using reverse ordered matrices M_composed = reduce(np.dot, np.flipud(Ms)) + absolutes = np.vstack((absolutes, np.ones_like(absolutes[0]))) + ordinates = np.vstack((ordinates, np.ones_like(ordinates[0]))) pier_correction = np.average( [reading.pier_correction for reading in readings], weights=weights ) + std, mae, mae_df, std_df = self.compute_metrics( + absolutes=absolutes, ordinates=ordinates, matrix=M_composed + ) return AdjustedMatrix( - matrix=M_composed, pier_correction=pier_correction, metrics=metrics + matrix=M_composed, + pier_correction=pier_correction, + metrics=[ + Metric(element="X", mae=mae[0], std=std[0]), + Metric(element="Y", mae=mae[1], std=std[1]), + Metric(element="Z", mae=mae[2], std=std[2]), + Metric(element="dF", mae=mae_df, std=std_df), + ], + ) + + def compute_metrics( + self, absolutes: List[float], ordinates: List[float], matrix: List[float] + ) -> Tuple[List[float], List[float], float, float]: + # expected values are absolutes + predicted = matrix @ ordinates + # mean absolute erros and standard deviations ignore the 4th row comprison, which is trivial + std = np.nanstd(predicted - absolutes, axis=1)[0:3] + mae = abs(np.nanmean(predicted - absolutes, axis=1))[0:3] + expected_f = ChannelConverter.get_computed_f_using_squares( + absolutes[0], absolutes[1], absolutes[2] + ) + predicted_f = ChannelConverter.get_computed_f_using_squares( + predicted[0], predicted[1], predicted[2] ) + df = ChannelConverter.get_deltaf(expected_f, predicted_f) + std_df = abs(np.nanstd(df)) + mae_df = abs(np.nanmean(df)) + return list(std), list(mae), std_df, mae_df def get_weights( self, diff --git a/geomagio/adjusted/Metric.py b/geomagio/adjusted/Metric.py new file mode 100644 index 0000000000000000000000000000000000000000..19ba74adab14143b05597e32a9cc095ab98666a2 --- /dev/null +++ b/geomagio/adjusted/Metric.py @@ -0,0 +1,7 @@ +from pydantic import BaseModel + + +class Metric(BaseModel): + element: str + mae: float + std: float diff --git a/test/adjusted_test/adjusted_test.py b/test/adjusted_test/adjusted_test.py index eaebd44714d93ffd41e12f1aa4fbe710cc98d260..26047ccccb1d6479dd50e7f2d841827a5f8e63bb 100644 --- a/test/adjusted_test/adjusted_test.py +++ b/test/adjusted_test/adjusted_test.py @@ -244,7 +244,6 @@ 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) @@ -276,7 +275,6 @@ 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) @@ -312,7 +310,6 @@ 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) @@ -342,7 +339,6 @@ 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) @@ -379,7 +375,6 @@ 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) @@ -417,7 +412,6 @@ 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) @@ -459,7 +453,6 @@ 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) @@ -497,6 +490,5 @@ 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)