import networkx as nx
import matplotlib.pyplot as plt
import numpy as np
import itertools
import math

MAX = 10
print("Launching with ", MAX, " iterations...")
import numpy as np

def write_to_file(P_sequences, filename="P_sequences.txt"):
    with open(filename, "w") as file:
        for idx, P_sequence in enumerate(P_sequences):
            file.write(f"P_sequence {idx + 1}:\n")
            np.savetxt(file, P_sequence, fmt="%.6f")
            file.write("\n")  # Add a blank line between matrices

def converge(v,p):
    i = 0
    while (np.max(np.abs(v - np.mean(v))) > 1):
        v = np.matmul(p, v)
        i += 1
    return i

def create_vector(max):
    np.random.seed(0)  # For reproducibility
    vector = np.random.randint(0, 1000, size=(max, 1))  # Random integers
    while np.any(np.abs(np.diff(vector.flatten())) < 100):  # Ensure the difference is at least 100
        vector = np.random.randint(0, 1000, size=(max, 1))
    return vector

# Creating graph
G = nx.Graph()
for i in range(MAX-1):
    G.add_edge(i+1, i+2)

# Initialize the matrix for each node
E_mat = {}
for i in G.nodes:
    # 0 either way for each value of the matrix
    E = np.zeros((MAX, MAX))


    for n in G.nodes:
        for p in G.nodes:
            # for each node -> n, then each node -> p
            # n is not a neightbor [v_n does not belong to N(v_i)]
            if n not in G.neighbors(i): 
                if p == n:
                    E[n - 1, p - 1] = 1
            # n is a neightbor [v_n belongs to N(v_i)]
            elif n in G.neighbors(i):
                if p == n or p == i:
                    E[n - 1, p - 1] = 0.5

    # Store the matrix for this node
    E_mat[i] = E

# Emission sequence
# For example sigma = sigma_1, sigma_2, ..., sigma_N we have :

# Creating multiple possible iterations
node_permutations = itertools.permutations(range(1, MAX + 1))

i = 0
P_sequences = []

# Eighenvalues high -> worst scenario, Eighenvalues low -> best scenario (lowest time)
lambda_2_values = []
convergence_times = []

for perm in node_permutations:
    # calculate sequence
    P_sequence = E_mat[perm[0]]
    for j in range(1,MAX):
        P_sequence = np.matmul(P_sequence, E_mat[perm[j]])
    
    # calculate eighenvalues and convergence time
    eigenvalues, _ = np.linalg.eig(P_sequence)
    # sort them then get best second value
    sorted_eigenvalues = np.sort(np.abs(eigenvalues))[::-1]
    lambda_2 = sorted_eigenvalues[1]
    # storing them by calculating the convergence time (1/1-lambda_2)
    lambda_2_values.append({'p_seq': P_sequence, 'val': lambda_2})
    #convergence_times.append(converge(vector_0,P_sequence))

    i = i+1

print("Theorical: ", math.factorial(MAX))
print("Final i: ",i)

# Create the vector v0
vector_0 = create_vector(MAX)

lambda_2_values.sort(key=lambda x: x['val'])
print("\n\nBest case: ", converge(vector_0,lambda_2_values[0]['p_seq']))
print("Worst case: ", converge(vector_0,lambda_2_values[len(lambda_2_values)-1]['p_seq']))
#write_to_file(P_sequences, filename="seq_"+str(MAX)+".txt")
    
# Plot
# plt.figure(figsize=(8, 6))
# plt.scatter(lambda_2_values, convergence_times, color="blue", alpha=0.7)
# plt.xlabel("Second Largest Eigenvalue (λ2)")
# plt.ylabel("Convergence Time (Cycles)")
# plt.title("Convergence Time vs. λ2")
# plt.grid(True)
# plt.show()