Skip to content

V2 enhancement: vectorized res & phase error calculation #177

@k-a-mendoza

Description

@k-a-mendoza

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

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type
    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions