Lately I've been digging into error propagation due to its relevance to inversion data weighting. Tracing back MtPy's routines, it seems many of the algorithms are a mix of cmath with lots of unnecessary index lookups. Here's a version of Mtpy.core.z.py (lines 102 - 129) and mtpy.utils.calculator (lines 347-389) which vectorizes the computation
in mtpy.utils.calculator
def z_error2r_phi_error(z, z_err):
"""
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.
:returns:
tuple containing relative error np.ndarray in resistivity, absolute error np.ndarray in phase (degrees)
:inputs:
z, np.ndarray complex valued impedance tensor
z_err, np.ndarray float representing the stdev error
"""
relative_z_err = self._z_err/np.abs(self._z)
relative_res_err = 2*relative_z_err
phi_err = np.degrees(np.arctan(relative_z_err))
phi_err[relative_res_err > 1.] = 90
return relative_res_err, phi_err
in mtpy.core.z
def compute_resistivity_phase(self, z_array=None, z_err_array=None,
freq=None):
... extra lines here which don't need to be changed....
self._resistivity = np.abs(self._z)** 2 / (5*self.freq[:,np.newaxis,np.newaxis])
self._phase = np.rad2deg(np.angle(self._z))
self._resistivity_err = np.zeros_like(self._resistivity, dtype=np.float)
self._phase_err = np.zeros_like(self._phase, dtype=np.float)
if self._z_err is not None:
res_err, phase_err = Mtcc.z_err2r_phi_err(self._z,self._z_err)
self._resistivity_err = res_err*self._resistivity
self._phase_err = phase_err
Lately I've been digging into error propagation due to its relevance to inversion data weighting. Tracing back MtPy's routines, it seems many of the algorithms are a mix of cmath with lots of unnecessary index lookups. Here's a version of Mtpy.core.z.py (lines 102 - 129) and mtpy.utils.calculator (lines 347-389) which vectorizes the computation
in mtpy.utils.calculator
in mtpy.core.z