docs: 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