# PageRank Algorithm

This notebook implements the PageRank algorithm, prepared as a homework in BLG202E - Numerical Methods in CE class at ITU, Spring 2020.


```python
!pip install mechanize
```

    Defaulting to user installation because normal site-packages is not writeable
    Requirement already satisfied: mechanize in /home/marif/.local/lib/python3.6/site-packages (0.4.5)
    Requirement already satisfied: html5lib>=0.999999999 in /home/marif/.local/lib/python3.6/site-packages (from mechanize) (1.0.1)
    Requirement already satisfied: six>=1.9 in /home/marif/.local/lib/python3.6/site-packages (from html5lib>=0.999999999->mechanize) (1.14.0)
    Requirement already satisfied: webencodings in /home/marif/.local/lib/python3.6/site-packages (from html5lib>=0.999999999->mechanize) (0.5.1)



```python
import numpy as np #For arrays and linear algebra routines
from mechanize import Browser #For getting the names of the papers from arXiv in part 3 results
                              #Not necessary, but to use it, above cell must be run to install mechanize library 
```


```python
eps = 1e-10
```

# Motivating Example: Academic Ranking System

Our motivating example uses the paper citation network of Arxiv High Energy Physics Theory category. We will rank the papers by their influence. Data is represented as a directed graph, and we have the list of edges in a text file. We will start by loading the text file that contains the directed graph into a NumPy array.


```python
array = np.genfromtxt("cit-HepTh.txt", dtype=int)
```

Let's see the first 10 edges.


```python
print(array[0:10])
```

    [[   1001 9304045]
     [   1001 9308122]
     [   1001 9309097]
     [   1001 9311042]
     [   1001 9401139]
     [   1001 9404151]
     [   1001 9407087]
     [   1001 9408099]
     [   1001 9501030]
     [   1001 9503124]]


We will create a transition matrix from the list that contains links between the nodes of the graph. A transition matrix shows transitions from nodes to other nodes and the probability that transition happens.

These weights or normalization factors are calculated by taking $\frac{1}{\textrm{(no of outgoing links from that node)}}$. Following code will define a function with parameter **arr**: multidimensional array of edges, and returns **labels**: array consisting of distinct nodes and **weights**: a dictionary that maps these nodes to their weights.


```python
def createWeights(arr):
    labels, counts = np.unique(arr[:,0], return_counts=True) ##Finds the unique entries in first column, returns their values 
                                                              #and their counts to calculate weights
    weight_dict = dict()
    for index, label in enumerate(labels):                   ##Creates a dictionary that holds the normalization factor,
        weight_dict[label] = 1 / counts[index]               # 1/(number of outgoing links) for every paper ID in labels

    return labels, weight_dict
```


```python
labels, weights = createWeights(array)
```

Following code will initialize transition matrix, finds the cells that represent transitions using a for loop that examines each node and links from that node, and fills in the cell with the weight specified in weights dictionary created above.


```python
def createTransitionMatrix(arr, labels, weight_dict):
    transitionMatrix = np.zeros((labels.size,labels.size))
    
    for i, label in enumerate(labels):
        links_from_label = arr[np.nonzero(label == arr[:,0]),1][0]
      #  if links_from_label.shape == (0,):                          ##For dangling nodes, nodes that have no outgoing links
      #      weight_dict[label] = 1 / labels.size                     
      #      links_to_label = labels                                  
        transitionMatrix[np.searchsorted(labels, links_from_label),i] = weight_dict[label]
    
    return transitionMatrix
```


```python
rankMatrix = createTransitionMatrix(array, labels, weights)
```

We expect our transition matrix to be column-stochastic, and we also expect it to be a sparse matrix it is mostly populated by zeros. Operations on sparse matrices can be done by some faster methods, so sparsity is an advantage in speed.

Now we will define and run two functions to test these attributes of our matrix. Results will show that our expectations hold.


```python
def checkStochastic(matrix):
    eps = 1e-10
    print("Sums of columns that do not add to 1 in a reasonable error margin {} will be shown.".format(eps))
    for i in matrix.sum(0):
        if not (i - 1 < eps):
            print(i)
    print("Calculation finished, average value of sums of all columns is {}.".format(np.mean(matrix.sum(0))))

def checkSparsity(matrix):
    zeros = np.count_nonzero(matrix==0)
    elements = matrix.size
    sparse_rate = np.divide(zeros,elements)
    print("There are {} elements in total, {} of them are zero.".format(elements, zeros))
    print("Sparsity rate of this matrix is %{}".format(sparse_rate*100))
    return sparse_rate
```


```python
checkStochastic(rankMatrix)
sparsity = checkSparsity(rankMatrix)
```

    Sums of columns that do not add to 1 in a reasonable error margin 1e-10 will be shown.
    Calculation finished, average value of sums of all columns is 0.9958169865008099.
    There are 627953481 elements in total, 627601470 of them are zero.
    Sparsity rate of this matrix is %99.94394314059069



