<commit_before>"""
Pixel-by-pixel linear rate (velocity) estimation using iterative weighted
least-squares method.

Based on original Matlab code by <NAME> and <NAME>, and Matlab 'lscov'
function.

.. codeauthor: <NAME> and <NAME>, GA
"""

from scipy.linalg import solve, cholesky, qr, inv
from numpy import nan, isnan, sqrt, diag, delete, ones, array, nonzero, float32
import numpy as np


def is_pos_def(x):
    return np.all(np.linalg.eigvals(x) > 1e-6)


def linear_rate(ifgs, vcm, pthr, nsig, maxsig, mst=None):
    """
    Pixel-by-pixel linear rate (velocity) estimation using iterative weighted least-squares method.

    :param ifgs: Sequence of ifg objs from which to extract observations
    :param vcm: Derived positive definite temporal variance covariance matrix
    :param pthr: Pixel threshold; minimum number of coherent observations for a pixel
    :param nsig: n-sigma ratio used to threshold 'model minus observation' residuals
    :param maxsig: Threshold for maximum allowable standard error
    :param mst: Pixel-wise matrix describing the minimum spanning tree network

    :return:
        python/matlab variable names
        rate/ifg_stack: stacked interferogram (i.e., rate map)
        error/std_stack: standard deviation of the stacked interferogram
                  (i.e., error map)
        samples/coh_sta: statistics of coherent pixels used for stacking
        demerror:  dem errors in metres, no implemented in python
    """

    rows, cols = ifgs[0].phase_data.shape

    # make 3D block of observations
    obs = array([np.where(isnan(x.phase_data), 0, x.phase_data) for x in ifgs])
    span = array([[x.time_span for x in ifgs]])

    # Update MST in case additional NaNs generated by APS filtering
    if mst is None:  # dummy mst if none is passed in
        mst = ~isnan(obs)
    else:
        mst[isnan(obs)] = 0

    # preallocate NaN arrays
    error = np.empty([rows, cols], dtype=float32) * np.nan
    rate = np.empty([rows, cols], dtype=float32) * np.nan
    samples = np.empty([rows, cols], dtype=float32) * np.nan

    # pixel-by-pixel calculation.
    # nested loops to loop over the 2 image dimensions
    for i in xrange(rows):
        for j in xrange(cols):

            # find the indices of independent ifgs for given pixel from MST
            ind = np.nonzero(mst[:, i, j])[0]  # only True's in mst are chosen

            # iterative loop to calculate 'robust' velocity for pixel

            while len(ind) >= pthr:
                # make vector of selected ifg observations
                ifgv = obs[ind, i, j]

                # form design matrix from appropriate ifg time spans
                B = span[:, ind]

                # Subset of full VCM matrix for selected observations
                vcm_temp = vcm[ind, np.vstack(ind)]

                """ start matlab lscov routine """
                
                # Get the lower triangle cholesky decomposition.
                # V must be positive definite (symmetrical and square)
                T = cholesky(vcm_temp, 1)

                # Incorporate inverse of VCM into the design matrix and observations vector
                A = solve(T, B.transpose())
                b = solve(T, ifgv.transpose())

                # Factor the design matrix, incorporate covariances or weights into the
                # system of equations, and transform the response vector.
                Q, R, P = qr(A, mode='economic', pivoting=True)
                z = Q.conj().transpose().dot(b)

                # Compute the Lstsq coefficient for the velocity
                v = solve(R, z)

                """end matlab lscov routine"""

                # Compute the model errors; added by <NAME>, 12/12/2011
                err1 = inv(vcm_temp).dot(B.conj().transpose())
                err2 = B.dot(err1)
                err = sqrt(diag(inv(err2)))

                # Compute the residuals (model minus observations)
                r = (B * v) - ifgv

                # determine the ratio of residuals and apriori variances (Danish method)
                w = cholesky(inv(vcm_temp))
                wr = abs(np.dot(w, r.transpose()))

                # test if maximum ratio is greater than user threshold.
                max_val = wr.max()

                if max_val > nsig:
                    # if yes, discard and re-do the calculation.
                    ind = delete(ind, wr.argmax())
                else:
                    # if no, save estimate, exit the while loop and go to next pixel
                    rate[i, j] = v
                    error[i, j] = err
                    samples[i, j] = ifgv.shape[0]
                    break


    # overwrite the data whose error is larger than the maximum sigma user threshold
    rate[error > maxsig] = nan
    error[error > maxsig] = nan
    samples[error > maxsig] = nan  # TODO: This step is missing in matlab?

    return rate, error, samples

if __name__ == "__main__":
    import os
    import shutil
    from subprocess import call<commit_msg>[linrate] WIP: parallel linrate main function<commit_after>"""
Pixel-by-pixel linear rate (velocity) estimation using iterative weighted
least-squares method.

Based on original Matlab code by <NAME> and <NAME>, and Matlab 'lscov'
function.

.. codeauthor: <NAME> and <NAME>, GA
"""

