#! /usr/bin/env python3
#
def power_rank ( A ):

#*****************************************************************************80
#
## power_rank() uses the power method for an eigenvector of a directed graph.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    10 April 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer A(N,N): the incidence matrix.
#
#  Output:
#
#    real X(N): the estimated eigenvector.
#
  from transition_from_incidence import transition_from_incidence
  import numpy as np

  steps = 200

  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 = transition_from_incidence ( A )
#
#  Set the starting vector.
#
  x = np.ones ( n ) / float ( n )
#
#  Carry out many iterations.
#  Normalization is not necessary because T is a transition matrix.
#
  for i in range ( 0, steps ):
    x = np.dot ( T, x )
#
#  Compare three successive iterates.
#
  Tx  = np.dot ( T, x )
  TTx = np.dot ( T, Tx )

  W = np.stack ( ( x, Tx, TTx ), axis = 1 )
  print ( '' )
  print ( 'x, T*x, T*T*x' )
  print ( W )

  return x

def power_rank_test ( ):

#*****************************************************************************80
#
## power_rank_test() tests power_rank().
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    10 April 2022
#
#  Author:
#
#    John Burkardt
#
  from moler_incidence import moler_incidence

  print ( '' )
  print ( 'power_rank_test():' )
  print ( '  power_rank() computes the power method rankings r()' )
  print ( '  for an incidence matrix A.' )
  print ( '' )
  print ( '  The Moler digraph has a sink node.' )
  print ( '  This means there will only be a single nonzero ranking!' )
 
  A = moler_incidence ( )
  print ( '' )
  print ( '  Incidence matrix A:' )
  print ( A )

  r = power_rank ( A )
  print ( '' )
  print ( '  power rankings r:' )
  print ( r )

  return

if ( __name__ == '__main__' ):
  power_rank_test ( )