---
Now that we have the required matrix, we can solve the equation
\begin{equation}
A x = x
\end{equation}
where A is the matrix, and x is the result vector that contains the rank. We will solve this by **power method**, by repeatedly multiplying an arbitrary vector* by our matrix until the difference in resulting vectors of two iterations is smaller than epsilon.

*While any arbitrary vector should work, it is better practice to use an all ones vector normalized by the size of itself, so initially every rank is equal and vector sums up to one. We will follow this practice.


```python
def solveRank(rankMatrix):
    eps = 1e-7
    v0 = np.ones(rankMatrix.shape[0]) / rankMatrix.shape[0]
#   v0 = np.random.random(rankMatrix.shape[0] / rankMatrix.shape[0])

    counter = 1
    while True:
        v = np.dot(rankMatrix, v0)
        v = v / np.linalg.norm(v)
        if (np.mean(np.abs(np.subtract(v,v0))) < 2*eps):
            break
        
#        print("Error: {}".format(np.mean(np.abs(np.subtract(v,v0)))))   ##Uncomment this line to print error in each step
                                                                          #If this function is taking too long, printing the error may be a good idea for debugging
        counter += 1
        v0 = v
    
    print("Appropriate vector found in {} iterations, final difference between two iteration result vectors was less than {}.".format(counter,eps))
    return v
```


```python
final = solveRank(rankMatrix)
```

    Appropriate vector found in 80 iterations, final difference between two iteration result vectors was less than 1e-07.


Finally, we will rank our papers according to the resulting vector, and show the first 10 papers.


```python
def rankPagesDescending(labels, final):
    return labels[final.argsort()][::-1]
```


```python
rankPagesDescending(labels, final)[0:10]
```




    array([9201015, 9207016, 9206003,  209015, 9205071, 9202067, 9201047,
           9205038, 9202018, 9205006])



Using the mechanize library, we can collect the information on these papers.


```python
ranking = rankPagesDescending(labels, final)
br = Browser()

for index,paper_id in enumerate(ranking[0:10]):
    str_id = str(paper_id)
    page_url = "https://arxiv.org/abs/hep-th/"
    while(len(str_id) < 7):
        str_id = '0' + str_id
    page_url += str_id
    br.open(page_url)
    paper_title = br.title()[17:]
    print("{}. paper ID is {}.".format(index + 1, paper_id))
    print("Name of the paper is: {}".format(paper_title))
    print(page_url)
```

    1. paper ID is 9201015.
    Name of the paper is: An Algorithm to Generate Classical Solutions for String Effective Action
    https://arxiv.org/abs/hep-th/9201015
    2. paper ID is 9207016.
    Name of the paper is: Noncompact Symmetries in String Theory
    https://arxiv.org/abs/hep-th/9207016
    3. paper ID is 9206003.
    Name of the paper is: From Form Factors to Correlation Functions: The Ising Model
    https://arxiv.org/abs/hep-th/9206003
    4. paper ID is 209015.
    Name of the paper is: Advances in String Theory in Curved Space Times
    https://arxiv.org/abs/hep-th/0209015
    5. paper ID is 9205071.
    Name of the paper is: Novel Symmetry of Non-Einsteinian Gravity in Two Dimensions
    https://arxiv.org/abs/hep-th/9205071
    6. paper ID is 9202067.
    Name of the paper is: Stringy Domain Walls and Other Stringy Topological Defects
    https://arxiv.org/abs/hep-th/9202067
    7. paper ID is 9201047.
    Name of the paper is: Duality-Invariant Gaugino Condensation and One-Loop Corrected Kahler Potentials in String Theory
    https://arxiv.org/abs/hep-th/9201047
    8. paper ID is 9205038.
    Name of the paper is: Recent Developments in Classical and Quantum Theories of Connections, Including General Relativity
    https://arxiv.org/abs/hep-th/9205038
    9. paper ID is 9202018.
    Name of the paper is: Jones Polynomials for Intersecting Knots as Physical States of Quantum Gravity
    https://arxiv.org/abs/hep-th/9202018
    10. paper ID is 9205006.
    Name of the paper is: Stabilized Quantum Gravity: Stochastic Interpretation and Numerical Simulation
    https://arxiv.org/abs/hep-th/9205006


In some edge cases, the transition matrix above may not abide by the conditions algorithm requires. To account for those situations, there is a **damping factor** defined for the algorithm, which takes the weighted average of our transition matrix with a matrix of all ones. While transition matrix represents probabilities of going from one node to another, adding these damping factor will give us a chance to randomly go from any node to another.


```python
def addDamping(transitionMatrix, labels, p = 0.15):
    if(p < 0 or p > 1):
        print("Please try again with a damping factor in interval [0,1].")
        return None
    rankMatrix = (1 - p) * transitionMatrix + p * ((1/labels.size) * np.ones((labels.size,labels.size)))
    return rankMatrix
```


