// quantum computing emulation: Quantum Walk Page Rank
//
// algorithm from "Google in a Quantum Network", Paparo and Martin-Delgado,
// Scientific Reports, 2012, https://doi.org/10.1038/srep00444
//
// Usage: QW-PR [-adE] [T] < in
//
//   T = number of time steps, default = 10
//  -a  compute and display cumulative average probabilties
//  -d  extra debug output, -dd for even more output
//  -E  show entanglement
//
// input data format: number_of_nodes node_1_to_list... 0 ... node_N_to_list... 0
//
// R. Perry, August 2019
//
#include <stdio.h>
#include <stdlib.h> // atoi()
#include <math.h>

#define REAL
#include "qce.h"

#define alpha 0.85 // damping parameter

unsigned int aflag = 0, navg = 0; double *avg; // for cumulative average probabilities

// print state vector in matrix form
//
void matprint( const char *msg, State g, unsigned int n, unsigned long N)
{
  unsigned long i, k, p;  printf( "%s:\n", msg);

  for( k = 0; k < N; ++k) // for each lower state index (to)
  {
    for( i = 0; i < N; ++i) // for each upper state index (from)
    {
      p = (i << n) | k; // state index = upper index | lower index
      printf( " %7.4g", g.a[p]);
    }
    putchar( '\n');
  }
}

// create g coefficients from input adjacency values
//
// input data refers to nodes 1...N but internally we use 0...(N-1)
//
void init_g( State g, unsigned int n, unsigned long N)
{
  unsigned long i, j, k, p, count;  double v, N1 = 1.0/N;  zero( g);

  for( i = 0; i < N; ++i) // for each node (from)
  {
    j = i << n; // upper state index
    count = 0;  // number of outgoing branches
    while(1)
    {
      if( scanf( "%lu", &k) != 1) error( "QW-PR: bad input");
      if( k == 0) break;
      ++count; if( count > N) error( "QW-PR: too many outgoing states");
      if( k > N) error( "QW-PR: value out of range");
      p = j | (k-1); // state index = upper index | lower index
      if( g.a[p] != 0) error( "QW-PR: duplicate to-node value");
      g.a[p] = 1;
    }
    for( k = 0; k < N; ++k) 
    {
      v = (count == 0) ? N1 : g.a[j|k]/count;    // patched connectivity
      g.a[j|k] = sqrt( alpha*v + (1-alpha)*N1 ); // sqrt(G)
    }
  }
}

// initialize q = g/sqrt(N)
//
void init_q( State q, State g, unsigned long N)
{
  double N1 = sqrt(1.0/N);  for( unsigned long i = 0; i < g.N; ++i) q.a[i] = N1*g.a[i];
}

// first half of time step
//
void step1( State q, State g, unsigned int n, unsigned long N)
{
  unsigned long i, j, k, p;  double v;

  for( i = 0; i < N; ++i) // for each node (from)
  {
    j = i << n;  v = 0;

    for( k = 0; k < N; ++k) { p = j|k; v += g.a[p]*q.a[p]; }

    for( k = 0; k < N; ++k) { p = j|k; q.a[p] = 2*g.a[p]*v-q.a[p]; }
  }
}

// second half of time step, using reversed upper/lower state indices on q
//
void step2( State q, State g, unsigned int n, unsigned long N)
{
  unsigned long i, j, k, p, r;  double v;

  for( i = 0; i < N; ++i) // for each node (from)
  {
    j = i << n;  v = 0;

    for( k = 0; k < N; ++k) { p = j|k; r = (k<<n)|i; v += g.a[p]*q.a[r]; }

    for( k = 0; k < N; ++k) { p = j|k; r = (k<<n)|i; q.a[r] = 2*g.a[p]*v-q.a[r]; }
  }
}

// display ranks at time t
//
void disp( State q, unsigned int t, unsigned int n, unsigned long N)
{
  unsigned long i, k, p;  double v;  printf( "%i", t);

  ++navg; // count for averages

  for( k = 0; k < N; ++k) // for each lower state index (to)
  {
    v = 0; for( i = 0; i < N; ++i) { p = (i<<n)|k; v += abs2(q.a[p]); }

    if( aflag) // update average
    {
      if( navg != 1) v = avg[k] + (v-avg[k])/navg;
      avg[k] = v;
    }

    printf( " %.8g", v);
  }

  putchar( '\n');
}

//---------------------------------------------------------------------

int main( int argc, char *argv[])
{
  int debug = 0, Eflag = 0;

 // options
 //
  if( argc > 1 && argv[1][0] == '-') 
  {
    for( int i = 1; argv[1][i]; ++i) switch( argv[1][i])
    {
      case 'a': ++aflag; break;
      case 'd': ++debug; break;
      case 'E': ++Eflag; break;
      default: error( "QW-PR: bad option");
    }
    --argc, ++argv;
  }

 // number of time steps
 //
  unsigned int T = 10;  if( argc > 1) T = atoi(argv[1]);

 // number of nodes: N <= 2**n
 //
  unsigned long N;  if( scanf( "%lu", &N) != 1 || N < 2) error( "QW-PR: bad N input");

  unsigned int n = 0;  while( (1LU << n) < N) ++n;

  State g = Q(2*n); // google coefficients

  printf( "QW-PR: N = %lu, n = %u, 2**n = %lu, g.n = %u, g.N = %lu\n",
    N, n, 1LU << n, g.n, g.N);

  init_g( g, n, N);  if( debug) matprint( "sqrt(G)", g, n, N);

  State q = Q(g.n);  init_q( q, g, N);

 // setup for average probabilities
 //
  if( aflag) avg = emalloc(N*sizeof(double));

  if( debug > 1) { print( "q", q, PRINT_BITS); if( Eflag) print( "E", q, PRINT_TANGLE); }

  if( debug) matprint( "q", q, n, N);

  disp( q, 0, n, N);

 // T time steps
 //
  for( unsigned int t = 1; t <= T; ++t)
  {
    step1( q, g, n, N);  if( debug) matprint( "step1", q, n, N);

    if( debug > 1) { print( "q", q, PRINT_BITS); if( Eflag) print( "E", q, PRINT_TANGLE); }

    step2( q, g, n, N);  if( debug) matprint( "step2", q, n, N);

    if( debug > 1) { print( "q", q, PRINT_BITS); if( Eflag) print( "E", q, PRINT_TANGLE); }

    disp( q, t, n, N);
  }

  return 0;
}