import numpy as np
import random
random.seed;
import scipy
from scipy import sparse
from scipy.sparse import linalg
import time
import math
from scipy import stats
from sklearn.cluster import KMeans
from random import shuffle

import sys
sys.path += ['/Package'] ### Specify the directory where the Package is
from clustering_more import *


'''This function outputs the directed edge list representation given the adjacency matrix'''

def edge_list(A): # given the adjacency matrix it return the directed edge list representation
    
    size = sum(np.sum(A,axis = 0))  # number of (directed) edges in the network
    edge_list = [[] for x in range(size)]
    k = 0
    
    for i in range(len(A)):
        for j in range(i+1,len(A)):
            
            if A[i,j] == 1:
                edge_list[k] = [i,j]
                edge_list[k+1] = [j,i]
                k += 2
                
    return edge_list  


'''This function generates the non-backtracking matrix from the adjacency matrix'''

def non_backtraking(A):
    
    el =  edge_list(A)
    
    m = len(el) # m is twice the number of the total edges
    B = np.zeros((m,m)) 
    
    for i in range(m):
        for j in range(m):
            if el[i][1] == el[j][0] and el[i][0] != el[j][1]:
                B[i][j] = 1
                
    return B




n_clusters = 4
n = 5000

c_in = 20
c_out = 5



fluctuation = 1.5


# generation of the C\Pi matrix

fractions = np.ones(n_clusters)
fractions = fractions/sum(fractions)

C_matrix = np.random.uniform(c_out-fluctuation,c_out+fluctuation,size = (n_clusters,n_clusters))
i_lower = np.tril_indices(n_clusters, -1)
C_matrix[i_lower] = C_matrix.T[i_lower]
C_matrix[0][0] = c_in

Pi = np.diag(fractions)
M = np.dot(C_matrix,Pi)
c = np.sum(M, axis = 0)[0]
for i in range(1,len(C_matrix)):
    M[i][i] += c - np.sum(M, axis = 0)[i]
    C_matrix[i][i] = M[i][i]/fractions[i]
    
values = np.linalg.eigvals(M)
values = np.sort(values)

# degree distribution

theta = np.random.uniform(3,13,size = n)**4
theta = theta/np.mean(theta)

t0 = time.time()

#adjacency matrix

A = adj(C_matrix,theta,fractions)
print('done matrix A')
print('\ntime:' + str(time.time()-t0))
	  

# non-backtracking matrix

B = non_backtraking(A)
print('done matrix B')
print('\ntime:' + str(time.time()-t0))

# diagonalization

val = np.linalg.eigvals(B)
x = val.real
y = val.imag
	  
print('done')
print('\ntime:' + str(time.time()-t0))
