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

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions