#! /usr/bin/env python3
#
def pagerank_test( ):

  import numpy as np

  print ( '' )
  print ( 'pagerank_test():' )
#
#  The power method for a random symmetric matrix.
#
  print ( '' )
  print ( '  power method for a random symmetric matrix:' )
  B = np.random.rand ( 5, 5 )
  B = B + np.transpose ( B )
  x = np.random.rand ( 5 )
  for k in range ( 0, 11 ):
    x = np.dot ( B, x )
    eigen = rayleigh ( B, x )
    print ( k, eigen )
    x = x / np.max ( np.abs ( x ) )

  print ( '  Rayleigh estimate:', eigen )
  eigen_all = np.linalg.eigvals ( B )
  print ( '  All eigenvalues from linalg.eigvals():' )
  print ( eigen_all )
#
#  Define A, the incidence matrix.
#
  A = np.array ( [ \
    [ 0, 1, 0, 0, 0, 1, 0 ], \
    [ 0, 0, 1, 0, 0, 0, 0 ], \
    [ 1, 0, 0, 1, 0, 0, 0 ], \
    [ 1, 0, 0, 0, 0, 0, 0 ], \
    [ 0, 0, 1, 0, 0, 0, 0 ], \
    [ 0, 0, 0, 0, 0, 0, 0 ], \
    [ 0, 0, 0, 0, 0, 0, 0 ] ], dtype = int )

  node_num = A.shape[0]

  print ( '' )
  print ( '  Incidence matrix A:' )
#
#  Plot the transition matrix.
#
  print ( '' )
  print ( '  Compute the transition matrix T:' )

  T = incidence_to_transition ( A )
#
#  Plot the transition matrix with and without labels.
#
  print ( '' )
  print ( '  Plot the transition matrix T with and without labels:' )

  node_labels = [ '0', '1', '2', '3', '4', '5', '6' ]

  transition_plot ( T, show_edge_labels = False, show_node_labels = False, 
    edge_labels = [], node_labels = node_labels, filename = 'transition'  )

  transition_plot ( T, show_edge_labels = True, show_node_labels = False, 
    edge_labels = T, node_labels = node_labels, filename = 'transition_probability' )
#
#  Use power method on T.
#
  print ( '' )
  print ( '  power method for T:' )
  x = np.random.rand ( node_num )
  for k in range ( 0, 21 ):
    x = np.dot ( T, x )
    eigen = rayleigh ( T, x )
    print ( k, eigen )
    x = x / np.max ( np.abs ( x ) )

  print ( '' )
  print ( '  Estimated eigenvector:' )
  print ( x )

  print ( '' )
  print ( '  Rayleigh estimate for eigenvalue:', eigen )
  eigen_all = np.linalg.eigvals ( T )
  print ( '  All eigenvalues from np.linalg.eigvals():' )
  print ( eigen_all )

  node_labels = np.zeros ( node_num )
  for i in range ( 0, node_num ):
    node_labels[i] = "{:.2f}".format(x[i])
 
  transition_plot ( T, show_edge_labels = True, show_node_labels = True, 
    edge_labels = T, node_labels = node_labels, filename = 'transition_ranking' )
#
#  Compute S
#
  S = incidence_to_surfer ( A )

  print ( '' )
  print ( '  power method for S:' )
  x = np.random.rand ( node_num )
  for k in range ( 0, 21 ):
    x = np.dot ( S, x )
    eigen = rayleigh ( S, x )
    print ( k, eigen )
    x = x / np.max ( np.abs ( x ) )

  print ( '' )
  print ( '  Estimated eigenvector:' )
  print ( x )

  print ( '' )
  print ( '  Rayleigh estimate for eigenvalue:', eigen )
  eigen_all = np.linalg.eigvals ( S )
  print ( '  All eigenvalues from np.linalg.eigvals():' )
  print ( eigen_all )

  node_labels = np.zeros ( node_num )
  for i in range ( 0, node_num ):
    node_labels[i] = "{:.2f}".format(x[i])
 
  transition_plot2 ( S, show_edge_labels = False, show_node_labels = True, 
    edge_labels = S, node_labels = node_labels, filename = 'surfer_ranking' )
  return

def incidence_to_transition ( A ):

#*****************************************************************************80
#
## incidence_to_transition() computes a transition matrix.
#
#  Input:
#
#    integer A(N,N), the incidence matrix.
#
#  Output:
#
#    real T(N,N): the transition matrix.
#
  import numpy as np
  n = A.shape[0]
  s = np.sum ( A, axis = 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
  T = np.transpose ( T )
  return T

def incidence_to_surfer ( A ):

#*****************************************************************************80
#
## incidence_to_surfer() computes a surfer matrix.
#
#  Input:
#
#    integer A(N,N), the incidence matrix.
#
#  Output:
#
#    real S(N,N): the surfer matrix.
#
  import numpy as np
  p = 0.15
  n = A.shape[0]
  s = np.sum ( A, axis = 1 )
  S = np.zeros ( [ n, n ] )
  for i in range ( 0, n ):
    if ( s[i] != 0.0 ):
      S[i,:] = p * A[i,:] / s[i] + ( 1.0 - p ) / n
    else:
      S[i,:] = 1.0 / n
  S = np.transpose ( S )
  return S

def rayleigh ( A, x ):

#*****************************************************************************80
#
## rayleigh() computes the Rayleigh estimate for an eigenvalue.
#
  import numpy as np
  value = np.dot ( x, np.dot ( A, x ) ) / np.dot ( x, x )
  return value

def transition_plot ( T, show_edge_labels, show_node_labels, 
  edge_labels = None, node_labels = None, filename = 'transition' ):

#*****************************************************************************80
#
## transition_plot() plots a digraph from its transition matrix.
#
  from graphviz import Digraph

  print ( node_labels )
  print ( show_node_labels )

  dot = Digraph ( format = 'png' )
  n = T.shape[0]
  for i in range ( 0, n ):
    code = str ( i )
    if ( show_node_labels ):
      dot.node ( str ( i ), str(i) + ': ' + str ( node_labels[i] ) )
    else:
      dot.node ( str ( i ), str ( i ) )

  for i in range ( 0, n ):
    for j in range ( 0, n ):
      if ( T[i,j] ):
        if ( show_edge_labels ):
          dot.edge ( str ( j ), str ( i ), str ( edge_labels[i,j] ) )
        else:
          dot.edge ( str ( j ), str ( i ) )

  dot.render ( filename + '.dot', view = True )
  print ( '  Graphics saved as "' + filename + '.dot.png' + '"' )

  return

def transition_plot2 ( T, show_edge_labels, show_node_labels, 
  edge_labels = None, node_labels = None, filename = 'transition' ):

#*****************************************************************************80
#
## transition_plot2() plots a graph from its transition matrix.
#
  from graphviz import Graph

  print ( node_labels )
  print ( show_node_labels )

  dot = Graph ( format = 'png' )
  n = T.shape[0]
  for i in range ( 0, n ):
    code = str ( i )
    if ( show_node_labels ):
      dot.node ( str ( i ), str(i) + ': ' + str ( node_labels[i] ) )
    else:
      dot.node ( str ( i ), str ( i ) )

  for i in range ( 0, n ):
    for j in range ( i, n ):
      if ( T[i,j] ):
        if ( show_edge_labels ):
          dot.edge ( str ( j ), str ( i ), str ( edge_labels[i,j] ) )
        else:
          dot.edge ( str ( j ), str ( i ) )

  dot.render ( filename + '.dot', view = True )
  print ( '  Graphics saved as "' + filename + '.dot.png' + '"' )

  return

pagerank_test ( )