#!/usr/bin/env python
# coding: utf-8

# # TP - PageRank
# 
# Le but de cette séance est d’implémenter l’algorithme de puissance itérée 
# pour pouvoir ensuite calculer la centralité de vecteur propre dans le cas 
# des graphes non orientés et également dans le des graphes orientés (Page Rank).
# Nous utiliserons le module `networkx` de Python.
# 
# Pour l'ensemble des graphes que nous allons manipuler dans ce TP, il faut 
# que les sommets soient obligatoirement numéroté de 0 à n-1.


import networkx as nx
import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
import time
import math


# ## Quelques fonctions utiles

# Return the adjacency matrix of G as list of lists
def adjacencyMatrix(G):
    return nx.to_numpy_array(G, nodelist=range(len(G))) 

# Return the matrix product of A and B
def multMatrix(A,B):
    return A.dot(B)

# Return the same graph as G whose vertices are relabeled from 0 to n-1
# and the mapping (old:new) names
def relabelGraph(G):
    mapping = dict(zip(G, range(0, len(G))))
    return nx.relabel_nodes(G, mapping), mapping 

# Return original name of the vertex from new name
def getOriginalName(mapping, v):
    return [k for k,u in mapping.items() if u == v][0]

# Draw the graph G. If a centrality measure is given in C, label the vertices. 
# Different layouts can be used: neato (default), fdp, patchwork, circo, ... 
# List of layout : https://graphviz.org/docs/layouts/
def drawGraph(G,C=None,prog="neato"):
    pos=nx.nx_agraph.pygraphviz_layout(G,prog=prog) # Position of vertices
    node_size=[300 for _ in G] # default size of nodes
    with_labels=True
    if C is not None: # if a centrality measure is given, label the vertices
        C = C.round(6)
        with_labels=False
        labels_nodes = {v:C[v][0] for v in G}
        max_C=max(C)[0]
        node_size=[math.exp(labels_nodes[v]*7/max_C) for v in labels_nodes] # size the vertices depending on the centrality measure
        nx.draw_networkx_labels(G,pos,labels_nodes,font_size=12,font_color="black") # draw labels
    nx.draw(G,pos=pos,node_size=node_size,with_labels = with_labels) # draw graph
    plt.show() 


# ## Quelques graphes pour faire des tests
# 
# 1. Un premier exemple de graphe non orienté à 4 sommets :

G1=nx.Graph()
G1.add_edge(0,1); G1.add_edge(0,3); G1.add_edge(1,2); G1.add_edge(2,0)
print("Matrice d'adjacence du graphe G1: \n", adjacencyMatrix(G1))
drawGraph(G1)

# 2. Un deuxième exemple de graphe généré aléatoirement sur le modèle Erdös-Rényi :
# 
# Rappel : Le graphe aléatoire G_{n,p} construit sur le modèle Erdos-Rényi est un 
# graphe à n sommets dans lequel il existe une arête entre chaque paire de sommets 
# avec une probabilité p.
# 
# Dans l'exemple suivant : n=20 et p=1/5.
# 
 
n=20; p=1/5
oriented=False
G2 = nx.generators.random_graphs.fast_gnp_random_graph(n,p,directed=oriented)
drawGraph(G2,prog="neato")

# 3. Un plus gros graphe : `Gcit`
# 
# Il s'agit d'un graphe de citations qui contient 27770 sommets et de 352807 arcs. 
# Chaque sommet représente un article scientifique et les arcs représentent les
# citations : si l'article A cite l'article B, alors il y aura un arc de A vers B.

file_name="cit-HepTh.txt"
Gcit = nx.DiGraph()
with open(file_name, 'r') as data:
        for line in data:
            p = line.split(" ")       
            Gcit.add_edge(int(p[0]),int(p[1]))
Gcit,mappingGcit = relabelGraph(Gcit) 

# # Exercice 1 : Implémentation de l'algorithme de puissance itérée
# 
# 1. Complétez la fonction `powerIteration` ci-dessous. Cette fonction
# implémente l’algorithme de puissance itérée ; elle prend
# en paramètre une matrice A, un nombre d'itérations maximum et un réel
# epsilon permettant de décider de l’arrêt. Cette fonction retourne un
# vecteur-colonne.
# 
# À chaque étape de calcul (y compris à l'étape initiale), faites en
# sorte que la somme des valeurs du vecteur-colonne soit égale à 1.
# 
# Pour info : un vecteur-colonne peut être obtenu avec
# `np.ndarray(shape=(n,1))` où n est la taille du vecteur.