from scipy.linalg import solve, cholesky, qr, inv
from numpy import nan, isnan, sqrt, diag, delete, ones, array, nonzero, float32
import numpy as np
import parmap


def is_pos_def(x):
    return np.all(np.linalg.eigvals(x) > 1e-6)


def linear_rate(ifgs, vcm, pthr, nsig, maxsig, mst=None):
    """
    Pixel-by-pixel linear rate (velocity) estimation using iterative weighted least-squares method.

    :param ifgs: Sequence of ifg objs from which to extract observations
    :param vcm: Derived positive definite temporal variance covariance matrix
    :param pthr: Pixel threshold; minimum number of coherent observations for a pixel
    :param nsig: n-sigma ratio used to threshold 'model minus observation' residuals
    :param maxsig: Threshold for maximum allowable standard error
    :param mst: Pixel-wise matrix describing the minimum spanning tree network

    :return:
        python/matlab variable names
        rate/ifg_stack: stacked interferogram (i.e., rate map)
        error/std_stack: standard deviation of the stacked interferogram
                  (i.e., error map)
        samples/coh_sta: statistics of coherent pixels used for stacking
        demerror:  dem errors in metres, no implemented in python
    """

    rows, cols = ifgs[0].phase_data.shape

    # make 3D block of observations
    obs = array([np.where(isnan(x.phase_data), 0, x.phase_data) for x in ifgs])
    span = array([[x.time_span for x in ifgs]])

    # Update MST in case additional NaNs generated by APS filtering
    if mst is None:  # dummy mst if none is passed in
        mst = ~isnan(obs)
    else:
        mst[isnan(obs)] = 0

    # preallocate NaN arrays
    error = np.empty([rows, cols], dtype=float32) * np.nan
    rate = np.empty([rows, cols], dtype=float32) * np.nan
    samples = np.empty([rows, cols], dtype=np.int16) * np.nan

    # pixel-by-pixel calculation.
    # nested loops to loop over the 2 image dimensions
    for i in xrange(rows):
        for j in xrange(cols):
            res = linear_rate_by_pixel(error, i, j, mst, nsig, obs, pthr, rate, samples, span, vcm)
            rate[i, j], error[i, j], samples[i, j] = res

    # overwrite the data whose error is larger than the maximum sigma user threshold
    rate[error > maxsig] = nan
    error[error > maxsig] = nan
    samples[error > maxsig] = nan  # TODO: This step is missing in matlab?

    return rate, error, samples


def linear_rate_by_pixel(error, i, j, mst, nsig, obs, pthr, rate, samples, span, vcm):
    # find the indices of independent ifgs for given pixel from MST
    ind = np.nonzero(mst[:, i, j])[0]  # only True's in mst are chosen
    # iterative loop to calculate 'robust' velocity for pixel

    while len(ind) >= pthr:
        # make vector of selected ifg observations
        ifgv = obs[ind, i, j]

        # form design matrix from appropriate ifg time spans
        B = span[:, ind]

        # Subset of full VCM matrix for selected observations
        vcm_temp = vcm[ind, np.vstack(ind)]

        """ start matlab lscov routine """

        # Get the lower triangle cholesky decomposition.
        # V must be positive definite (symmetrical and square)
        T = cholesky(vcm_temp, 1)

        # Incorporate inverse of VCM into the design matrix and observations vector
        A = solve(T, B.transpose())
        b = solve(T, ifgv.transpose())

        # Factor the design matrix, incorporate covariances or weights into the
        # system of equations, and transform the response vector.
        Q, R, _ = qr(A, mode='economic', pivoting=True)
        z = Q.conj().transpose().dot(b)

        # Compute the Lstsq coefficient for the velocity
        v = solve(R, z)

        """end matlab lscov routine"""

        # Compute the model errors; added by <NAME>, 12/12/2011
        err1 = inv(vcm_temp).dot(B.conj().transpose())
        err2 = B.dot(err1)
        err = sqrt(diag(inv(err2)))

        # Compute the residuals (model minus observations)
        r = (B * v) - ifgv

        # determine the ratio of residuals and apriori variances (Danish method)
        w = cholesky(inv(vcm_temp))
        wr = abs(np.dot(w, r.transpose()))

        # test if maximum ratio is greater than user threshold.
        max_val = wr.max()
        if max_val > nsig:
            # if yes, discard and re-do the calculation.
            ind = delete(ind, wr.argmax())
        else:
            # if no, save estimate, exit the while loop and go to next pixel
            return v[0], err[0], ifgv.shape[0]
    # dummy return for no change
    return np.nan, np.nan, np.nan


if __name__ == "__main__":
    import os
    import shutil
    from subprocess import call
{'fraction_non_alphanumeric': 0.08973460190285427, 'fraction_numerical': 0.006209313970956434, 'mean_word_length': 2.3498825897349884, 'pattern_counts': {'":': 2, '<': 13, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '462345', 'n_tokens_mistral': 3077, 'n_tokens_neox': 2850, 'n_words': 1099}