```python
rankMatrix = addDamping(rankMatrix, labels)
checkStochastic(rankMatrix)
sparsity = checkSparsity(rankMatrix)
```

    Sums of columns that do not add to 1 in a reasonable error margin 1e-10 will be shown.
    Calculation finished, average value of sums of all columns is 0.9964444385253646.
    There are 627953481 elements in total, 0 of them are zero.
    Sparsity rate of this matrix is %0.0


We see that while stochastic structure of the matrix holds, it is no longer a sparse matrix after adding damping factor. This is because now every node has the probability to go to any other node, even if this is a very small probability. So we do not have any 0 cells in our matrix anymore.

Now, we will solve out matrix again.


```python
final = solveRank(rankMatrix)
ranking = rankPagesDescending(labels, final)
br = Browser()

for index,paper_id in enumerate(ranking[0:10]):
    str_id = str(paper_id)
    page_url = "https://arxiv.org/abs/hep-th/"
    while(len(str_id) < 7):
        str_id = '0' + str_id
    page_url += str_id
    br.open(page_url)
    paper_title = br.title()[17:]
    print("{}. paper ID is {}.".format(index + 1, paper_id))
    print("Name of the paper is: {}".format(paper_title))
    print(page_url)
```

    Appropriate vector found in 30 iterations, final difference between two iteration result vectors was less than 1e-07.
    1. paper ID is 9201015.
    Name of the paper is: An Algorithm to Generate Classical Solutions for String Effective Action
    https://arxiv.org/abs/hep-th/9201015
    2. paper ID is 9207016.
    Name of the paper is: Noncompact Symmetries in String Theory
    https://arxiv.org/abs/hep-th/9207016
    3. paper ID is 9205071.
    Name of the paper is: Novel Symmetry of Non-Einsteinian Gravity in Two Dimensions
    https://arxiv.org/abs/hep-th/9205071
    4. paper ID is 209015.
    Name of the paper is: Advances in String Theory in Curved Space Times
    https://arxiv.org/abs/hep-th/0209015
    5. paper ID is 9202067.
    Name of the paper is: Stringy Domain Walls and Other Stringy Topological Defects
    https://arxiv.org/abs/hep-th/9202067
    6. paper ID is 9201047.
    Name of the paper is: Duality-Invariant Gaugino Condensation and One-Loop Corrected Kahler Potentials in String Theory
    https://arxiv.org/abs/hep-th/9201047
    7. paper ID is 9202018.
    Name of the paper is: Jones Polynomials for Intersecting Knots as Physical States of Quantum Gravity
    https://arxiv.org/abs/hep-th/9202018
    8. paper ID is 9205006.
    Name of the paper is: Stabilized Quantum Gravity: Stochastic Interpretation and Numerical Simulation
    https://arxiv.org/abs/hep-th/9205006
    9. paper ID is 9205038.
    Name of the paper is: Recent Developments in Classical and Quantum Theories of Connections, Including General Relativity
    https://arxiv.org/abs/hep-th/9205038
    10. paper ID is 9206048.
    Name of the paper is: Conformally Exact Results for SL(2,R)\times SO(1,1)^{d-2}/SO(1,1) Coset Models
    https://arxiv.org/abs/hep-th/9206048


---

Since our dataset is very big, it is not easy to say if our algorithm correctly chooses the rankings. So we will test the same functions on a smaller dataset below, and we will compare the results with our expectation that will be computed by hand mathematically.



```python
def calculateRankVector(data, damping=False):
    labels, weights = createWeights(data)
    rankMatrix = createTransitionMatrix(data, labels, weights)
    if damping:
        rankMatrix = addDamping(rankMatrix, labels)
    resultVector = solveRank(rankMatrix)
    return resultVector
```


```python
filename = "testset.txt" #Change the filename to test it with another file 
data = np.genfromtxt(filename, dtype = str)
resultVector = calculateRankVector(data, damping=True)
ranking = rankPagesDescending(labels, resultVector)
print(ranking)
print(resultVector.reshape(resultVector.size,1))
```

    Appropriate vector found in 20 iterations, final difference between two iteration result vectors was less than 1e-07.
    ['A' 'D' 'C' 'B']
    [[0.63395272]
     [0.34267715]
     [0.444879  ]
     [0.53175087]]


a9535e3ff61bf761c668565600181c90ec44bfef
ipynb
Jupyter Notebook
pagerank.ipynb
marifdemirtas/pagerank
0d5796c5720a35aa84b2aa1ef98343a28016d390
MIT
pagerank.ipynb
marifdemirtas/pagerank
0d5796c5720a35aa84b2aa1ef98343a28016d390
MIT
pagerank.ipynb
marifdemirtas/pagerank
0d5796c5720a35aa84b2aa1ef98343a28016d390
MIT
Qwen/Qwen-72B
1. YES
2. YES