diff --git a/geomagio/adjusted/transform/LeastSq.py b/geomagio/adjusted/transform/LeastSq.py index f33eace1564b7bac087a7e27de70e37c932f2d2a..0caa4b1a03f47311354aff3a56434cdf11c762cb 100644 --- a/geomagio/adjusted/transform/LeastSq.py +++ b/geomagio/adjusted/transform/LeastSq.py @@ -1,6 +1,6 @@ import numpy as np import scipy.linalg as spl -from typing import List, Optional, Tuple +from typing import List, Optional, Tuple, Union from .Transform import Transform @@ -8,21 +8,49 @@ from .Transform import Transform class LeastSq(Transform): """Intance of Transform. Applies least squares to generate matrices""" - def get_stacked_values(self, absolutes, ordinates, weights=None): - # 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 calculate( + self, + ordinates: Tuple[List[float], List[float], List[float]], + absolutes: Tuple[List[float], List[float], List[float]], + weights: Optional[List[float]] = None, + ) -> np.array: + """Calculates matrix with least squares and accompanying methods + Defaults to least squares calculation with no constraints + """ + abs_stacked, ord_stacked = self.get_stacked_values( + absolutes, ordinates, weights + ) + 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 + matrix, res, rank, sigma = spl.lstsq(ord_stacked.T, abs_stacked.T) + if self.valid(rank): + return self.get_matrix(matrix, absolutes, ordinates, weights) + print("Poorly conditioned or singular matrix, returning NaNs") + return np.nan * np.ones((4, 4)) + + def get_matrix( + self, + matrix: List[List[float]], + absolutes: Optional[Tuple[List[float], List[float], List[float]]] = None, + ordinates: Optional[Tuple[List[float], List[float], List[float]]] = None, + weights: Optional[List[float]] = None, + ) -> np.array: + """Returns matrix formatted for no constraints + NOTE: absolutes, ordinates, and weights are only used by QRFactorization's child function + """ + return np.array( + [ + [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], + ] + ) - def get_stacked_absolutes(self, absolutes): + def get_stacked_absolutes( + self, absolutes: Tuple[List[float], List[float], List[float]] + ) -> List[float]: """Formats absolutes for least squares method Attributes @@ -35,7 +63,10 @@ class LeastSq(Transform): """ return np.vstack([absolutes[0], absolutes[1], absolutes[2]]).T.ravel() - def get_stacked_ordinates(self, ordinates): + def get_stacked_ordinates( + self, ordinates: Tuple[List[float], List[float], List[float]] + ) -> List[List[float]]: + """ Formats ordinates for least squares method """ # (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)) @@ -54,17 +85,34 @@ class LeastSq(Transform): return ord_stacked - def valid(self, rank): - if rank < self.ndims: - return False - return True + def get_stacked_values( + self, + absolutes: Tuple[List[float], List[float], List[float]], + ordinates: Tuple[List[float], List[float], List[float]], + weights: Optional[List[float]] = None, + ) -> Tuple[List[float], List[List[float]]]: + """Gathers stacked stacked absolutes/ordinates + NOTE: weights are only used in QRFactorization's child function + """ + # 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_weighted_values( self, values: Tuple[List[float], List[float], List[float]], - weights: Optional[List[float]], - ) -> Tuple[List[float], List[float], List[float]]: - """Application of weights for least squares methods, which calls for square rooting of weights + weights: Optional[List[float]] = None, + ) -> Union[List[float], List[List[float]]]: + """Application of weights for least squares methods, which calls for square roots Attributes ---------- @@ -81,31 +129,8 @@ class LeastSq(Transform): weights = np.vstack((weights, weights, weights)).T.ravel() 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: - """Calculates affine with no constraints using least squares.""" - abs_stacked, ord_stacked = self.get_stacked_values( - absolutes, ordinates, weights - ) - 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 - matrix, res, rank, sigma = spl.lstsq(ord_stacked.T, abs_stacked.T) - if self.valid(rank): - return self.get_matrix(matrix, absolutes, ordinates, weights) - print("Poorly conditioned or singular matrix, returning NaNs") - return np.nan * np.ones((4, 4)) - - def get_matrix(self, matrix, absolutes=None, ordinates=None, weights=None): - return np.array( - [ - [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], - ] - ) + def valid(self, rank: float) -> bool: + """ validates whether or not a matrix can reliably transform the method's number of dimensions """ + if rank < self.ndims: + return False + return True diff --git a/geomagio/adjusted/transform/QRFactorization.py b/geomagio/adjusted/transform/QRFactorization.py index b633ffb6602f02170f758878899799a0293751e7..598000e1a5c3548708de308ab0d2937468e39ba0 100644 --- a/geomagio/adjusted/transform/QRFactorization.py +++ b/geomagio/adjusted/transform/QRFactorization.py @@ -7,43 +7,19 @@ from .SVD import SVD class QRFactorization(LeastSq): - """Calculates affine using singular value decomposition with QR factorization""" + """Calculates affine using least squares with QR factorization""" - ndims = 2 - svd = SVD(ndims=ndims) + ndims: int = 2 + svd: SVD = SVD(ndims=ndims) - def get_weighted_values( + def get_matrix( self, - values: Tuple[List[float], List[float], List[float]], - weights: Optional[List[float]], - ): - """ Applies least squares weights in two dimensions(X and Y)""" - if weights is None: - return values - weights = np.sqrt(weights) - return np.array([values[i] * weights for i in range(self.ndims)]) - - def get_stacked_values(self, absolutes, ordinates, weights): - weighted_absolutes = self.svd.get_weighted_values( - values=absolutes, weights=weights - ) - weighted_ordinates = self.svd.get_weighted_values( - values=ordinates, weights=weights - ) - # LHS, or dependent variables - abs_stacked = self.svd.get_stacked_values( - values=absolutes, - weighted_values=weighted_absolutes, - ) - - # RHS, or independent variables - ord_stacked = self.svd.get_stacked_values( - values=ordinates, - weighted_values=weighted_ordinates, - ) - return abs_stacked, ord_stacked - - def get_matrix(self, matrix, absolutes, ordinates, weights): + matrix: List[List[float]], + absolutes: Tuple[List[float], List[float], List[float]], + ordinates: Tuple[List[float], List[float], List[float]], + weights: List[float], + ) -> np.array: + """ performs QR factorization steps and formats result within the returned matrix """ # QR fatorization # NOTE: forcing the diagonal elements of Q to be positive # ensures that the determinant is 1, not -1, and is @@ -78,3 +54,40 @@ class QRFactorization(LeastSq): ], [0.0, 0.0, 0.0, 1.0], ] + + def get_stacked_values( + self, + absolutes: Tuple[List[float], List[float], List[float]], + ordinates: Tuple[List[float], List[float], List[float]], + weights: Optional[List[float]] = None, + ) -> Tuple[List[List[float]], List[List[float]]]: + """ stacks and weights absolutes/ordinates """ + weighted_absolutes = self.svd.get_weighted_values( + values=absolutes, weights=weights + ) + weighted_ordinates = self.svd.get_weighted_values( + values=ordinates, weights=weights + ) + # LHS, or dependent variables + abs_stacked = self.svd.get_stacked_values( + values=absolutes, + weighted_values=weighted_absolutes, + ) + + # RHS, or independent variables + ord_stacked = self.svd.get_stacked_values( + values=ordinates, + weighted_values=weighted_ordinates, + ) + return abs_stacked, ord_stacked + + def get_weighted_values( + self, + values: Tuple[List[float], List[float], List[float]], + weights: Optional[List[float]] = None, + ) -> List[List[float]]: + """ Applies least squares weights in two dimensions(X and Y)""" + if weights is None: + return values + weights = np.sqrt(weights) + return np.array([values[i] * weights for i in range(self.ndims)]) diff --git a/geomagio/adjusted/transform/Rescale3D.py b/geomagio/adjusted/transform/Rescale3D.py index c74a9bb7258d67378b8e7bcc740da021f9066466..331b940a9018ee51968ed2fafddbc7b1a24dec47 100644 --- a/geomagio/adjusted/transform/Rescale3D.py +++ b/geomagio/adjusted/transform/Rescale3D.py @@ -1,13 +1,29 @@ import numpy as np import scipy.linalg as spl -from typing import List, Tuple +from typing import List, Optional, Tuple from .LeastSq import LeastSq class Rescale3D(LeastSq): """Calculates affine using using least squares, constrained to re-scale each axis""" - def get_stacked_ordinates(self, ordinates): + def get_matrix( + self, + matrix: List[List[float]], + absolutes: Optional[Tuple[List[float], List[float], List[float]]] = None, + ordinates: Optional[Tuple[List[float], List[float], List[float]]] = None, + weights: Optional[List[float]] = None, + ) -> np.array: + return [ + [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], + ] + + def get_stacked_ordinates( + self, ordinates: Tuple[List[float], List[float], List[float]] + ) -> List[List[float]]: # (reduces degrees of freedom by 13: # - 2 for making x independent of y,z; # - 2 for making y,z independent of x; @@ -20,11 +36,3 @@ class Rescale3D(LeastSq): ord_stacked[1, 1::3] = ordinates[1] ord_stacked[2, 2::3] = ordinates[2] return ord_stacked - - def get_matrix(self, matrix, absolutes=None, ordinates=None, weights=None): - return [ - [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/RotationTranslationXY.py b/geomagio/adjusted/transform/RotationTranslationXY.py index f140a38450e7905646b97643020956754bea4668..e6a352a15e03147083063b2eae623c18a9f7e420 100644 --- a/geomagio/adjusted/transform/RotationTranslationXY.py +++ b/geomagio/adjusted/transform/RotationTranslationXY.py @@ -9,12 +9,15 @@ class RotationTranslationXY(SVD): constrained to rotation and translation in XY(no scale or shear), and only translation in Z""" - ndims = 2 + ndims: int = 2 - def get_rotation_matrix(self, U, Vh): - return np.dot(Vh.T, np.dot(np.diag([1, np.linalg.det(np.dot(Vh.T, U.T))]), U.T)) - - def get_matrix(self, R, T, weighted_absolutes, weighted_ordinates): + def get_matrix( + self, + R: List[List[float]], + T: List[List[float]], + weighted_absolutes: Tuple[List[float], List[float], List[float]], + weighted_ordinates: Tuple[List[float], List[float], List[float]], + ) -> np.array: return [ [R[0, 0], R[0, 1], 0.0, T[0]], [R[1, 0], R[1, 1], 0.0, T[1]], @@ -27,3 +30,8 @@ class RotationTranslationXY(SVD): ], [0.0, 0.0, 0.0, 1.0], ] + + def get_rotation_matrix( + self, U: List[List[float]], Vh: List[List[float]] + ) -> List[List[float]]: + return np.dot(Vh.T, np.dot(np.diag([1, np.linalg.det(np.dot(Vh.T, U.T))]), U.T)) diff --git a/geomagio/adjusted/transform/SVD.py b/geomagio/adjusted/transform/SVD.py index f9ff490560f8a8766b73fb4c1b9ff85fef16ecf1..9972009cb357c475ec0906e809de4ed7f74d991d 100644 --- a/geomagio/adjusted/transform/SVD.py +++ b/geomagio/adjusted/transform/SVD.py @@ -7,7 +7,37 @@ from .Transform import Transform class SVD(Transform): """Instance of Transform. Applies singular value decomposition to generate matrices""" - def get_covariance_matrix(self, absolutes, ordinates, 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: + """Calculates matrix with singular value decomposition and accompanying methods + Defaults to singular value decomposition constrained for 3D rotation/translation + """ + weighted_ordinates = self.get_weighted_values(values=ordinates, weights=weights) + weighted_absolutes = self.get_weighted_values(values=absolutes, weights=weights) + # generate weighted "covariance" matrix + H = self.get_covariance_matrix(absolutes, ordinates, weights) + # 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) + if self.valid(S): + R = self.get_rotation_matrix(U, Vh) + # now get translation using weighted centroids and R + T = self.get_translation_matrix(R, weighted_absolutes, weighted_ordinates) + return self.get_matrix(R, T, weighted_absolutes, weighted_ordinates) + print("Poorly conditioned or singular matrix, returning NaNs") + return np.nan * np.ones((4, 4)) + + def get_covariance_matrix( + self, + absolutes: Tuple[List[float], List[float], List[float]], + ordinates: Tuple[List[float], List[float], List[float]], + weights: List[float], + ) -> List[List[float]]: + """ calculate covariance matrix with weighted absolutes/ordinates """ weighted_ordinates = self.get_weighted_values(values=ordinates, weights=weights) weighted_absolutes = self.get_weighted_values(values=absolutes, weights=weights) # generate weighted "covariance" matrix @@ -26,7 +56,40 @@ class SVD(Transform): ) return H - def get_stacked_values(self, values, weighted_values) -> np.array: + def get_matrix( + self, + R: List[List[float]], + T: List[List[float]], + weighted_absolutes: Optional[ + Tuple[List[float], List[float], List[float]] + ] = None, + weighted_ordinates: Optional[ + Tuple[List[float], List[float], List[float]] + ] = None, + ) -> np.array: + """Returns matrix formatted for 3D rotation/translation + NOTE: weighted absolutes/ordinates are only used by RotationTranslationXY's child function + """ + return [ + [R[0, 0], R[0, 1], R[0, 2], T[0]], + [R[1, 0], R[1, 1], R[1, 2], T[1]], + [R[2, 0], R[2, 1], R[2, 2], T[2]], + [0.0, 0.0, 0.0, 1.0], + ] + + def get_rotation_matrix( + self, U: List[List[float]], Vh: List[List[float]] + ) -> List[List[float]]: + """ computes rotation matrix from products of singular value decomposition """ + return np.dot( + Vh.T, np.dot(np.diag([1, 1, np.linalg.det(np.dot(Vh.T, U.T))]), U.T) + ) + + def get_stacked_values( + self, + values: Tuple[List[float], List[float], List[float]], + weighted_values: Tuple[List[float], List[float], List[float]], + ) -> np.array: """Supports intermediate mathematical steps by differencing and shaping values for SVD Attributes @@ -41,11 +104,22 @@ class SVD(Transform): """ return np.vstack([[values[i] - weighted_values[i]] for i in range(self.ndims)]) + def get_translation_matrix( + self, + R: List[List[float]], + weighted_absolutes: Tuple[List[float], List[float], List[float]], + weighted_ordinates: Tuple[List[float], List[float], List[float]], + ) -> List[List[float]]: + """ computes translation matrix from rotation matrix and weighted absolutes/ordinates """ + return np.array([weighted_absolutes[i] for i in range(self.ndims)]) - np.dot( + R, [weighted_ordinates[i] for i in range(self.ndims)] + ) + def get_weighted_values( self, values: Tuple[List[float], List[float], List[float]], weights: Optional[List[float]], - ) -> Tuple[List[float], List[float], List[float]]: + ) -> Tuple[float, float, float]: """Application of weights for SVD methods, which call for weighted averages""" if weights is None: weights = np.ones_like(values[0]) @@ -55,52 +129,8 @@ class SVD(Transform): np.average(values[2], weights=weights), ) - def get_rotation_matrix(self, U, Vh): - return np.dot( - Vh.T, np.dot(np.diag([1, 1, np.linalg.det(np.dot(Vh.T, U.T))]), U.T) - ) - - def get_translation_matrix(self, R, weighted_absolutes, weighted_ordinates): - return np.array([weighted_absolutes[i] for i in range(self.ndims)]) - np.dot( - R, [weighted_ordinates[i] for i in range(self.ndims)] - ) - - def valid(self, singular_values): + def valid(self, singular_values: List[float]) -> bool: + """ validates whether or not a matrix can reliably transform the method's number of dimensions """ if np.sum(singular_values) < self.ndims: return False return True - - def calculate( - self, - ordinates: Tuple[List[float], List[float], List[float]], - absolutes: Tuple[List[float], List[float], List[float]], - weights: List[float], - ) -> np.array: - weighted_ordinates = self.get_weighted_values(values=ordinates, weights=weights) - weighted_absolutes = self.get_weighted_values(values=absolutes, weights=weights) - # generate weighted "covariance" matrix - H = self.get_covariance_matrix(absolutes, ordinates, weights) - # 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) - if self.valid(S): - R = self.get_rotation_matrix(U, Vh) - # now get translation using weighted centroids and R - T = self.get_translation_matrix(R, weighted_absolutes, weighted_ordinates) - return self.get_matrix(R, T, weighted_absolutes, weighted_ordinates) - print("Poorly conditioned or singular matrix, returning NaNs") - return np.nan * np.ones((4, 4)) - - def get_matrix( - self, - R, - T, - weighted_absolutes=None, - weighted_ordinates=None, - ): - return [ - [R[0, 0], R[0, 1], R[0, 2], T[0]], - [R[1, 0], R[1, 1], R[1, 2], T[1]], - [R[2, 0], R[2, 1], R[2, 2], T[2]], - [0.0, 0.0, 0.0, 1.0], - ] diff --git a/geomagio/adjusted/transform/ShearYZ.py b/geomagio/adjusted/transform/ShearYZ.py index 7873c7b2934c1e167a5fa8c1c806f5b25c459967..afe3377a6317e85107d5a692178d2121f6e29d79 100644 --- a/geomagio/adjusted/transform/ShearYZ.py +++ b/geomagio/adjusted/transform/ShearYZ.py @@ -1,6 +1,6 @@ import numpy as np import scipy.linalg as spl -from typing import List, Tuple +from typing import List, Optional, Tuple from .LeastSq import LeastSq @@ -8,7 +8,23 @@ from .LeastSq import LeastSq class ShearYZ(LeastSq): """Calculates affine using least squares, constrained to shear y and z, but not x.""" - def get_stacked_ordinates(self, ordinates): + def get_matrix( + self, + matrix: List[List[float]], + absolutes: Optional[Tuple[List[float], List[float], List[float]]] = None, + ordinates: Optional[Tuple[List[float], List[float], List[float]]] = None, + weights: Optional[List[float]] = None, + ) -> np.array: + return [ + [1.0, 0.0, 0.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], + ] + + def get_stacked_ordinates( + self, ordinates: Tuple[List[float], List[float], List[float]] + ) -> List[List[float]]: # (reduces degrees of freedom by 13: # - 2 for making x independent of y,z; # - 1 for making y independent of z; @@ -22,11 +38,3 @@ class ShearYZ(LeastSq): ord_stacked[2, 1::3] = ordinates[1] ord_stacked[2, 2::3] = 1.0 return ord_stacked - - def get_matrix(self, matrix, absolutes=None, ordinates=None, weights=None): - return [ - [1.0, 0.0, 0.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/Transform.py b/geomagio/adjusted/transform/Transform.py index a365d9f9f040be5392ae346033be80b834590340..4d1f0c75bd92c657005f1ceda395930087d6c574 100644 --- a/geomagio/adjusted/transform/Transform.py +++ b/geomagio/adjusted/transform/Transform.py @@ -16,7 +16,23 @@ class Transform(BaseModel): acausal: bool = False memory: Optional[float] = None - ndims = 3 + ndims: int = 3 + + def calculate( + self, + ordinates: Tuple[List[float], List[float], List[float]], + absolutes: Tuple[List[float], List[float], List[float]], + weights: List[float], + ) -> np.array: + """Type skeleton inherited by any instance of Transform + + Attributes + ---------- + ordinates: H, E and Z ordinates + absolutes: X, Y and Z absolutes(NOTE: absolutes must be rotated from original H, E and Z values) + weights: time weights to apply during calculations of matrices + """ + return def get_weights(self, times: UTCDateTime, time: int = None) -> List[float]: """ @@ -49,19 +65,3 @@ class Transform(BaseModel): weights[times > time] = 0.0 return 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: - """Type skeleton inherited by any instance of Transform - - Attributes - ---------- - ordinates: H, E and Z ordinates - absolutes: X, Y and Z absolutes(NOTE: absolutes must be rotated from original H, E and Z values) - weights: time weights to apply during calculations of matrices - """ - return diff --git a/geomagio/adjusted/transform/TranslateOrigins.py b/geomagio/adjusted/transform/TranslateOrigins.py index fbb08b8f9f573f5eee5e182d7143c1078e71bf43..c10d16f4e16d74d55ef70dd21a77ebd4e5cadf1a 100644 --- a/geomagio/adjusted/transform/TranslateOrigins.py +++ b/geomagio/adjusted/transform/TranslateOrigins.py @@ -8,7 +8,35 @@ from .LeastSq import LeastSq class TranslateOrigins(LeastSq): """Calculates affine using using least squares, constrained to tanslate origins""" - def get_stacked_values(self, absolutes, ordinates, weights=None): + def get_matrix( + self, + matrix: List[List[float]], + absolutes: Optional[Tuple[List[float], List[float], List[float]]] = None, + ordinates: Optional[Tuple[List[float], List[float], List[float]]] = None, + weights: Optional[List[float]] = None, + ) -> np.array: + return [ + [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], + ] + + def get_stacked_ordinates( + self, ordinates: Tuple[List[float], List[float], List[float]] + ) -> List[List[float]]: + 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 + + def get_stacked_values( + self, + absolutes: Tuple[List[float], List[float], List[float]], + ordinates: Tuple[List[float], List[float], List[float]], + weights: Optional[List[float]] = None, + ) -> Tuple[List[float], List[List[float]]]: # LHS, or dependent variables abs_stacked = self.get_stacked_absolutes(absolutes) # subtract ords from abs to force simple translation @@ -23,7 +51,7 @@ class TranslateOrigins(LeastSq): self, values: Tuple[List[float], List[float], List[float]], weights: Optional[List[float]], - ): + ) -> List[float]: """Weights are applied after matrix creation steps, requiring weights to be stacked similar to ordinates and absolutes""" if weights is not None: @@ -32,18 +60,3 @@ class TranslateOrigins(LeastSq): else: weights = 1 return values * weights - - 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 - - def get_matrix(self, matrix, absolutes=None, ordinates=None, weights=None): - return [ - [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 4e56a6e835871b2457c08b5d94e7827c0e48ef91..1dc7724287bec0c2eec2c73e1c973ad3abd9a737 100644 --- a/geomagio/adjusted/transform/ZRotationHScale.py +++ b/geomagio/adjusted/transform/ZRotationHScale.py @@ -1,6 +1,6 @@ import numpy as np import scipy.linalg as spl -from typing import List, Tuple +from typing import List, Optional, Tuple from .LeastSq import LeastSq @@ -9,7 +9,23 @@ class ZRotationHscale(LeastSq): """Calculates affine using least squares, constrained to rotate about the Z axis and apply uniform horizontal scaling.""" - def get_stacked_ordinates(self, ordinates): + def get_matrix( + self, + matrix: List[List[float]], + absolutes: Optional[Tuple[List[float], List[float], List[float]]] = None, + ordinates: Optional[Tuple[List[float], List[float], List[float]]] = None, + weights: Optional[List[float]] = None, + ) -> np.array: + return [ + [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], + ] + + def get_stacked_ordinates( + self, ordinates: Tuple[List[float], List[float], List[float]] + ) -> List[List[float]]: # (reduces degrees of freedom by 10: # - 2 for making x,y independent of z; # - 2 for making z independent of x,y @@ -25,11 +41,3 @@ class ZRotationHscale(LeastSq): ord_stacked[4, 2::3] = ordinates[2] ord_stacked[5, 2::3] = 1.0 return ord_stacked - - def get_matrix(self, matrix, absolutes=None, ordinates=None, weights=None): - return [ - [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 13393ac38937db544986bdd8bc2483733adf6a04..cb72a6ed597a77b280970e553cb4cf735c2bd652 100644 --- a/geomagio/adjusted/transform/ZRotationHScaleZBaseline.py +++ b/geomagio/adjusted/transform/ZRotationHScaleZBaseline.py @@ -1,6 +1,6 @@ import numpy as np import scipy.linalg as spl -from typing import List, Tuple +from typing import List, Optional, Tuple from .LeastSq import LeastSq @@ -9,16 +9,23 @@ 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 get_stacked_values(self, absolutes, ordinates, weights=None): - # 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_matrix( + self, + matrix: List[List[float]], + absolutes: Optional[Tuple[List[float], List[float], List[float]]] = None, + ordinates: Optional[Tuple[List[float], List[float], List[float]]] = None, + weights: Optional[List[float]] = None, + ) -> np.array: + return [ + [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], + ] - def get_stacked_ordinates(self, ordinates): + def get_stacked_ordinates( + self, ordinates: Tuple[List[float], List[float], List[float]] + ) -> List[List[float]]: # (reduces degrees of freedom by 13: # - 2 for making x,y independent of z; # - 2 for making z independent of x,y; @@ -34,10 +41,16 @@ class ZRotationHscaleZbaseline(LeastSq): ord_stacked[2, 2::3] = 1.0 return ord_stacked - def get_matrix(self, matrix, absolutes=None, ordinates=None, weights=None): - return [ - [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], - ] + def get_stacked_values( + self, + absolutes: Tuple[List[float], List[float], List[float]], + ordinates: Tuple[List[float], List[float], List[float]], + weights: Optional[List[float]] = None, + ) -> Tuple[List[float], List[float]]: + # 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 diff --git a/geomagio/adjusted/transform/ZRotationShear.py b/geomagio/adjusted/transform/ZRotationShear.py index d1872c0e60e0f94c60afb549478668ad9cc84dbb..cb7adb3af21d4a6ab504444c06db95c74434c0c3 100644 --- a/geomagio/adjusted/transform/ZRotationShear.py +++ b/geomagio/adjusted/transform/ZRotationShear.py @@ -1,12 +1,30 @@ import numpy as np import scipy.linalg as spl -from typing import List, Tuple +from typing import List, Optional, Tuple from .LeastSq import LeastSq class ZRotationShear(LeastSq): - def get_stacked_ordinates(self, ordinates): + """Calculates affine using least squares, constrained to rotate about the Z axis """ + + def get_matrix( + self, + matrix: List[List[float]], + absolutes: Optional[Tuple[List[float], List[float], List[float]]] = None, + ordinates: Optional[Tuple[List[float], List[float], List[float]]] = None, + weights: Optional[List[float]] = None, + ) -> np.array: + return [ + [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], + ] + + def get_stacked_ordinates( + self, ordinates: Tuple[List[float], List[float], List[float]] + ) -> List[List[float]]: # (reduces degrees of freedom by 8: # - 2 for making x,y independent of z; # - 2 for making z independent of x,y @@ -21,11 +39,3 @@ class ZRotationShear(LeastSq): ord_stacked[6, 2::3] = ordinates[2] ord_stacked[7, 2::3] = 1.0 return ord_stacked - - def get_matrix(self, matrix, absolutes=None, ordinates=None, weights=None): - return [ - [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], - ]