calculate the sensitivity matrix for entry j of left or right eigenvector k given transition matrix T.
    Parameters
    ----------
    T : numpy.ndarray shape = (n, n)
        Transition matrix
    k : int
        eigenvector index ordered with descending eigenvalues
    j : int
        entry of eigenvector k for which the sensitivity is to be computed
    right : boolean (default: True)
        If set to True (default) the right eigenvectors are considered, otherwise the left ones

    Returns
    -------
    x : ndarray, shape=(n, n)
        Sensitivity matrix for entry index around transition matrix T. Reversibility is not assumed.

    Remarks
    -------
    Eigenvectors can naturally be scaled and so will their sensitivity depend on their size.
    For that reason we need to agree on a normalization for which the sensitivity is computed.
    Here we use the natural norm(vector) = 1 condition which is different from the results, e.g.
    the rdl_decomposition returns.
    This is especially important for the stationary distribution, which is the first left eigenvector.
    For this reason this function return a different sensitivity for the first left eigenvector
    than the function stationary_distribution_sensitivity and this function should not be used in this
    case!

def eigenvector_sensitivity(T, k, j, right=True):
    """
    calculate the sensitivity matrix for entry j of left or right eigenvector k given transition matrix T.
    Parameters
    ----------
    T : numpy.ndarray shape = (n, n)
        Transition matrix
    k : int
        eigenvector index ordered with descending eigenvalues
    j : int
        entry of eigenvector k for which the sensitivity is to be computed
    right : boolean (default: True)
        If set to True (default) the right eigenvectors are considered, otherwise the left ones

    Returns
    -------
    x : ndarray, shape=(n, n)
        Sensitivity matrix for entry index around transition matrix T. Reversibility is not assumed.

    Remarks
    -------
    Eigenvectors can naturally be scaled and so will their sensitivity depend on their size.
    For that reason we need to agree on a normalization for which the sensitivity is computed.
    Here we use the natural norm(vector) = 1 condition which is different from the results, e.g.
    the rdl_decomposition returns.
    This is especially important for the stationary distribution, which is the first left eigenvector.
    For this reason this function return a different sensitivity for the first left eigenvector
    than the function stationary_distribution_sensitivity and this function should not be used in this
    case!
    """

    n = len(T)

    if not right:
        T = numpy.transpose(T)

    eValues, rightEigenvectors = numpy.linalg.eig(T)
    leftEigenvectors = numpy.linalg.inv(rightEigenvectors)
    perm = numpy.argsort(eValues)[::-1]

    eValues = eValues[perm]
    rightEigenvectors = rightEigenvectors[:, perm]
    leftEigenvectors = leftEigenvectors[perm]

    rEV = rightEigenvectors[:, k]
    lEV = leftEigenvectors[k]
    eVal = eValues[k]

    vecA = numpy.zeros(n)
    vecA[j] = 1.0

    matA = T - eVal * numpy.identity(n)
    # Use here rEV as additional condition, means that we assume the vector to be
    # orthogonal to rEV
    matA = numpy.concatenate((matA, [rEV]))

    phi = numpy.linalg.lstsq(numpy.transpose(matA), vecA, rcond=-1)

    phi = numpy.delete(phi[0], -1)

    sensitivity = -numpy.outer(phi, rEV) + numpy.dot(phi, rEV) * numpy.outer(lEV, rEV)

    if not right:
        sensitivity = numpy.transpose(sensitivity)

    return sensitivity