diff --git a/geomagio/adjusted/transform/LeastSq.py b/geomagio/adjusted/transform/LeastSq.py
index 7f0148737b20060104ae537451f39e3bdc42796c..45c6cb8bc8260b6b4d4b0939ab8c7a33e829c479 100644
--- a/geomagio/adjusted/transform/LeastSq.py
+++ b/geomagio/adjusted/transform/LeastSq.py
@@ -4,10 +4,24 @@ from typing import List, Optional, Tuple
from .Transform import Transform
-# TODO: GET_STACKED_ORDINATES SO METHOD CAN BE SHARED
+
class LeastSq(Transform):
"""Intance of Transform. Applies least squares to generate matrices"""
+ def get_stacked_values(self, absolutes, ordinates):
+ # LHS, or dependent variables
+ # [A[0,0], A[1,0], A[2,0], A[0,1], A[1,1], A[2,1], ...]
+ abs_stacked = self.get_stacked_absolutes(absolutes)
+ # RHS, or independent variables
+ # [
+ # [o[0,0], 0, 0, o[0,1], 0, 0, ...],
+ # [0, o[1,0], 0, 0, o[1,1], 0, ...],
+ # [0, 0, o[2,0], 0, 0, o[2,1], ...],
+ # ...
+ # ]
+ ord_stacked = self.get_stacked_ordinates(ordinates)
+ return abs_stacked, ord_stacked
+
def get_stacked_absolutes(self, absolutes):
"""Formats absolutes for least squares method
@@ -21,6 +35,25 @@ class LeastSq(Transform):
"""
return np.vstack([absolutes[0], absolutes[1], absolutes[2]]).T.ravel()
+ def get_stacked_ordinates(self, ordinates):
+ # (reduces degrees of freedom by 4:
+ # - 4 for the last row of zeros and a one)
+ ord_stacked = np.zeros((12, len(ordinates[0]) * 3))
+ ord_stacked[0, 0::3] = ordinates[0]
+ ord_stacked[1, 0::3] = ordinates[1]
+ ord_stacked[2, 0::3] = ordinates[2]
+ ord_stacked[3, 0::3] = 1.0
+ ord_stacked[4, 1::3] = ordinates[0]
+ ord_stacked[5, 1::3] = ordinates[1]
+ ord_stacked[6, 1::3] = ordinates[2]
+ ord_stacked[7, 1::3] = 1.0
+ ord_stacked[8, 2::3] = ordinates[0]
+ ord_stacked[9, 2::3] = ordinates[1]
+ ord_stacked[10, 2::3] = ordinates[2]
+ ord_stacked[11, 2::3] = 1.0
+
+ return ord_stacked
+
def get_weighted_values(
self,
values: Tuple[List[float], List[float], List[float]],
@@ -50,46 +83,22 @@ class LeastSq(Transform):
weights: List[float],
) -> np.array:
"""Calculates affine with no constraints using least squares."""
- # LHS, or dependent variables
- #
- # [A[0,0], A[1,0], A[2,0], A[0,1], A[1,1], A[2,1], ...]
- abs_stacked = self.get_stacked_absolutes(absolutes)
- # RHS, or independent variables
- # (reduces degrees of freedom by 4:
- # - 4 for the last row of zeros and a one)
- # [
- # [o[0,0], 0, 0, o[0,1], 0, 0, ...],
- # [0, o[1,0], 0, 0, o[1,1], 0, ...],
- # [0, 0, o[2,0], 0, 0, o[2,1], ...],
- # ...
- # ]
- ord_stacked = np.zeros((12, len(ordinates[0]) * 3))
- ord_stacked[0, 0::3] = ordinates[0]
- ord_stacked[1, 0::3] = ordinates[1]
- ord_stacked[2, 0::3] = ordinates[2]
- ord_stacked[3, 0::3] = 1.0
- ord_stacked[4, 1::3] = ordinates[0]
- ord_stacked[5, 1::3] = ordinates[1]
- ord_stacked[6, 1::3] = ordinates[2]
- ord_stacked[7, 1::3] = 1.0
- ord_stacked[8, 2::3] = ordinates[0]
- ord_stacked[9, 2::3] = ordinates[1]
- ord_stacked[10, 2::3] = ordinates[2]
- ord_stacked[11, 2::3] = 1.0
-
+ abs_stacked, ord_stacked = self.get_stacked_values(absolutes, ordinates)
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)
+ matrix, res, rank, sigma = spl.lstsq(ord_stacked.T, abs_stacked.T)
if rank < 3:
print("Poorly conditioned or singular matrix, returning NaNs")
return np.nan * np.ones((4, 4))
+ return self.format_matrix(matrix)
+
+ def format_matrix(self, matrix):
return np.array(
[
- [M_r[0], M_r[1], M_r[2], M_r[3]],
- [M_r[4], M_r[5], M_r[6], M_r[7]],
- [M_r[8], M_r[9], M_r[10], M_r[11]],
+ [matrix[0], matrix[1], matrix[2], matrix[3]],
+ [matrix[4], matrix[5], matrix[6], matrix[7]],
+ [matrix[8], matrix[9], matrix[10], matrix[11]],
[0.0, 0.0, 0.0, 1.0],
]
)
diff --git a/geomagio/adjusted/transform/Rescale3D.py b/geomagio/adjusted/transform/Rescale3D.py
index 5fdf02ae126605fec05c134998c31b517c91ea05..c55ca8f1e45003296c60b85a59c106cb553f2261 100644
--- a/geomagio/adjusted/transform/Rescale3D.py
+++ b/geomagio/adjusted/transform/Rescale3D.py
@@ -7,14 +7,7 @@ from .LeastSq import LeastSq
class Rescale3D(LeastSq):
"""Calculates affine using using least squares, constrained to re-scale each axis"""
- def calculate(
- self,
- ordinates: Tuple[List[float], List[float], List[float]],
- absolutes: Tuple[List[float], List[float], List[float]],
- weights: List[float],
- ) -> np.array:
- abs_stacked = self.get_stacked_absolutes(absolutes=absolutes)
- # RHS, or independent variables
+ def get_stacked_ordinates(self, ordinates):
# (reduces degrees of freedom by 13:
# - 2 for making x independent of y,z;
# - 2 for making y,z independent of x;
@@ -26,19 +19,12 @@ class Rescale3D(LeastSq):
ord_stacked[0, 0::3] = ordinates[0]
ord_stacked[1, 1::3] = ordinates[1]
ord_stacked[2, 2::3] = ordinates[2]
+ return ord_stacked
- 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)
- if rank < 3:
- print("Poorly conditioned or singular matrix, returning NaNs")
- return np.nan * np.ones((4, 4))
-
+ def format_matrix(self, matrix):
return [
- [M_r[0], 0.0, 0.0, 0.0],
- [0.0, M_r[1], 0.0, 0.0],
- [0.0, 0.0, M_r[2], 0.0],
+ [matrix[0], 0.0, 0.0, 0.0],
+ [0.0, matrix[1], 0.0, 0.0],
+ [0.0, 0.0, matrix[2], 0.0],
[0.0, 0.0, 0.0, 1.0],
]
diff --git a/geomagio/adjusted/transform/ShearYZ.py b/geomagio/adjusted/transform/ShearYZ.py
index b413ff58362a6dac9426da8b5c058cd931ca6b9b..dc2e43c213bdb30c976655a7503b79d98157d14a 100644
--- a/geomagio/adjusted/transform/ShearYZ.py
+++ b/geomagio/adjusted/transform/ShearYZ.py
@@ -8,13 +8,7 @@ from .LeastSq import LeastSq
class ShearYZ(LeastSq):
"""Calculates affine using least squares, constrained to shear y and z, but not x."""
- def calculate(
- self,
- ordinates: Tuple[List[float], List[float], List[float]],
- absolutes: Tuple[List[float], List[float], List[float]],
- weights: List[float],
- ) -> np.array:
- # RHS, or independent variables
+ def get_stacked_ordinates(self, ordinates):
# (reduces degrees of freedom by 13:
# - 2 for making x independent of y,z;
# - 1 for making y independent of z;
@@ -27,18 +21,12 @@ class ShearYZ(LeastSq):
ord_stacked[2, 0::3] = ordinates[0]
ord_stacked[2, 1::3] = ordinates[1]
ord_stacked[2, 2::3] = 1.0
- 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:
- print("Poorly conditioned or singular matrix, returning NaNs")
- return np.nan * np.ones((4, 4))
+ return ord_stacked
+ def format_matrix(self, matrix):
return [
[1.0, 0.0, 0.0, 0.0],
- [M_r[0], 1.0, 0.0, 0.0],
- [M_r[1], M_r[2], 1.0, 0.0],
+ [matrix[0], 1.0, 0.0, 0.0],
+ [matrix[1], matrix[2], 1.0, 0.0],
[0.0, 0.0, 0.0, 1.0],
]
diff --git a/geomagio/adjusted/transform/TranslateOrigins.py b/geomagio/adjusted/transform/TranslateOrigins.py
index 7b416fa7ecb5e6f0cf3f1744c243336368223f92..69c21915fb007c0bcea79e7c0c6dcb958752a40a 100644
--- a/geomagio/adjusted/transform/TranslateOrigins.py
+++ b/geomagio/adjusted/transform/TranslateOrigins.py
@@ -8,6 +8,17 @@ from .LeastSq import LeastSq
class TranslateOrigins(LeastSq):
"""Calculates affine using using least squares, constrained to tanslate origins"""
+ def get_stacked_values(self, absolutes, ordinates):
+ # LHS, or dependent variables
+ abs_stacked = self.get_stacked_absolutes(absolutes)
+ # subtract ords from abs to force simple translation
+ abs_stacked[0::3] = absolutes[0] - ordinates[0]
+ abs_stacked[1::3] = absolutes[1] - ordinates[1]
+ abs_stacked[2::3] = absolutes[2] - ordinates[2]
+ # RHS, or independent variables
+ ord_stacked = self.get_stacked_ordinates(ordinates)
+ return abs_stacked, ord_stacked
+
def get_weighted_values(
self,
values: Tuple[List[float], List[float], List[float]],
@@ -22,38 +33,17 @@ class TranslateOrigins(LeastSq):
weights = 1
return values * weights
- def calculate(
- self,
- ordinates: Tuple[List[float], List[float], List[float]],
- absolutes: Tuple[List[float], List[float], List[float]],
- weights: List[float],
- ) -> np.array:
-
- abs_stacked = self.get_stacked_absolutes(absolutes)
+ def get_stacked_ordinates(self, ordinates):
ord_stacked = np.zeros((3, len(ordinates[0]) * 3))
-
ord_stacked[0, 0::3] = 1.0
ord_stacked[1, 1::3] = 1.0
ord_stacked[2, 2::3] = 1.0
+ return ord_stacked
- # subtract ords from abs to force simple translation
- abs_stacked[0::3] = absolutes[0] - ordinates[0]
- abs_stacked[1::3] = absolutes[1] - ordinates[1]
- abs_stacked[2::3] = absolutes[2] - ordinates[2]
-
- # apply weights
- 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:
- print("Poorly conditioned or singular matrix, returning NaNs")
- return np.nan * np.ones((4, 4))
-
+ def format_matrix(self, matrix):
return [
- [1.0, 0.0, 0.0, M_r[0]],
- [0.0, 1.0, 0.0, M_r[1]],
- [0.0, 0.0, 1.0, M_r[2]],
+ [1.0, 0.0, 0.0, matrix[0]],
+ [0.0, 1.0, 0.0, matrix[1]],
+ [0.0, 0.0, 1.0, matrix[2]],
[0.0, 0.0, 0.0, 1.0],
]
diff --git a/geomagio/adjusted/transform/ZRotationHScale.py b/geomagio/adjusted/transform/ZRotationHScale.py
index bc08bfc7d744a0f28c746717d71e38be1795802c..a5536bce8202b0a25f40cc6622bef37ec0c181cf 100644
--- a/geomagio/adjusted/transform/ZRotationHScale.py
+++ b/geomagio/adjusted/transform/ZRotationHScale.py
@@ -9,16 +9,7 @@ class ZRotationHscale(LeastSq):
"""Calculates affine using least squares, constrained to rotate about the Z axis
and apply uniform horizontal scaling."""
- def calculate(
- self,
- ordinates: Tuple[List[float], List[float], List[float]],
- absolutes: Tuple[List[float], List[float], List[float]],
- weights: List[float],
- ) -> np.array:
- # LHS, or dependent variables
- #
- abs_stacked = self.get_stacked_absolutes(absolutes)
- # RHS, or independent variables
+ def get_stacked_ordinates(self, ordinates):
# (reduces degrees of freedom by 10:
# - 2 for making x,y independent of z;
# - 2 for making z independent of x,y
@@ -33,19 +24,12 @@ class ZRotationHscale(LeastSq):
ord_stacked[3, 1::3] = 1.0
ord_stacked[4, 2::3] = ordinates[2]
ord_stacked[5, 2::3] = 1.0
+ return ord_stacked
- 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)
- if rank < 3:
- print("Poorly conditioned or singular matrix, returning NaNs")
- return np.nan * np.ones((4, 4))
-
+ def format_matrix(self, matrix):
return [
- [M_r[0], M_r[1], 0.0, M_r[2]],
- [-M_r[1], M_r[0], 0.0, M_r[3]],
- [0.0, 0.0, M_r[4], M_r[5]],
+ [matrix[0], matrix[1], 0.0, matrix[2]],
+ [-matrix[1], matrix[0], 0.0, matrix[3]],
+ [0.0, 0.0, matrix[4], matrix[5]],
[0.0, 0.0, 0.0, 1.0],
]
diff --git a/geomagio/adjusted/transform/ZRotationHScaleZBaseline.py b/geomagio/adjusted/transform/ZRotationHScaleZBaseline.py
index 0be5b249f609dc467516f2b6cdfc374fa1ada8dc..569786afcbf5be6f58704d2f9ab50a80158b2c13 100644
--- a/geomagio/adjusted/transform/ZRotationHScaleZBaseline.py
+++ b/geomagio/adjusted/transform/ZRotationHScaleZBaseline.py
@@ -9,13 +9,16 @@ class ZRotationHscaleZbaseline(LeastSq):
"""Calculates affine using least squares, constrained to rotate about the Z axis,
apply a uniform horizontal scaling, and apply a baseline shift for the Z axis."""
- def calculate(
- self,
- ordinates: Tuple[List[float], List[float], List[float]],
- absolutes: Tuple[List[float], List[float], List[float]],
- weights: List[float],
- ) -> np.array:
+ def get_stacked_values(self, absolutes, ordinates):
+ # LHS, or dependent variables
+ abs_stacked = self.get_stacked_absolutes(absolutes)
+ # subtract z_o from z_a to force simple z translation
+ abs_stacked[2::3] = absolutes[2] - ordinates[2]
# RHS, or independent variables
+ ord_stacked = self.get_stacked_ordinates(ordinates)
+ return abs_stacked, ord_stacked
+
+ def get_stacked_ordinates(self, ordinates):
# (reduces degrees of freedom by 13:
# - 2 for making x,y independent of z;
# - 2 for making z independent of x,y;
@@ -29,23 +32,12 @@ class ZRotationHscaleZbaseline(LeastSq):
ord_stacked[1, 0::3] = ordinates[1]
ord_stacked[1, 1::3] = -ordinates[0]
ord_stacked[2, 2::3] = 1.0
+ return ord_stacked
- # subtract z_o from z_a to force simple z translation
- abs_stacked = self.get_stacked_absolutes(absolutes)
- abs_stacked[2::3] = absolutes[2] - 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)
- if rank < 3:
- print("Poorly conditioned or singular matrix, returning NaNs")
- return np.nan * np.ones((4, 4))
-
+ def format_matrix(self, matrix):
return [
- [M_r[0], M_r[1], 0.0, 0.0],
- [-M_r[1], M_r[0], 0.0, 0.0],
- [0.0, 0.0, 1.0, M_r[2]],
+ [matrix[0], matrix[1], 0.0, 0.0],
+ [-matrix[1], matrix[0], 0.0, 0.0],
+ [0.0, 0.0, 1.0, matrix[2]],
[0.0, 0.0, 0.0, 1.0],
]
diff --git a/geomagio/adjusted/transform/ZRotationShear.py b/geomagio/adjusted/transform/ZRotationShear.py
index 322eb8dbf226c9bf789be0d0cd6334c4005bbf24..693cbb4152b5081cdc8806349dd58e6e02c6c286 100644
--- a/geomagio/adjusted/transform/ZRotationShear.py
+++ b/geomagio/adjusted/transform/ZRotationShear.py
@@ -6,17 +6,7 @@ from .LeastSq import LeastSq
class ZRotationShear(LeastSq):
- def calculate(
- self,
- ordinates: Tuple[List[float], List[float], List[float]],
- absolutes: Tuple[List[float], List[float], List[float]],
- weights: List[float],
- ) -> np.array:
- """Calculates affine using least squares, constrained to rotate about the Z axis."""
- # LHS, or dependent variables
- #
- abs_stacked = self.get_stacked_absolutes(absolutes)
- # RHS, or independent variables
+ def get_stacked_ordinates(self, ordinates):
# (reduces degrees of freedom by 8:
# - 2 for making x,y independent of z;
# - 2 for making z independent of x,y
@@ -30,19 +20,12 @@ class ZRotationShear(LeastSq):
ord_stacked[5, 1::3] = 1.0
ord_stacked[6, 2::3] = ordinates[2]
ord_stacked[7, 2::3] = 1.0
+ return ord_stacked
- 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)
- if rank < 3:
- print("Poorly conditioned or singular matrix, returning NaNs")
- return np.nan * np.ones((4, 4))
-
+ def format_matrix(self, matrix):
return [
- [M_r[0], M_r[1], 0.0, M_r[2]],
- [M_r[3], M_r[4], 0.0, M_r[5]],
- [0.0, 0.0, M_r[6], M_r[7]],
+ [matrix[0], matrix[1], 0.0, matrix[2]],
+ [matrix[3], matrix[4], 0.0, matrix[5]],
+ [0.0, 0.0, matrix[6], matrix[7]],
[0.0, 0.0, 0.0, 1.0],
]