//
//  Peter Richter (prichter@princeton.edu)
//  Princeton University Undergraduate Mathematics Laboratory
//  May 11, 2001
//
//  walk.cpp - compute the second-largest eigenvalue of each bipartite graph
//             of order N and regularity K along a random walk of length M.
//

#include <iostream.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define FALSE    0
#define TRUE     1

#define K        3       /*   number of edges meeting each vertex      */
#define N      500       /*   even number of vertices in the graph     */
#define M     1000       /*   number of steps in the random walk       */
#define EPS   1e-8       /*   tolerance for power method convergence   */

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

double error_bound(double e, double x[N], int g[N][K])
{
  /* bound the error in the eigenvalue e using the eigenvector x of e^2 */

  // compute the 2-norm of x
  double x_norm = 0;
  for (int i = 0; i < N; i++)
    x_norm += x[i] * x[i];
  x_norm = sqrt(x_norm);

  // normalize x
  double y[N];
  for (int i = 0; i < N; i++)
    y[i] = x[i] / x_norm;

  // apply the adjacency matrix to y
  double t[N];
  for (int i = 0; i < N; i++) {
    t[i] = 0;
    for (int j = 0; j < K; j++)
      t[i] += y[g[i][j]];
  }

  // apply the adjacency matrix to t
  double z[N];
  for (int i = 0; i < N; i++) {
    z[i] = 0;
    for (int j = 0; j < K; j++)
      z[i] += t[g[i][j]];
  }

  // compute the error vector zeta
  double zeta[N];
  for (int i = 0; i < N; i++)
    zeta[i] = z[i] - (e * e) * y[i];

  // compute the 2-norm of zeta
  double zeta_norm = 0;
  for (int i = 0; i < N; i++)
    zeta_norm += zeta[i] * zeta[i];
  zeta_norm = sqrt(zeta_norm);

  // return the error bound on e
  return sqrt(zeta_norm);
}

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

double eig_two(int g[N][K], double p)
{
  /* compute the second largest eigenvalue of the adjacency matrix for g */

  // choose an initial vector x if this routine has not been called before
  static int first_call = TRUE;
  static double x[N];
  if (first_call == TRUE) {
    for (int i = 0; i < N; i++)
      x[i] = ((double) rand()) / RAND_MAX;
    first_call = FALSE;
  }

  // calculate the second eigenvalue using the power method
  double new_approx = -K-1, old_approx = -K-2;
  while (fabs(new_approx - old_approx) >= EPS ||
         error_bound(new_approx, x, g) >= fabs(new_approx - p)) {

    // compute the orthogonality coefficients
    double c = 0, d = 0;
    for (int i = 0; i < N; i++) {
      c += x[i];
      d += (i < N/2) ? x[i] : -x[i];
    }
    c /= N;
    d /= N;

    // project out the two principal eigenvectors from x
    double y[N];
    for (int i = 0; i < N; i++)
      y[i] = (i < N/2) ? x[i]-c-d : x[i]-c+d;

    // compute the 2-norm of y
    double y_norm = 0;
    for (int i = 0; i < N; i++)
      y_norm += y[i] * y[i];
    y_norm = sqrt(y_norm);

    // normalize y
    double z[N];
    for (int i = 0; i < N; i++)
      z[i] = y[i] / y_norm;

    // apply the adjacency matrix to z
    double t[N];
    for (int i = 0; i < N; i++) {
      t[i] = 0;
      for (int j = 0; j < K; j++)
        t[i] += z[g[i][j]];
    }

    // apply the adjacency matrix to t
    for (int i = 0; i < N; i++) {
      x[i] = 0;
      for (int j = 0; j < K; j++)
        x[i] += t[g[i][j]];
    }

    // compute the square root of (x,z)
    double lambda = 0;
    for (int i = 0; i < N; i++)
      lambda += x[i] * z[i];
    lambda = sqrt(lambda);

    // store this estimate
    old_approx = new_approx;
    new_approx = lambda;
  }
  return new_approx;
}

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

void shuffle(int a[])
{
  /* permute randomly the elements of the array a */

  // do the fisher-yates shuffle!
  for (int i = (N/2)-1; i >= 0; i--) {
    int j = rand() % (i+1);
    int t = a[i];
    a[i] = a[j];
    a[j] = t;
  }
}

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

void walk_init(int g[N][K])
{
  /* begin the random walk along the space of bipartite K-regular graphs */

  // seed the random number generator
  long seed = time(NULL);
  srand(seed);

  int okay;              // flag to indicate the graph is simple and connected
  int v[K][N/2];         // arrays to hold K permutations of vertices N/2...N-1

  do {

    // reset flag to true
    okay = TRUE;

    // fill each array v[i] with the independent set of vertices N/2...N-1
    for (int i = 0; i < K; i++)
      for (int j = 0; j < N/2; j++)
        v[i][j] = N/2 + j;

    // shuffle each array v[i] randomly; each (j, v[i][j]) represents an edge
    for (int i = 0; i < K; i++)
      shuffle(v[i]);

    // check that the graph is simple
    for (int i = 0; i < N/2; i++)
      for (int j = 0; j < K; j++)
        for (int r = 0; r < K; r++)
          if (r != j && v[j][i] == v[r][i])
            okay = FALSE;

    // re-shuffle if the graph is not simple
    if (okay == FALSE)
      continue;

    // store the graph in the N-by-K table g
    for (int i = 0; i < K; i++)
      for (int j = 0; j < N/2; j++)
        g[j][i] = v[i][j];
    for (int i = N/2; i < N; i++)
      for (int j = 0, s = 0; j < N/2; j++)
        for (int r = 0; r < K; r++)
          if (g[j][r] == i)
            g[i][s++] = j;

    // check that the graph is connected; loop again if it is not
    if (fabs(eig_two(g, 1/EPS) - K) < sqrt(EPS))
      okay = FALSE;

  } while (okay == FALSE);
}

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

void walk_step(int g[N][K])
{
  /* step to the next graph in the random walk */

  int okay;              // flag to indicate the next step preserves structure

  do {

    // reset flag to true
    okay = TRUE;

    // pick two distinct vertices from independent set of vertices 0...(N/2)-1
    int v2, v1 = rand() % (N/2);
    do {
      v2 = rand() % (N/2);
    } while (v2 == v1);

    // pick distinct vertices connected to v1 (but not v2) and v2 (but not v1)
    int i1 = rand() % K, w1 = g[v1][i1];
    int i2 = rand() % K, w2 = g[v2][i2];
    for (int i = 0; i < K; i++)
      if (w1 == g[v2][i] || w2 == g[v1][i])
        okay = FALSE;

    // loop again if this was not accomplished
    if (okay == FALSE)
      continue;

    // make the next step in the walk
    g[v1][i1] = w2;
    g[v2][i2] = w1;
    for (int i = 0; i < K; i++) {
      if (g[w2][i] == v2)
        g[w2][i] = v1;
      if (g[w1][i] == v1)
        g[w1][i] = v2;
    }

  } while (okay == FALSE);
}

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

void main()
{
  int g[N][K];            // the graph; each (i, g[i][j]) represents an edge
  double p = 2*sqrt(K-1); // critical value for characterizing ramanujan graphs
  double e;               // second largest eigenvalue of the adjacency matrix

  // perform the random walk, recording the eigenvalue of each matrix
  walk_init(g);
  for (int i = 0; i < M; i++) {
    if ((e = eig_two(g, p)) <= p)              // <- g is ramanujan
      cout << e << '\t' << '*' << endl;
    else                                       // <- g is not ramanujan
      cout << e << endl;
    walk_step(g);
  }
}

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