import numpy as np

def shifting_redirection(M, eigenvalue, eigenvector):
    """
    Apply shifting redirection to the matrix to compute next eigenpair: M = M-lambda v
    """
    return(M-eigenvalue*np.matmul(eigenvector.T, eigenvector))

def power_method(M, epsilon = 0.0001, max_iter = 10000):
    """
    This function computes the principal component of M by using the power method with parameters:
    - epsilon: (float) Termination criterion to stop the power method when changes in the solution is marginale
    - max_iter: (int) Hard termination criterion
    Notes:
    - I added another condition based on the dot product of two consecutive solutions
    """
    # Initialization
    x = [None]*int(max_iter)
    x[0] = np.random.rand(M.shape[0])
    x[1] = np.matmul(M, x[0])
    count = 0

    # Compute eigenvector
    while((np.linalg.norm(x[count] - x[count-1]) > epsilon) & (count < max_iter)):
        # Actual computations
        x[count+1] = np.matmul(M, x[count])/np.linalg.norm(np.matmul(M, x[count]))
        count += 1

    # Compute eigenvalue
    eigenvalue = np.matmul(np.matmul(x[count].T, M), x[count])

    return(x[count], eigenvalue)

def eigenpairs(M, epsilon = 0.00001, max_iter = 10e2, plot = True):
    # Initialization
    eigenvectors = [None]*M.shape[0]
    eigenvalues = [None]*M.shape[0]

    for i in range(0, M.shape[0]):
        # Actual computing
        eigenvectors[i], eigenvalues[i] = power_method(M, epsilon, max_iter, iteration = i+1)
        M = shifting_redirection(M, eigenvalues[i], eigenvectors[i])

    return(eigenvectors, eigenvalues)