Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

V2 enhancement: vectorized res & phase error calculation #177

Open
k-a-mendoza opened this issue Oct 26, 2022 · 0 comments
Open

V2 enhancement: vectorized res & phase error calculation #177

k-a-mendoza opened this issue Oct 26, 2022 · 0 comments

Comments

@k-a-mendoza
Copy link

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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant