#!/usr/bin/env python
"""
mtpy/utils/calculator.py
Helper functions for standard calculations, e.g. error propagation
@UofA, 2013
(LK)
"""
# =================================================================
import cmath
import math
import numpy as np
import mtpy.utils.exceptions as MTex
# =================================================================
# define uncertainty for differences between time steps
epsilon = 1e-9
# magnetic permeability in free space in H/m (=Vs/Am)
mu0 = 4e-7 * np.pi
# =================================================================
[docs]
def centre_point(xarray, yarray):
"""Get the centre point of arrays of x and y values."""
return (xarray.max() + xarray.min()) / 2.0, (yarray.max() + yarray.min()) / 2.0
[docs]
def roundsf(number, sf):
"""Round a number to a specified number of significant figures (sf)."""
# can't have < 1 s.f.
sf = max(sf, 1.0)
rounding = int(np.ceil(-np.log10(number) + sf - 1.0))
return np.round(number, rounding)
[docs]
def get_period_list(
period_min, period_max, periods_per_decade, include_outside_range=True
):
"""Get a list of values (e.g. periods), evenly spaced in log space and
including values on multiples of 10
:return s: Numpy array containing list of values.
"""
log_period_min = np.log10(period_min)
log_period_max = np.log10(period_max)
# check if log_period_min is a whole number
if log_period_min % 1 > 0:
# list of periods, around the minimum period, that will be present in specified
# periods per decade
aligned_logperiods_min = np.linspace(
np.floor(log_period_min),
np.ceil(log_period_min),
periods_per_decade + 1,
)
lpmin_diff = log_period_min - aligned_logperiods_min
# index of starting period, smallest value > 0
if include_outside_range:
spimin = np.where(lpmin_diff > 0)[0][-1]
else:
spimin = np.where(lpmin_diff < 0)[0][0]
start_period = aligned_logperiods_min[spimin]
else:
start_period = log_period_min
if log_period_max % 1 > 0:
# list of periods, around the maximum period, that will be present in specified
# periods per decade
aligned_logperiods_max = np.linspace(
np.floor(log_period_max),
np.ceil(log_period_max),
periods_per_decade + 1,
)
lpmax_diff = log_period_max - aligned_logperiods_max
# index of starting period, smallest value > 0
if include_outside_range:
spimax = np.where(lpmax_diff < 0)[0][0]
else:
spimax = np.where(lpmax_diff > 0)[0][-1]
stop_period = aligned_logperiods_max[spimax]
else:
stop_period = log_period_max
return np.logspace(
start_period,
stop_period,
int((stop_period - start_period) * periods_per_decade + 1),
)
[docs]
def nearest_index(val, array):
"""Find the index of the nearest value in the array.
:param val: The value to search for.
:param array: The array to search in.
:return: Index: integer describing position of nearest value in array.
"""
# absolute difference between value and array
diff = np.abs(array - val)
return np.where(diff == min(diff))[0][0]
[docs]
def make_log_increasing_array(z1_layer, target_depth, n_layers, increment_factor=0.999):
"""Create depth array with log increasing cells, down to target depth,
inputs are z1_layer thickness, target depth, number of layers (n_layers)
"""
# make initial guess for maximum cell thickness
max_cell_thickness = target_depth
# make initial guess for log_z
log_z = np.logspace(np.log10(z1_layer), np.log10(max_cell_thickness), num=n_layers)
counter = 0
while np.sum(log_z) > target_depth:
max_cell_thickness *= increment_factor
log_z = np.logspace(
np.log10(z1_layer), np.log10(max_cell_thickness), num=n_layers
)
counter += 1
if counter > 1e6:
break
return log_z
[docs]
def invertmatrix_incl_errors(inmatrix, inmatrix_error=None):
"""Invertmatrix incl errors."""
if inmatrix is None:
raise MTex.MTpyError_input_arguments("Matrix must be defined")
if (inmatrix_error is not None) and (inmatrix.shape != inmatrix_error.shape):
raise MTex.MTpyError_input_arguments(
"Matrix and err-matrix shapes do not match: %s - %s"
% (str(inmatrix.shape), str(inmatrix_error.shape))
)
if inmatrix.shape[-2] != inmatrix.shape[-1]:
raise MTex.MTpyError_input_arguments("Matrices must be square!")
if (inmatrix_error is not None) and (
inmatrix_error.shape[-2] != inmatrix_error.shape[-1]
):
raise MTex.MTpyError_input_arguments("Matrices must be square!")
dim = inmatrix.shape[-1]
det = np.linalg.det(inmatrix)
if det == 0:
raise MTex.MTpyError_input_arguments(
"Matrix is singular - I cannot invert that!"
)
inv_matrix = np.zeros_like(inmatrix)
if dim != 2:
raise MTex.MTpyError_input_arguments("Only 2D matrices supported yet")
inv_matrix = np.linalg.inv(inmatrix)
inv_matrix_err = None
if inmatrix_error is not None:
inmatrix_error = np.real(inmatrix_error)
inv_matrix_err = np.zeros_like(inmatrix_error)
for i in range(2):
for j in range(2):
# looping over the entries of the error matrix
err = 0.0
for k in range(2):
for l in range(2):
# each entry has 4 summands
err += np.abs(
-inv_matrix[i, k] * inv_matrix[l, j] * inmatrix_error[k, l]
)
inv_matrix_err[i, j] = err
return inv_matrix, inv_matrix_err
[docs]
def rhophi2z(rho, phi, freq):
"""Convert impedance-style information given in Rho/Phi format into
complex valued Z.
Input:
rho - 2x2 array (real) - in Ohm m
phi - 2x2 array (real) - in degrees
freq - scalar - frequency in Hz
Output:
Z - 2x2 array (complex)
"""
try:
if rho.shape != (2, 2) or phi.shape != (2, 2):
raise
if not (rho.dtype in ["float", "int"] and phi.dtype in ["float", "int"]):
raise
except:
raise MTex.MTpyError_input_arguments(
"ERROR - arguments must be two 2x2 arrays (real)"
)
z = np.zeros((2, 2), "complex")
for i in range(2):
for j in range(2):
abs_z = np.sqrt(5 * freq * rho[i, j])
z[i, j] = cmath.rect(abs_z, math.radians(phi[i, j]))
return z
[docs]
def compute_determinant_error(z_array, z_err_array, method="theoretical", repeats=1000):
"""Compute the error of the determinant of z using a stochastic method
seed random z arrays with a normal distribution around the input array
:param repeats:
Defaults to 1000.
:param z_array: Z (impedance) array containing real and imaginary values.
:param z_err_array: Impedance error array containing real values,
in MT we assume the real and imag errors are the same.
:param method: Method to use, theoretical calculation or stochastic, defaults to "theoretical".
:return: Error: array of real values with same shape as z_err_array
representing the error in the determinant of Z.
:return: Error_sqrt: array of real values with same shape as z_err_array
representing the error in the (determinant of Z)**0.5.
"""
if method == "stochatic":
arraylist = []
for r in range(repeats):
errmag = np.random.normal(loc=0, scale=z_err_array, size=z_array.shape)
arraylist = np.append(arraylist, z_array + errmag * (1.0 + 1j))
arraylist = arraylist.reshape(repeats, z_array.shape[0], 2, 2)
detlist = np.linalg.det(arraylist)
error = np.std(detlist, axis=0)
else:
error = np.abs(
z_err_array[:, 0, 0] * np.abs(z_array[:, 1, 1])
+ z_err_array[:, 1, 1] * np.abs(z_array[:, 0, 0])
- z_err_array[:, 0, 1] * np.abs(z_array[:, 1, 0])
- z_err_array[:, 1, 0] * np.abs(z_array[:, 0, 1])
)
return error
[docs]
def propagate_error_polar2rect(r, r_error, phi, phi_error):
"""Find error estimations for the transformation from polar to cartesian
coordinates.
Uncertainties in polar representation define a section of an annulus.
Find the 4 corners of this section and additionally the outer boundary
point, which is defined by phi = phi0, rho = rho0 + sigma rho.
The cartesian "box" defining the uncertainties in x,y is the outer
bound around the annulus section, defined by the four outermost
points. So check the four corners as well as the outer boundary edge
of the section to find the extrema in x znd y. These give you the
sigma_x/y.
"""
corners = [
(
np.real(cmath.rect(r - r_error, phi - phi_error)),
np.imag(cmath.rect(r - r_error, phi - phi_error)),
),
(
np.real(cmath.rect(r + r_error, phi - phi_error)),
np.imag(cmath.rect(r + r_error, phi - phi_error)),
),
(
np.real(cmath.rect(r + r_error, phi + phi_error)),
np.imag(cmath.rect(r + r_error, phi + phi_error)),
),
(
np.real(cmath.rect(r - r_error, phi + phi_error)),
np.imag(cmath.rect(r - r_error, phi + phi_error)),
),
(
np.real(cmath.rect(r + r_error, phi)),
np.imag(cmath.rect(r + r_error, phi)),
),
]
lo_x = [i[0] for i in corners]
lo_y = [i[1] for i in corners]
point = (np.real(cmath.rect(r, phi)), np.imag(cmath.rect(r, phi)))
lo_xdiffs = [abs(point[0] - i) for i in lo_x]
lo_ydiffs = [abs(point[1] - i) for i in lo_y]
xerr = max(lo_xdiffs)
yerr = max(lo_ydiffs)
return xerr, yerr
[docs]
def propagate_error_rect2polar(x, x_error, y, y_error):
"""Propagate error rect2polar."""
# x_error, y_error define a rectangular uncertainty box
# rho error is the difference between the closest and furthest point of the box (w.r.t. the origin)
# approximation: just take corners and midpoint of edges
lo_points = [
(x + x_error, y),
(x - x_error, y),
(x, y - y_error),
(x, y + y_error),
(x - x_error, y - y_error),
(x + x_error, y - y_error),
(x + x_error, y + y_error),
(x - x_error, y + y_error),
]
# check, if origin is within the box:
origin_in_box = False
if x_error >= np.abs(x) and y_error >= np.abs(y):
origin_in_box = True
try:
lo_polar_points = [cmath.polar(complex(*i)) for i in lo_points]
except OverflowError:
lo_polar_points = []
for i in lo_points:
real, imag = i
if abs(real) > 1e32:
real = 1e32
if abs(imag) > 1e32:
imag = 1e32
lo_polar_points.append(cmath.polar(complex(real, imag)))
lo_rho = [i[0] for i in lo_polar_points]
lo_phi = [math.degrees(i[1]) % 360 for i in lo_polar_points]
rho_err = 0.5 * (max(lo_rho) - min(lo_rho))
phi_err = 0.5 * (max(lo_phi) - min(lo_phi))
if (270 < max(lo_phi) < 360) and (0 < min(lo_phi) < 90):
tmp1 = [i for i in lo_phi if (0 < i < 90)]
tmp4 = [i for i in lo_phi if (270 < i < 360)]
phi_err = 0.5 * ((max(tmp1) - min(tmp4)) % 360)
if phi_err > 180:
phi_err = (-phi_err) % 360
if origin_in_box is True:
# largest rho:
rho_err = 2 * rho_err + min(lo_rho)
# maximum angle uncertainty:
phi_err = 180.0
return rho_err, phi_err
[docs]
def z_error2r_phi_error(z_real, z_imag, error):
"""Error estimation from rectangular to polar coordinates.
By standard error propagation, relative error in resistivity is
2*relative error in z amplitude.
Uncertainty in phase (in degrees) is computed by defining a circle around
the z vector in the complex plane. The uncertainty is the absolute angle
between the vector to (x,y) and the vector between the origin and the
tangent to the circle.
:return s: Tuple containing relative error in resistivity, absolute error in phase.
"""
z_real = np.array(z_real)
z_imag = np.array(z_imag)
error = np.array(error)
z_amp = np.abs(z_real + 1j * z_imag)
if (z_amp == 0).all():
z_rel_err = np.array([0.0])
else:
z_rel_err = error / z_amp
res_rel_err = 2.0 * z_rel_err
# if the relative error of the amplitude is >=100% that means that the relative
# error of the resistivity is 200% - that is then equivalent to an uncertainty
# in the phase angle of 90 degrees:
phi_err = np.degrees(np.arctan(z_rel_err))
phi_err[res_rel_err > 1.0] = 90.0
return res_rel_err, phi_err
# rotation:
# 1. rotation positive in clockwise direction
# 2. orientation of new X-axis X' given by rotation angle
# 3. express contents of Z/tipper (points P) in this new system (points P')
# 4. rotation for points calculated as P' = ([cos , sin ],[-sin, cos]) * P <=> P' = R * P
# 5. => B' = R * B and E' = R * E
# (Rt is the inverse rotation matrix)
# 6. E = Z * B => Rt * E' = Z * Rt * B' => E' = (R*Z*Rt) * B' => Z' = (R*Z*Rt)
# 7. Bz = T * B => Bz = T * Rt * B' => T' = (T * Rt)
# (Since Tipper is a row vector)
# 7.a for general column vectors v: v' = R * v
# Rotation of the uncertainties:
# a) rotate Z into Z'
# b) use propagation of errors on Z' to obtain the rotated Z'err
# That is NOT the same as the rotated error matrix Zerr (although the result is similar)
[docs]
def get_rotation_matrix(angle, clockwise=False):
"""
get the rotation matrix for the proper rotation
:param angle: angle in degrees to rotate by
:type angle: float
:param clockwise: [ True | False ] if False counterclockwise rotation
matrix is returned
:type clockwise: bool
"""
angle = np.deg2rad(angle)
cphi = np.cos(angle)
sphi = np.sin(angle)
if clockwise:
return np.array([[cphi, sphi], [-sphi, cphi]])
else:
return np.array([[cphi, -sphi], [sphi, cphi]])
[docs]
def rotate_matrix_errors(error, angle):
"""
Rotate errors of a matrix
:param error: error matrix (2 x 2)
:type error: np.array
:param angle: rotation angle in degrees
:type angle: float
:return: propogated errors
:rtype: None or np.array
"""
try:
phi = np.deg2rad(float(angle) % 360)
except TypeError:
raise MTex.MTpyError_inputarguments('"Angle" must be a float')
cphi = np.cos(phi)
sphi = np.sin(phi)
error_matrix = None
if error is not None:
err_orig = np.real(error)
error_matrix = np.zeros_like(error)
# standard propagation of errors:
error_matrix[0, 0] = np.sqrt(
(cphi**2 * err_orig[0, 0]) ** 2
+ (cphi * sphi * err_orig[0, 1]) ** 2
+ (cphi * sphi * err_orig[1, 0]) ** 2
+ (sphi**2 * err_orig[1, 1]) ** 2
)
error_matrix[0, 1] = np.sqrt(
(cphi**2 * err_orig[0, 1]) ** 2
+ (cphi * sphi * err_orig[1, 1]) ** 2
+ (cphi * sphi * err_orig[0, 0]) ** 2
+ (sphi**2 * err_orig[1, 0]) ** 2
)
error_matrix[1, 0] = np.sqrt(
(cphi**2 * err_orig[1, 0]) ** 2
+ (cphi * sphi * err_orig[1, 1]) ** 2
+ (cphi * sphi * err_orig[0, 0]) ** 2
+ (sphi**2 * err_orig[0, 1]) ** 2
)
error_matrix[1, 1] = np.sqrt(
(cphi**2 * err_orig[1, 1]) ** 2
+ (cphi * sphi * err_orig[0, 1]) ** 2
+ (cphi * sphi * err_orig[1, 0]) ** 2
+ (sphi**2 * err_orig[0, 0]) ** 2
)
return error_matrix
[docs]
def rotate_matrix_with_errors(in_matrix, angle, error=None, clockwise=True):
"""Rotate a matrix including errors given an angle in degrees.
:param in_matrix: A n x 2 x 2 matrix to rotate.
:type in_matrix: np.ndarray
:param angle: Angle to rotate by assuming clockwise positive from
0 = north.
:type angle: float
:param error: A n x 2 x 2 matrix of associated errors,, defaults to None.
:type error: np.ndarray, optional
:param clockwise: rotate clockwise [True] or counter-clockwise [False]
:type clockwise: bool
:raises MTex: If input array is incorrect.
:return: Rotated matrix.
:rtype: np.ndarray
:return: Rotated matrix errors.
:rtype: np.ndarray
"""
if in_matrix is None:
raise MTex.MTpyError_inputarguments("Matrix must be defined")
if (error is not None) and (in_matrix.shape != error.shape):
msg = "matricies are not the same shape in_matrix={0}, err={1}".format(
in_matrix.shape, error.shape
)
raise MTex.MTpyError_inputarguments(msg)
rot_mat = get_rotation_matrix(angle, clockwise)
# jacobian rotation or similarity transformation is R A RT
# rotated_matrix = np.dot(np.dot(rot_mat, in_matrix), rot_mat.T)
rotated_matrix = rot_mat @ in_matrix @ rot_mat.T
error_matrix = rotate_matrix_errors(error, angle)
return rotated_matrix, error_matrix
[docs]
def rotate_vector_with_errors(in_vector, angle, error=None, clockwise=True):
"""Rotate a vector including errors given an angle in degrees.
:param in_vector:
:param in_matrix: A n x 1 x 2 vector to rotate.
:type in_matrix: np.ndarray
:param angle: Angle to rotate by assuming clockwise positive from
0 = north.
:type angle: float
:param error: A n x 1 x 2 vector of associated errors,, defaults to None.
:type error: np.ndarray, optional
:raises MTex: If input array is incorrect.
:return: Rotated vector.
:rtype: np.ndarray
:return: Rotated vector errors.
:rtype: np.ndarray
"""
if in_vector is None:
raise MTex.MTpyError_inputarguments(
"Vector AND error-vector must" + " be defined"
)
if (error is not None) and (in_vector.shape != error.shape):
msg = "matricies are not the same shape in_vector={0}, err={1}".format(
in_vector.shape, error.shape
)
raise MTex.MTpyError_inputarguments(msg)
rot_mat = get_rotation_matrix(angle, clockwise)
if in_vector.shape == (1, 2):
rotated_vector = np.matmul(in_vector, np.linalg.inv(rot_mat))
else:
rotated_vector = np.matmul(rot_mat, in_vector)
err_vector = None
if error is not None:
err_vector = np.zeros_like(error)
if error.shape == (1, 2):
err_vector = np.matmul(error, np.abs(np.linalg.inv(rot_mat)))
else:
err_vector = np.matmul(np.abs(rot_mat), error)
return rotated_vector, err_vector
[docs]
def multiplymatrices_incl_errors(
inmatrix1, inmatrix2, inmatrix1_error=None, inmatrix2_error=None
):
"""Multiplymatrices incl errors."""
if inmatrix1 is None or inmatrix2 is None:
raise MTex.MTpyError_input_arguments("ERROR - two 2x2 arrays needed as input")
if inmatrix1.shape != inmatrix2.shape:
raise MTex.MTpyError_input_arguments(
"ERROR - two 2x2 arrays with same dimensions needed as input"
)
prod = np.matmul(np.asarray(inmatrix1), np.asarray(inmatrix2))
if (inmatrix1_error is None) or (inmatrix1_error is None):
return prod, None
var = np.zeros((2, 2))
var[0, 0] = (
(inmatrix1_error[0, 0] * inmatrix2[0, 0]) ** 2
+ (inmatrix1_error[0, 1] * inmatrix2[1, 0]) ** 2
+ (inmatrix2_error[0, 0] * inmatrix1[0, 0]) ** 2
+ (inmatrix2_error[1, 0] * inmatrix1[0, 1]) ** 2
)
var[0, 1] = (
(inmatrix1_error[0, 0] * inmatrix2[0, 1]) ** 2
+ (inmatrix1_error[0, 1] * inmatrix2[1, 1]) ** 2
+ (inmatrix2_error[0, 1] * inmatrix1[0, 0]) ** 2
+ (inmatrix2_error[1, 1] * inmatrix1[0, 1]) ** 2
)
var[1, 0] = (
(inmatrix1_error[1, 0] * inmatrix2[0, 0]) ** 2
+ (inmatrix1_error[1, 1] * inmatrix2[1, 0]) ** 2
+ (inmatrix2_error[0, 0] * inmatrix1[1, 0]) ** 2
+ (inmatrix2_error[1, 0] * inmatrix1[1, 1]) ** 2
)
var[1, 1] = (
(inmatrix1_error[1, 0] * inmatrix2[0, 1]) ** 2
+ (inmatrix1_error[1, 1] * inmatrix2[1, 1]) ** 2
+ (inmatrix2_error[0, 1] * inmatrix1[1, 0]) ** 2
+ (inmatrix2_error[1, 1] * inmatrix1[1, 1]) ** 2
)
return prod, np.sqrt(var)
[docs]
def reorient_data2D(x_values, y_values, x_sensor_angle=0, y_sensor_angle=90):
"""Re-orient time series data of a sensor pair, which has not been in default (x=0, y=90) orientation.
Input:
- x-values - Numpy array
- y-values - Numpy array
Note: same length for both! - If not, the shorter length is taken
Optional:
- Angle of the x-sensor - measured in degrees, clockwise from North (0)
- Angle of the y-sensor - measured in degrees, clockwise from North (0)
Output:
- corrected x-values (North)
- corrected y-values (East)
"""
x_values = np.array(x_values)
y_values = np.array(y_values)
try:
if x_values.dtype not in ["complex", "float", "int"]:
raise
if len(x_values) != len(y_values):
raise
except:
raise MTex.MTpyError_input_arguments(
"ERROR - both input arrays must be of same length"
)
if len(x_values) != len(y_values):
l = min(len(x_values), len(y_values))
x_values = x_values[:l]
y_values = y_values[:l]
in_array = np.zeros((len(x_values), 2), x_values.dtype)
in_array[:, 0] = x_values
in_array[:, 1] = y_values
try:
x_angle = math.radians(x_sensor_angle)
y_angle = math.radians(y_sensor_angle)
except:
raise MTex.MTpyError_input_arguments(
"ERROR - both angles must be of type int or float"
)
T = np.matrix(
[
[np.real(cmath.rect(1, x_angle)), np.imag(cmath.rect(1, x_angle))],
[np.real(cmath.rect(1, y_angle)), np.imag(cmath.rect(1, y_angle))],
]
)
try:
new_array = np.dot(in_array, T.I)
except:
raise MTex.MTpyError_input_arguments(
"ERROR - angles must define independent axes to span 2D"
)
# print new_array.shape
return new_array[:, 0], new_array[:, 1]