#%%
import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(seed=1)
n = 5

maxIter = 1001 # numarul maxim de iteratii
toleranta = 0.0001 # nivelul de toleranta

y = np.random.rand(n) # se genereaza random un vector pentru vectorul propriu
ySave = y # salvez y initial pentru a putea fi folosit la ambele metode

def powerMethod(M):
    '''
        Metoda primeste parametrii M si y\n
        M - matricea pe care se lucreaza\n
        y - vector propriu asociat valorii proprii dominante
    '''

    y = np.random.rand(n) # se genereaza random un vector pentru vectorul propriu

    # Se initializeaza un vector gol pentru erorile pe care le vom calcula mai tarziu
    erori = np.zeros(maxIter)

    y = y / np.linalg.norm(y)

    i = 0 # numarul curent de iteratii
    e = 1 # eroarea initiala

    # Se calculeaza vectorul propriu unitar asociat valorii proprii a matricei
    while e > toleranta:
        if i > maxIter:
            print("S-a atins numarul maxim de iteratii fara a se fi obtinut nivelul prescris al tolerantei")
            break

        z = np.dot(M, y)
        z = z / np.linalg.norm(z)

        zTrans = np.transpose(z)
        e = abs(1 - abs(zTrans @ y)) # se calculeaza eroarea
        erori[i] = e # se face push in vectorul de erori la eraoarea curenta

        y = z # se actualizeaza vectorul propriu

        i += 1

    print("Vectorul propriu:\n", y)

    return erori


def inversePowerMethod(M):
    '''
        Metoda primeste parametrii M si y\n
        M - matricea pe care se lucreaza\n
        y - vector propriu asociat valorii proprii dominante
    '''

    y = ySave / np.linalg.norm(ySave)
    i = 0
    e = 1

    erori = np.zeros(maxIter)

    while e > toleranta:
        if i > maxIter:
            print("S-a atins numarul maxim de iteratii fara a se fi obtinut nivelul prescris al tolerantei")
            break

        miu = np.transpose(y) @ M @ y

        In = np.eye(n)
        z = np.linalg.solve(miu * In - M, y)

        z = z / np.linalg.norm(z)

        zTrans = np.transpose(z)
        e = abs(1 - abs(zTrans @ y)) # se calculeaza eroarea
        erori[i] = e # se face push in vectorul de erori la eraoarea curenta

        y = z # se actualizeaza vectorul propriu
        i += 1

    print("Vectorul propriu:\n", y)

    return erori



M1 = rng.integers(0, 20, size=(n, n))
print("Matricea generata:\n", M1)

w, v = np.linalg.eig(M1)
print("V:\n", v)

erori1 = powerMethod(M1)
plt.plot(erori1)
plt.title("Eroarea pentru metoda puterii")
plt.show()

erori1 = inversePowerMethod(M1)
print("V:\n", v)
plt.plot(erori1)
plt.title("Eroarea pentru metoda puterii inverse")
plt.show()

M2 = [
    [1.556, 0.6, 0.8, 6.3, 4.4],
    [0, 1.344, 0.8, 6.1, 8.2],
    [0, 0, 1.345, 0.7, 5.3],
    [0, 0, 0, 1.356, 1.25],
    [0, 0, 0, 0, 1.681]
]

w, v = np.linalg.eig(M2)
print("V:\n", v)

erori2 = powerMethod(M2)
plt.plot(erori2)
plt.title("Eroarea pentru metoda puterii")
plt.show()

erori2 = inversePowerMethod(M2)
print("V:\n", v)
plt.plot(erori2)
plt.title("Eroarea pentru metoda puterii inverse")
plt.show()

# %%