#! /usr/bin/env python3
#
def digraph_plot ( A, header ):
#*****************************************************************************80
#
## digraph_plot() plots a directed graph.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 August 2022
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer A(N,N): the incidence matrix.
#
# string header: the prefix to be used for the comment,
# dot filename, and png filename.
#
from graphviz import Digraph
dot = Digraph ( comment = header, format = 'png' )
#
# Specify the nodes, giving each an internal code, and a label.
#
n = A.shape[0]
for i in range ( 0, n ):
code = str ( i )
label = str ( i )
dot.node ( code, label )
#
# Specify the edges as connections between two nodes.
#
for i in range ( 0, n ):
for j in range ( 0, n ):
if ( A[i,j] ):
dot.edge ( str ( i ), str ( j ) )
print ( dot.source )
#
# Save graph to a file, and optionally display an image to the screen.
#
filename = header + '.dot'
dot.render ( filename, view = False )
filename = header + '.png'
print ( ' Graphics saved as "' + filename + '"' )
return
def google_matrix ( A, p = 0.15 ):
#*****************************************************************************80
#
## google_matrix() converts an incidence matrix to a Google transition matrix.
#
# Discussion:
#
# If the incidence matrix has a node I with no connectivity,
# (A(I,1:N) = 0) then we artificially set G(1:N,I)=1/N, so that
# the transition matrix property of having unit column sums is preserved.
#
# For nodes with connectivity, we assume that 85 percent of the time
# we will take a link at random, and 15 percent of the time, we will
# jump to an arbitrary link.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 August 2022
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer A(N,N): the incidence matrix.
#
# real P: the probability that the user will simply jump to a random
# web page. Default = 0.15
#
# Output:
#
# real G(N,N): the Google matrix.
#
import numpy as np
n = A.shape[0]
#
# Compute the row sums of A.
# s(i) is the number of links on web page i.
#
s = np.sum ( A, axis = 1 )
#
# Allocate G.
#
G = np.zeros ( [ n, n ] )
#
# Fill any zero row with 1/n.
#
for i in range ( 0, n ):
if ( s[i] == 0.0 ):
G[i,:] = 1.0 / n
else:
G[i,:] = ( 1.0 - p ) * A[i,:] / s[i] + p / n
#
# Transpose so the matrix is a transition/stochastic matrix.
#
G = np.transpose ( G )
return G
def google_rank ( A, it_num ):
#*****************************************************************************80
#
## google_rank() uses the power method on the Google matrix.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 August 2022
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer A(N,N): the incidence matrix.
#
# integer it_num: the number of iterations to take.
#
import numpy as np
print ( '' )
print ( 'google_rank():' )
print ( ' Given an NxN incidence matrix A,' )
print ( ' compute the Google matrix G,' )
print ( ' Then start with a vector of N values 1/N,' )
print ( ' and repeatedly compute x <= G*x' )
print ( '' )
print ( ' After many steps, compare last three iterates.' )
print ( ' If they are close, we are probably at an eigenvector' )
print ( ' associated with the eigenvalue 1.' )
n = A.shape[0]
#
# Compute the Google matrix.
#
G = incidence_to_google ( A )
#
# Carry out many iterations.
# Normalization is not necessary because G is a transition matrix.
#
for i in range ( 0, it_num + 1 ):
if ( i == 0 ):
x = np.ones ( n ) / n
else:
x = np.matmul ( G, x )
#
# Compare three successive iterates.
#
Gx = np.matmul ( G, x )
GGx = np.matmul ( G, Gx )
Output = np.column_stack ( [ x, Gx, GGx ] )
print ( '' )
print ( ' x, G*x, G*G*x' )
print ( Output )
return
def incidence_to_google ( A ):
#*****************************************************************************80
#
## incidence_to_google() converts an incidence matrix to a Google transition matrix.
#
# Discussion:
#
# If the input incidence matrix has a node I with no connectivity,
# (A(I,1:N) = 0) then we artificially set T(1:N,I)=1/N, so that
# the transition matrix property of having unit column sums is preserved.
#
# For nodes with connectivity, we assume that 85 percent of the time
# we will take a link at random, and 15 percent of the time, we will
# jump to an arbitrary link.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 August 2022
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer A(N,N): the incidence matrix.
#
# Output:
#
# real T(N,N): the transition matrix.
#
import numpy as np
n = A.shape[0]
T = np.zeros ( [ n, n ] )
s = np.sum ( A, axis = 1 )
for i in range ( 0, n ):
if ( s[i] == 0.0 ):
T[i,:] = 1.0 / n
else:
T[i,:] = 0.85 * A[i,:] / s[i] + 0.15 / n
T = np.transpose ( T )
return T
def incidence_to_transition ( A ):
#*****************************************************************************80
#
## incidence_to_transition() converts an incidence matrix to a transition matrix.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 August 2022
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer A(N,N), the incidence matrix.
#
# Output:
#
# real T(N,N): the transition matrix.
#
import numpy as np
n = A.shape[0]
#
# Get the row sums.
#
s = np.sum ( A, axis = 1 )
#
# Normalize each row so it sums to 1.
#
T = np.zeros ( [ n, n ] )
for i in range ( 0, n ):
if ( s[i] != 0.0 ):
T[i,:] = A[i,:] / s[i]
else:
T[i,i] = 1.0
#
# Transpose the matrix.
#
T = np.transpose ( T )
return T
def mac1_inc ( ):
#*****************************************************************************80
#
## mac1_inc() returns the incidence matrix for MacCormick's example 1.
#
# Discussion:
#
# This matrix appears on page 35 of the reference.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 August 2022
#
# Author:
#
# John Burkardt
#
# Reference:
#
# John MacCormick,
# Nine Algorithms That Changed the Future: The Ingenious Ideas that
# Drive Today's Computers,
# Princeton University Press,
# ISBN-13: 978-0691158198.
#
# Output:
#
# integer A(5,5), the incidence matrix.
#
import numpy as np
A = np.zeros ( [ 5, 5 ] )
A[0,1] = 1.0
A[1,4] = 1.0
A[2,0] = 1.0
A[3,0] = 1.0
A[4,0] = 1.0
return A
def mac2_inc ( ):
#*****************************************************************************80
#
## mac_inc2() returns the incidence matrix for MacCormick's example 2.
#
# Discussion:
#
# This matrix appears on page 33 of the reference.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 August 2022
#
# Author:
#
# John Burkardt
#
# Reference:
#
# John MacCormick,
# Nine Algorithms That Changed the Future: The Ingenious Ideas that
# Drive Today's Computers,
# Princeton University Press,
# ISBN-13: 978-0691158198.
#
# Output:
#
# integer A(16,16), the incidence matrix.
#
import numpy as np
A = np.zeros ( [ 16, 16 ] )
A[0,1] = 1.0
A[0,2] = 1.0
A[0,5] = 1.0
A[1,3] = 1.0
A[2,3] = 1.0
A[2,4] = 1.0
A[3,4] = 1.0
A[4,6] = 1.0
A[4,7] = 1.0
A[5,4] = 1.0
A[5,6] = 1.0
A[6,14] = 1.0
A[7,6] = 1.0
A[7,8] = 1.0
A[7,9] = 1.0
A[8,9] = 1.0
A[9,6] = 1.0
A[9,10] = 1.0
A[9,11] = 1.0
A[9,12] = 1.0
A[9,13] = 1.0
A[10,14] = 1.0
A[11,14] = 1.0
A[12,14] = 1.0
A[13,14] = 1.0
A[14,15] = 1.0
A[15,0] = 1.0
return A
def moler_inc ( ):
#*****************************************************************************80
#
## moler_inc() returns the incidence matrix for a network of Moler.
#
# Discussion:
#
# This matrix appears on page 5 of the reference.
#
# We added a self-link for node 5, which otherwise has no outlinks.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 August 2022
#
# Author:
#
# John Burkardt
#
# Reference:
#
# Cleve Moler,
# Experiments with Matlab,
# Chapter 7: Google PageRank,
# https://www.mathworks.com/moler/exm/chapters/pagerank.pdf
#
# Output:
#
# integer A(6,6), the incidence matrix.
#
import numpy as np
A = np.zeros ( [ 6, 6 ] )
A[0,1] = 1.0
A[0,5] = 1.0
A[1,2] = 1.0
A[1,3] = 1.0
A[2,3] = 1.0
A[2,4] = 1.0
A[2,5] = 1.0
A[3,0] = 1.0
A[4,4] = 1.0
A[5,0] = 1.0
return A
def pagerank_test ( ):
#*****************************************************************************80
#
## pagerank_test() tests pagerank().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 August 2022
#
# Author:
#
# John Burkardt
#
import platform
print ( '' )
print ( 'pagerank_test():' )
print ( ' Python version: ' + platform.python_version ( ) )
print ( ' pagerank() demonstrates several page rank algorithms.' )
for test in range ( 0, 5 ):
print ( '' )
print ( ' Test', test )
if ( test == 0 ):
label = '5 node MacCormick example with a cycle'
A = mac1_inc ( )
header = 'mac1'
elif ( test == 1 ):
label = '16 node MacCormick example'
n = 16
A = mac2_inc ( )
header = 'mac2'
elif ( test == 2 ):
label = '6 node Moler example'
A = moler_inc ( )
header = 'moler'
elif ( test == 3 ):
label = '6 node example'
A = six_inc ( )
header = 'six'
elif ( test == 4 ):
label = '15 node Sauer example'
A = sauer_inc ( )
header = 'sauer'
print ( label )
print ( ' Incidence matrix A:' )
print ( A )
T = incidence_to_transition ( A )
print ( ' Transition matrix T:' )
print ( T )
#
# Plot the directed graph.
#
digraph_plot ( A, header )
#
# Compute the rankings.
#
it_num = 100
google_rank ( A, it_num )
it_num = 100
power_rank ( A, it_num )
it_num = 10000
surf_rank ( A, it_num )
#
# Terminate.
#
print ( '' )
print ( 'pagerank_test():' )
print ( ' Normal end of execution.' )
return
def power_rank ( A, it_num ):
#*****************************************************************************80
#
## power_rank() uses the power method to seek an eigenvector of a directed graph.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 August 2022
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer A(N,N): the incidence matrix.
#
# integer it_num: the number of iterations to take.
#
import numpy as np
print ( '' )
print ( 'power_rank():' )
print ( ' Given an NxN incidence matrix A,' )
print ( ' compute the transition matrix T,' )
print ( ' Then start with a vector of N values 1/N,' )
print ( ' and repeatedly compute x <= T*x' )
print ( '' )
print ( ' After many steps, compare last three iterates.' )
print ( ' If they are close, we are probably at an eigenvector' )
print ( ' associated with the eigenvalue 1.' )
n = A.shape[0]
#
# Compute the transition matrix.
#
T = incidence_to_transition ( A )
#
# Carry out many iterations.
# Normalization is not necessary because T is a transition matrix.
#
for i in range ( 0, it_num + 1 ):
if ( i == 0 ):
x = np.ones ( n ) / n
else:
x = np.matmul ( T, x )
#
# Compare three successive iterates.
#
Tx = np.matmul ( T, x )
TTx = np.matmul ( T, Tx )
Output = np.column_stack ( [ x, Tx, TTx ] )
print ( '' )
print ( ' x, T*x, T*T*x' )
print ( Output )
return
def sauer_inc ( ):
#*****************************************************************************80
#
## sauer_inc() returns the incidence matrix for the Sauer example.
#
# Discussion:
#
# This matrix appears on page 564 of the reference.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 August 2022
#
# Author:
#
# John Burkardt
#
# Reference:
#
# Timothy Sauer,
# Numerical Analysis,
# Pearson, 2006,
# ISBN: 0-321-26898-9,
# LC: QA297.S348
#
# Output:
#
# integer A(15,15): the incidence matrix.
#
import numpy as np
A = np.zeros ( [ 15, 15 ] )
A[0,1] = 1
A[0,8] = 1
A[1,2] = 1
A[1,4] = 1
A[1,6] = 1
A[2,1] = 1
A[2,5] = 1
A[2,7] = 1
A[3,2] = 1
A[3,11] = 1
A[4,0] = 1
A[4,9] = 1
A[5,9] = 1
A[5,10] = 1
A[6,9] = 1
A[6,10] = 1
A[7,3] = 1
A[7,10] = 1
A[8,4] = 1
A[8,5] = 1
A[8,9] = 1
A[9,12] = 1
A[10,14] = 1
A[11,6] = 1
A[11,7] = 1
A[11,10] = 1
A[12,8] = 1
A[12,13] = 1
A[13,9] = 1
A[13,10] = 1
A[13,12] = 1
A[13,14] = 1
A[14,11] = 1
A[14,13] = 1
return A
def six_inc ( ):
#*****************************************************************************80
#
## six_inc() returns the incidence matrix for a 6-node example.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 August 2022
#
# Author:
#
# John Burkardt
#
# Output:
#
# integer A(6,6), the incidence matrix.
#
import numpy as np
A = np.zeros ( [ 6, 6 ] )
A[0,1] = 1.0
A[1,2] = 1.0
A[1,4] = 1.0
A[2,0] = 1.0
A[2,3] = 1.0
A[3,0] = 1.0
A[4,1] = 1.0
A[4,5] = 1.0
A[5,3] = 1.0
return A
def surf_rank ( A, it_num ):
#*****************************************************************************80
#
## surf_rank() uses the "random surfer" method for node rankings of a digraph.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 07 August 2022
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer A(N,N): the incidence matrix.
#
# integer it_num: the number of iterations to take.
#
from numpy.random import default_rng
import numpy as np
rng = default_rng ( )
print ( '' )
print ( 'surf_rank():' )
print ( ' Given an NxN incidence matrix A,' )
print ( ' compute the Google matrix T,' )
print ( ' and then repeatedly visit nodes, keeping track of' )
print ( ' how often you visited.' )
print ( '' )
print ( ' 15 per cent of the time, simply take a random jump to a node.' )
print ( ' The rest of the time, follow a random link from the current node.' )
print ( '' )
print ( ' When done, the node weight is the number of visits' )
print ( ' normalized by the total number of visits.' )
n = A.shape[0]
v = np.zeros ( n )
#
# Get the transition matrix from the incidence matrix.
#
T = incidence_to_transition ( A )
#
# Carry out many moves.
#
for k in range ( 0, it_num ):
if ( k == 0 ):
j = rng.integers ( low = 0, high = n, size = 1, endpoint = False )
else:
p1 = rng.random ( )
q1 = 0.15
if ( p1 <= q1 ):
j = rng.integers ( low = 0, high = n, size = 1, endpoint = False )
else:
p2 = rng.random ( )
q2 = 0.0
for i in range ( 0, n ):
q2 = q2 + T[i,j]
if ( p2 <= q2 ):
j = i
break
v[j] = v[j] + 1
#
# Normalize the counts.
#
v = v / it_num
print ( '' )
print ( ' Node weights:' )
print ( v )
return
def timestamp ( ):
#*****************************************************************************80
#
## timestamp() prints the date as a timestamp.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 August 2019
#
# Author:
#
# John Burkardt
#
import time
t = time.time ( )
print ( time.ctime ( t ) )
return
if ( __name__ == '__main__' ):
timestamp ( )
pagerank_test ( )
timestamp ( )