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 mpl_toolkits.mplot3d import Axes3D
from random import shuffle

import sys 
sys.path += ['/Package'] # Put here your directory
from clustering_more import *


n = 5000 # size of the network
n_blocks_c = 8 # number of points on the k axis
n_blocks_n = 8 # number of points on the alpha axis
n_average = 10 # number of averages

number_clusters_vector = np.arange(2,2+n_blocks_n)
n_blocks_n = len(number_clusters_vector)
c_out = 3
c_in_v = np.linspace(4,40,n_blocks_c) # c_in vector

n_cycles = 8
max_n_classes = max(number_clusters_vector)
percentage = 'no'

# degree distribution

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



t0 = time.time()

Over = np.zeros((n_blocks_c,n_blocks_n))
Over_s = np.zeros((n_blocks_c,n_blocks_n))

Detected = np.zeros((n_blocks_c,n_blocks_n)) # number of classes detected according to our algorithm
Detectable = np.zeros((n_blocks_c,n_blocks_n)) # number of classes theoretically detectable

phi = np.mean(theta**2)

for i in range(n_blocks_c):
    
    c_in = c_in_v[i]

    for j in range(n_blocks_n):

        n_clusters = number_clusters_vector[j]
       
        ### generation of the adjacency matrix

        fluctuation = 2/n_clusters
        fractions = np.ones(n_clusters)/n_clusters

        C_matrix = np.random.uniform(c_out-fluctuation,c_out+fluctuation,size = (n_clusters,n_clusters)) # set the off diagonal elements
        i_lower = np.tril_indices(n_clusters, -1)
        C_matrix[i_lower] = C_matrix.T[i_lower] # symmetrize the matrix
        C_matrix[0][0] = c_in # fix C_11

        Pi = np.diag(fractions)
        M = np.dot(C_matrix,Pi) # M = C\Pi
        c = np.sum(M, axis = 1)[0]
        for l in range(1,len(C_matrix)):
            M[l][l] += c - np.sum(M, axis = 1)[l]
            C_matrix[l][l] = M[l][l]/fractions[l]

        
        val, vec = np.linalg.eig(M)
        detectable = np.sum((val > np.sqrt(c/phi))*1) # number of detecteble classes
        v,vv = np.linalg.eig(M)
        v = np.sort(v)
        r_vector = v[-1]/v[::-1]
        r_vector = r_vector[1:]


        for k in range(n_average):
            
            new_n = int(n/n_clusters)*n_clusters # change the value of n to the closest integer such that n/n_clusters is still an integer
            theta_new = theta[:new_n]
            real_classes = np.zeros(new_n)
            for m in range(n_clusters):
                real_classes[int(new_n*np.sum(fractions[:m])):int(new_n*np.sum(fractions[:m+1]))] = m

            real_classes = real_classes.astype(int)
            real_classes[-1] = m
            
            A = adj(C_matrix,theta_new,fractions)  # generation of the adjacency matrix
            A = sparse.coo_matrix(A)
            A = A.astype(float)
            classes, X, nc  = BH_with_r(A,n_cycles, n_clusters, n_clusters, percentage, r_vector)  # proposed algorithm
            classes_s, X_s  = saade(A, n_clusters) # saade's algorithm

            classes, ov = overlap(real_classes, classes)
            classes_s, ov_s = overlap(real_classes, classes_s)

            Over[i][j] += ov/n_average
            Over_s[i][j] += ov_s/n_average

            d = np.sum(A.A, axis = 0)
            r_max = np.sqrt(np.sum(d**2)/np.sum(d))
            r = r_max

            eigenvalues, eigenvectors = Hessian(A,d,r,n_clusters)
            n_detected = sum((eigenvalues < 0)*1)

            Detected[i][j] += n_detected/n_average
            Detectable[i][j] += detectable/n_average

            OUT = 'Completion : i = ' + str(i+1) + '/'+str(n_blocks_c) + ', j = ' + str(j+1) + '/'+str(n_blocks_n) +', k = ' + str(k+1) + '/'+str(n_average) + ' Time : ' + str(time.time()-t0)
            sys.stdout.write('\r%s' % OUT)
            sys.stdout.flush()