def powerIteration(A, max_iterations=100, epsilon=1e-6):
    n = A.shape[0] # nombre de lignes de A = nombre de sommets du graphe

    # 1. Initialiser un vecteur-colonne C de taille n avec des valeurs uniformes (1/n)
    C = # TODO

    for i in range(max_iterations):
        # 2. Conserver une copie du vecteur C de l'itération précédente.
        C_prec = # TODO

        # 3. Appliquer une itération de la méthode de la puissance
        C = # TODO

        # 4. Normaliser le vecteur C pour que la somme de ses composantes soit égale à 1.
        # TODO
        
        # 5. Calculer la différence maximale en valeur absolue entre C et C_prec.
        # Si cette différence est inférieure à epsilon, arrêter la boucle.
        if # TODO :
            return C
        
    return C

# 2. *Normalisation* : Pourquoi l'étape de normalisation
# à /chaque/ itération est-elle cruciale ? Que risquerait-il de se
# passer si on ne normalisait pas, ou seulement à la fin ?

# TODO

# 3. *Critère d'arrêt* : Vous avez utilisé la plus grande différence en
# valeur absolue entre les composantes de C et C_prec. Une autre
# approche consiste à utiliser la somme des différences.
# Selon vous, lequel de ces deux critères est le plus
# strict ? Justifiez votre réponse.

# TODO

# 4. *Test* : Testez votre fonction sur une matrice simple comme A =
# np.array([[0.5, 0.5], [0.5, 0.5]]). Quel est le vecteur propre
# attendu ? Votre fonction le trouve-t-elle ?

# TODO

# # Exercice 2 : Centralité de vecteur propre (cas non orienté)

# 1. *Testez* : Utilisez la fonction `powerIteration` écrite à l'exercice 1 pour calculer la centralité de 
# vecteur propre **dans le cas non orienté** (utilisez le graphe `G1` défini plus haut). Vous 
# fixerez un réel epsilon permettant d’indiquer la précision du calcul. Faites appel à la 
# fonction `drawGraph` pour voir le résultat de vos calculs s’afficher sur les sommets du graphe.

# TODO

# 2. *Expérimentation* : Modifiez votre fonction powerItération pour qu'elle
# retourne aussi le nombre d'itérations. Faites varier ~epsilon~ (1e-3,
# 1e-6, 1e-9) et observez l'impact sur le nombre d'itérations. Quel
# compromis cela met-il en évidence ?

# TODO

# # Exercice 3 : Page Rank (cas orienté)
# 
# 1. Écrire une fonction `matrixPageRank` permettant de construire la matrice de passage du 
# Page Rank qui nous permettra de calculer la centralité de vecteur propre dans 
# le **cas orienté** . Cette fonction prend en paramètre un graphe orienté et un réel alpha 
# correspondant au facteur de zap (dumping factor). 

def matrixPageRank(G,alpha=0.85):
  # TODO: Implémenter la construction de la matrice P comme vu en cours
    pass  

# 2. Utilisez la fonction `matrixPageRank` pour calculer la matrice de passage du graphe 
# disponible sur Moodle (que vous devez construire auparavant), puis calculez la valeur du 
# Page Rank de chaque sommet à l'aide de la fonction `powerIteration`.  
# 
# La valeur du Page Rank de chaque sommet est la suivante (arrondie à 2 décimales) :
# 
# `[[0.31]
#  [0.21]
#  [0.26]
#  [0.21]]` 

# TODO

# 3. Testez votre fonction sur le graphe de citations scientifiques
# `GCit`. Quel noeud est le plus important ?

# TODO

# 4. *Rôle d'alpha* : Recalculez le PageRank de G3 pour alpha = 1 et
# alpha = 0. Que se passe-t-il pour ~alpha = 1~ ? Le résultat
# est-il celui attendu ? Que représente le PageRank quand alpha =
# 0 ? Le résultat est-il logique ? 

# TODO

# ## Exercice 4
# 
# La construction de la matrice de passage P dans matrixPageRank
# est coûteuse, notamment à cause des produits matriciels. L'objectif
# est de construire cette matrice P de manière plus efficace, en
# travaillant ligne par ligne.

# Pour un sommet v donné, la ligne v de la matrice P dépend de deux cas :
# 1. Si v n'est pas un puits (d+(v) > 0) : la ligne v de P est
# alpha * (1/d+(v)) * A_v + (1-alpha)/n * u', où A_v est la ligne
# v de la matrice d'adjacence et u' un vecteur ligne de 1.
# 2. Si v est un puits (d+(v) = 0) : la ligne v de P est 1/n * u'.

# 1. Écrivez une nouvelle fonction matrixPageRankFast(G, alpha=0.85)
# qui implémente cette approche.

def matrixPageRankFast(G,alpha=0.85):
    # TODO: Implémenter la construction optimisée de la matrice P
    pass


# 2. *Validation* : Utilisez np.allclose(P_slow, P_fast) pour
# vérifier que votre nouvelle fonction produit bien la même matrice
# que matrixPageRank sur un petit graphe comme G3.

# TODO

# 3. *Performance* : Utilisez le module time pour comparer le temps
# d'exécution de matrixPageRank(Gcit) et
# matrixPageRankFast(Gcit). Calculez le gain de performance (ex: "la
# version rapide est X fois plus rapide").

# TODO