/*************************************************************/
/* Search Algorithms for Unstructured Networks               */
/*         Pawel Pralat (01-Dec-2005)                        */
/*************************************************************/
/* 10-Jun-2005: first version (G(n,p) generating, flooding)  */
/* 08-Jul-2005: generating connected G(n,p)                  */
/* 26-Jul-2005: generating connected d-regular graph         */
/* 06-Sep-2005: generating clustered network                 */
/* 16-Sep-2005: generating power law graphs (protean graphs) */
/* 26-Sep-2005: flooding-walking algorithm                   */
/* 07-Nov-2005: normalized flooding                          */
/* 01-Dec-2005: k-walkers                                    */
/*************************************************************/

/*************************************************************/
/* type of searching                                         */
/*************************************************************/

#define s_flooding                         /* flooding */
// #define s_rand_walk                        /* random walk */
// #define s_fwf                              /* flooding, walking */
// #define s_k_walk                           /* k-walkers starting from the same vertex */
// #define s_norm_flood                       /* normalized flooding -- you have to define nf_n */
// int nf_n = 5;                              // number of neighbours

/*************************************************************/
/* type of graph                                             */
/*************************************************************/

#define g_gnp                             /* G(n,p) */
// #define g_d_reg                           /* d-regular */
// #define g_p2p                             /* P2P network */
// #define g_prot                            /* power law graph (protean) -- you have to define eta */
// double eta=0.91;                                // eta = 0.91 then gamma = 2.1
// double eta=0.59;                                // eta = 0.59 then gamma = 2.7

#include <stdio.h>
#include <math.h>

struct node
{
  long v;
  struct node * next;
};

long n;                     /* number of vertices */
int ex;                     /* expected degree */
long length;                /* number of querying */
struct node ** e;           /* edges (lists of neighbours) */
long ** edge;               /* edges (table of neighbours) */
long * deg;                 /* degrees */    

/* defs for rand2 */
#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)
/*end ran2 defs */

long myseed = -131;

/* random number generator from Press et al                                  */
/* Call with idum a negative integer to initialize; thereafter, do not alter */
/* idum between successive deviates in a sequence.                           */
/* RNMX should approximate the largest floating value that is less than 1.   */

float ran2(long *idum)
{
  int j;
  long k;
  static long idum2=123456789;
  static long iy=0;
  static long iv[NTAB];
  float temp;
    
  if(*idum <= 0)
    {                            /* initialise  */
      if (-(*idum)<1) *idum=1; /* preventing idum=0 */
      else *idum = -(*idum);
      idum2 = (*idum);
      for (j=NTAB+7;j>=0;j--)
        {                        /*load the shuffle table (after 8 warm-ups) */
	  k = (*idum)/IQ1;
	  *idum = IA1*(*idum-k*IQ1)-k*IR1;
	  if (*idum < 0) *idum += IM1;
	  if (j< NTAB) iv[j] = *idum;
	}
      iy = iv[0];
    }
  
  k=(*idum)/IQ1;            /* start here when not initialising */
  *idum = IA1*(*idum - k*IQ1)-k*IR1; /* compute idum = (IA1*idum) %IM1 without   */      
  if (*idum < 0) *idum += IM1;       /*   overflows by Schrage's method          */

  k = idum2/IQ2;
  idum2 = IA2*(idum2-k*IQ2) - k*IR2;  /* compute idum2 = (IA2*idum) %IM2 lokewise */
  if (idum2 < 0) idum2 += IM2;
  j = iy/NDIV;                       /* will be in range 0 .. NTAB - 1 */
  iy = iv[j] - idum2;                /* Here idum is shuffled, idum and idum2 are */
  iv[j] = *idum;                     /* combined to generate output */
  if(iy < 1) iy += IMM1;
  if ((temp=AM*iy) > RNMX) return RNMX; /* users don't expect endpoint value */
  else return temp;
}

long RAN(long *n)
{
  return floor(ran2(&myseed)*(*n));
}

int main(int argc, char * argv[])
{

#if defined s_fwf

  if (argc < 7)
    {
      printf("Usage: program <number of vertices> <E(deg)> <number of querying> <number of iterations> <output file> <number of runners> [seed] \n\n");
      return 0;
    }

  long fwf= atol(argv[6]);

  if (argc == 8)
    myseed = atol(argv[7]);

#elif defined s_k_walk

  if (argc < 7)
    {
      printf("Usage: program <number of vertices> <E(deg)> <number of querying> <number of iterations> <output file> <number of runners> [seed] \n\n");
      return 0;
    }

  long k_walk= atol(argv[6]);

  if (argc == 8)
    myseed = atol(argv[7]);

#else

  if (argc < 6)
    {
      printf("Usage: program <number of vertices> <E(deg)> <number of querying> <number of iterations> <output file> [seed] \n\n");
      return 0;
    }

  if (argc == 7)
    myseed = atol(argv[6]);

#endif

  long iters,it,i,j,first,last,vertex;

  /* program parameters */
  n = atol(argv[1]);
  ex = atol(argv[2]);
  length = atol(argv[3]);
  iters = atol(argv[4]);
  FILE * out = fopen(argv[5],"w");
  if (out == NULL) printf("FILE OPEN ERROR!");

  if(myseed>0) myseed=-myseed;

  /* memory allocation */
  e = (struct node **) malloc(n*sizeof(struct node *));
  edge = (long **) malloc(n*sizeof(long *));
  deg = (long *) malloc(n*sizeof(long));
  int * visited = (int *) malloc(n*sizeof(int));
  long * queue = (long *) malloc(n*sizeof(long));
  long * prev = (long *) malloc(n*sizeof(long));
  long * vector = (long *) malloc((n+1)*sizeof(long));
  double * final = (double *) malloc((n+1)*sizeof(double));

  for (i=0 ; i<n ; i++)
    {
      e[i]=NULL;
      deg[i]=0;
    }
  
  long max;

/* G(n,p) (WITH MIN DEGREE > 0) GRAPH GENERATION - BEGIN */
#ifdef g_gnp

  max = ((unsigned long) -1)/3;
  double coin;

  it=0;
  for (i=0 ; i<n ; i++)
    {
      for (j=i+1 ; j<n ; j++)
	{
	  coin = 1.0*RAN(&max)/max;    /* random number from (0,1) */
	  if (n*coin < ex)
	    {
	      deg[i]++;
	      struct node * new_node = (struct node *) malloc(sizeof(struct node));
	      new_node->v = j;
	      new_node->next = e[i];
	      e[i] = new_node;
	    
	      deg[j]++;
	      new_node = (struct node *) malloc(sizeof(struct node));
	      new_node->v = i;
	      new_node->next = e[j];
	      e[j] = new_node;
	    }
	}

      if (deg[i] == 0)   /* no isolated vertices */
	{
	  it++;
	  do {
	    j = RAN(&n);
	  } while (j == i);

	  deg[i]++;
	  struct node * new_node = (struct node *) malloc(sizeof(struct node));
	  new_node->v = j;
	  new_node->next = e[i];
	  e[i] = new_node;

	  deg[j]++;
	  new_node = (struct node *) malloc(sizeof(struct node));
	  new_node->v = i;
	  new_node->next = e[j];
	  e[j] = new_node;
	}
    }
    fprintf(out, "%dI", it);

#endif
/* G(n,p) (WITH MIN DEGREE > 0) GRAPH GENERATION - END */

/* d-REGULAR GRAPH GENERATION (d-process) - BEGIN */
#ifdef g_d_reg

    long remain, count, q_i, q_j;
    long big=100;
    int found=1;
    struct node * tmp;
    
    it=0;
    do 
      {
	it++;
	remain = n;
	
	for (i=0 ; i<n ; i++)
	  queue[i]=i;
      
	found=1;
	while ((remain>0) && (found==1))
	{
	  count=0;
	  while ((count<big) && (found==1))
	    {
	      count++;
	      q_i=RAN(&remain);
	      q_j=q_i; 
	      while (q_i==q_j)
		q_j=RAN(&remain);
	      i=queue[q_i];
	      j=queue[q_j];
	      
	      tmp = e[i];
	      found=0;
	      while (tmp)
		{
		  if (tmp->v==j)
		    found=1;
		  tmp=tmp->next;
		}
	      
	      if (found==0)
		{
		  deg[i]++;
		  struct node * new_node = (struct node *) malloc(sizeof(struct node));
		  new_node->v = j;
		  new_node->next = e[i];
		  e[i] = new_node;
		  if (deg[i]==ex)
		    queue[q_i]=queue[--remain];
		  
		  deg[j]++;
		  new_node = (struct node *) malloc(sizeof(struct node));
		  new_node->v = i;
		  new_node->next = e[j];
		  e[j] = new_node;
		  if (deg[j]==ex)
		    queue[q_j]=queue[--remain];
		}
	    }

	  if ((count==big) || (remain==1))
	    {
	      found = 0;

	      /* we have to come back to empty graph */
	      for (i=0 ; i<n ; i++)
		{
		  while (e[i] != NULL)
		    {
		      struct node * tmp_node = e[i];
		      e[i] = e[i]->next;
		      free(tmp_node);
		    }
		}
	      
	      for (i=0 ; i<n ; i++)
		deg[i]=0;
	    }
	  else
	    {
	      found = 1;
	    }
	}
    } while (found==0);	       
	      
    fprintf(out, "%dG", it-1);

#endif
/* d-REGULAR GRAPH GENERATION (d-process) - END */

/* CLUSTERED GRAPH GENERATION - BEGIN */
#ifdef g_p2p

    int size=100;                        /* size of dense graphs */

    // connecting supernodes (generating 3-regular graph)

    long remain, count, q_i, q_j;
    long big=100;
    int found=1;
    struct node * tmp;
    
    it=0;
    do 
      {
	it++;
	remain = n/size;
	
	for (i=0 ; i<n ; i++)
	  queue[i]=i*size;
      
	found=1;
	while ((remain>0) && (found==1))
	{
	  count=0;
	  while ((count<big) && (found==1))
	    {
	      count++;
	      q_i=RAN(&remain);
	      q_j=q_i; 
	      while (q_i==q_j)
		q_j=RAN(&remain);
	      i=queue[q_i];
	      j=queue[q_j];
	      
	      tmp = e[i];
	      found=0;
	      while (tmp)
		{
		  if (tmp->v==j)
		    found=1;
		  tmp=tmp->next;
		}
	      
	      if (found==0)
		{
		  deg[i]++;
		  struct node * new_node = (struct node *) malloc(sizeof(struct node));
		  new_node->v = j;
		  new_node->next = e[i];
		  e[i] = new_node;
		  if (deg[i]==3)
		    queue[q_i]=queue[--remain];
		  
		  deg[j]++;
		  new_node = (struct node *) malloc(sizeof(struct node));
		  new_node->v = i;
		  new_node->next = e[j];
		  e[j] = new_node;
		  if (deg[j]==3)
		    queue[q_j]=queue[--remain];
		}
	    }

	  if ((count==big) || (remain==1))
	    {
	      found = 0;

	      /* we have to come back to empty graph */
	      for (i=0 ; i<n ; i+=size)
		{
		  while (e[i] != NULL)
		    {
		      struct node * tmp_node = e[i];
		      e[i] = e[i]->next;
		      free(tmp_node);
		    }
		}
	      
	      for (i=0 ; i<n ; i+=size)
		deg[i]=0;
	    }
	  else
	    {
	      found = 1;
	    }
	}
    } while (found==0);	       
	      
    fprintf(out, "%dG", it-1);

  // generating dense graphs

  max = ((unsigned long) -1)/3;
  double coin;

  it=0;
  int k;
  for (k=0 ; k<n/size ; k++)
    {
      for (i=k*size ; i<(k+1)*size ; i++)
	{
	  for (j=i+1 ; j<(k+1)*size ; j++)
	    {
	      coin = 1.0*RAN(&max)/max;    /* random number from (0,1) */
	      if (size*coin < ex)
		{
		  deg[i]++;
		  struct node * new_node = (struct node *) malloc(sizeof(struct node));
		  new_node->v = j;
		  new_node->next = e[i];
		  e[i] = new_node;
	    
		  deg[j]++;
		  new_node = (struct node *) malloc(sizeof(struct node));
		  new_node->v = i;
		  new_node->next = e[j];
		  e[j] = new_node;
		}
	    }

	  if (deg[i] == 0)   /* no isolated vertices */
	    {
	      it++;
	      do {
		j = RAN(&n);
	      } while (j == i);

	      deg[i]++;
	      struct node * new_node = (struct node *) malloc(sizeof(struct node));
	      new_node->v = j;
	      new_node->next = e[i];
	      e[i] = new_node;
	      
	      deg[j]++;
	      new_node = (struct node *) malloc(sizeof(struct node));
	      new_node->v = i;
	      new_node->next = e[j];
	      e[j] = new_node;
	    }
	}
    }
    fprintf(out, "%dI", it);

#endif
/* CLUSTERED GRAPH GENERATION - END */

/* POWER LAW GRAPH GENERATION - BEGIN */
#ifdef g_prot

    max = ((unsigned long) -1)/3;
    double coin;

    double sum=0;
    long nvisited=n;
    for (i=0 ; i<n ; i++)
	sum += exp(-eta*log(i+1));
    final[0]=0;
    visited[0]=0;
    queue[0]=0;
    for (i=1 ; i<n ; i++)
      {
	final[i]=final[i-1]+exp(-eta*log(i))/sum;
	visited[i]=0;
	queue[i]=i;
      }
    final[n]=final[n-1]+exp(-eta*log(n))+1;

    long tt=0;
    long curr, curr_ind;
    while (nvisited)
      {
	tt++;
	
	// choose random vertex
	curr_ind = RAN(&n);
	curr = queue[curr_ind];
	if (visited[curr]==0)
	  {
	    nvisited--;
	    visited[curr]=1;
	  }

	// move this vertex to the end of permutation
	for (i=curr_ind ; i<n-1 ; i++) queue[i]=queue[i+1];
	queue[n-1]=curr;

	// delete all edges
	while (e[curr] != NULL)
	  {
	    deg[curr]--;
	    struct node * tmp_node = e[curr];
	    e[curr] = e[curr]->next;
	    i = tmp_node->v;
	    free(tmp_node);

            deg[i]--;
	    if (e[i]->v==curr)
	      {
		tmp_node = e[i];
		e[i] = e[i]->next;
		free(tmp_node);
	      }
	    else
	      {
		tmp_node = e[i];
		while (tmp_node->next->v!=curr)
		  tmp_node = tmp_node->next;
		struct node * del_node = tmp_node->next;
		tmp_node->next = tmp_node->next->next;
		free(del_node);
	      }
	  }

	// generate ex edges
	for (i=0 ; i<ex ; i++)
	  {
	    int found;
	    do
	      {
		coin = 1.0*RAN(&max)/max;    /* random number from (0,1) */
		long i_down = 0;
		long i_up = n-1;
		long i_av = (i_down+i_up)/2;
		while ((final[i_av] > coin) || (final[i_av+1] <= coin))
		  {
		    if (final[i_av] > coin)
		      {
			i_up = i_av-1;
			i_av = (i_down+i_up)/2;
		      }
		    else
		      {
			i_down = i_av+1;
			i_av = (i_down+i_up)/2;
		      }
		  }
		
		j=queue[i_av]; 

		// there are no parallel edges
		found=0;
		struct node * tmp = e[curr];
		while ((tmp) && (found==0))
		  {
		    if (tmp->v==j)
		      found=1;
		    tmp=tmp->next;
		  }
	      } while (found);

	    deg[curr]++;
	    struct node * new_node = (struct node *) malloc(sizeof(struct node));
	    new_node->v = j;
	    new_node->next = e[curr];
	    e[curr] = new_node;

	    deg[j]++;
	    new_node = (struct node *) malloc(sizeof(struct node));
	    new_node->v = curr;
	    new_node->next = e[j];
	    e[j] = new_node;
	  }
      }

    it=0;
    for (i=0 ; i<n ; i++)
      if (deg[i] == 0)   /* no isolated vertices */
	{
	  it++;
	  do {
	    j = RAN(&n);
	  } while (j == i);
	  
	  deg[i]++;
	  struct node * new_node = (struct node *) malloc(sizeof(struct node));
	  new_node->v = j;
	  new_node->next = e[i];
	  e[i] = new_node;
	  
	  deg[j]++;
	  new_node = (struct node *) malloc(sizeof(struct node));
	  new_node->v = i;
	  new_node->next = e[j];
	  e[j] = new_node;
	}
    fprintf(out, "%dI", it);

    for (i=0 ; i<10 ; i++)
      printf("nubmer of neighbours (%d) = %d\n", i, deg[queue[i]]);

#endif
/* POWER LAW GRAPH GENERATION - END */

  it=0;
  do {                                                  // we need connected graph 
    visited[0]=1;
    for (i=1 ; i<n ; i++)
      visited[i]=0;

    vertex = 0;
    first = last = 0;
    queue[last++] = vertex;

    while (first < last)
      {
	vertex = queue[first++];

	struct node * tmp = e[vertex];
	while (tmp)
	  {
	    if (visited[tmp->v]==0)
	      {
		visited[tmp->v]=1;
		queue[last++]=tmp->v;
	      }
	    tmp=tmp->next;
	  }
      }

    if (first!=n)                                                                  // graph is disconnected!
      {
	it++;
	
	j=queue[first-1];
	i=0;
	while(visited[i])
	  i++;
                                              	// add edge (i,j) between components
	deg[i]++;                      
	struct node * new_node = (struct node *) malloc(sizeof(struct node));
	new_node->v = j;
	new_node->next = e[i];
	e[i] = new_node;

	deg[j]++;
	new_node = (struct node *) malloc(sizeof(struct node));
	new_node->v = i;
	new_node->next = e[j];
	e[j] = new_node;
      }
  } while (first != n);

  fprintf(out, "%dD ", it);

  // we have a graph! let's build the table of neighbours

  for (i=0 ; i<n ; i++)
    {
      edge[i] = (long *) malloc(deg[i]*sizeof(long));
      struct node * tmp = e[i];
      j=0;
      while (tmp)
	{
	  edge[i][j]=tmp->v;
	  tmp=tmp->next;
	  j++;
	}      
    }

  // we are ready to test the sampling algorithms

  for (i=0 ; i<=n ; i++)
    final[i]=0;

  double nr_queries = 0;   // the average number of queries for flooding (FW-algorithm)

  max = 0;
  for (it=0 ; it<iters ; it++)
    {
      for (i=0 ; i<=n ; i++)
	vector[i]=0;

      for (i=0 ; i<n ; i++)
	visited[i]=0;

      vertex = it;

/* FLOODING - BEGIN */
#ifdef s_flooding
 
      first = 0;
      last = 0;
      prev[last] = -1;
      queue[last++] = vertex;
      visited[vertex]++;

      while (last < length)
	{
	  long pr = prev[first];
	  vertex = queue[first++];

	  struct node * tmp = e[vertex];
	  while ((tmp) && (last < length))
	    {
	      if (tmp->v != pr)
		{
		  visited[tmp->v]++;
		  prev[last] = vertex;
		  queue[last++] = tmp->v;
		}
	      tmp=tmp->next;
	    }
	}
#endif
/* FLOODING - END */

/* NORMALIZED FLOODING - BEGIN */
#ifdef s_norm_flood

      first = 0;
      last = 0;
      prev[last] = -1;
      queue[last++] = vertex;
      visited[vertex]++;

      while (last < length)
        {
          long pr = prev[first];
          vertex = queue[first++];

	  if (deg[vertex]-1 <= nf_n)
	    {                                     // we have to use all neighbours
	      struct node * tmp = e[vertex];

	      while ((tmp) && (last < length))
		{
		  if (tmp->v != pr)
		    {
		      visited[tmp->v]++;
		      prev[last] = vertex;
		      queue[last++] = tmp->v;
		    }
		  tmp=tmp->next;
		}
	    }
	  else
	    {                             // we have to use nf_s random neighbours
	      long tmp;
	      for (i = deg[vertex] ; i > deg[vertex]-nf_n ; i--)
		if (last < length)
		  {
		    do {
		      j = RAN(&i);
		    } while (edge[vertex][j] == pr);

		    visited[edge[vertex][j]]++;
		    prev[last] = vertex;
		    queue[last++] = edge[vertex][j];
		    
		    tmp = edge[vertex][j];
		    edge[vertex][j] = edge[vertex][i-1];
		    edge[vertex][i-1] = tmp;
		  }
	    }
        }
#endif
/* NORMALIZED FLOODING - END */
 
/* RANDOM WALK - BEGIN */
#ifdef s_rand_walk

      long pr = -1;
      visited[vertex]++;

      for (i=1 ; i<length ; i++)
	{
	  if (deg[vertex] == 1)
	    {
	      visited[edge[vertex][0]]++;
	      pr = vertex;
	      vertex = edge[vertex][0];
	    }
	  else
	    {
	      long nr;
	      do {
		nr = RAN(&deg[vertex]);
	      } while (edge[vertex][nr] == pr);

	      visited[edge[vertex][nr]]++;
	      pr = vertex;
              vertex = edge[vertex][nr];   
            }
        }

#endif
/* RANDOM WALK - END */

/* K RANDOM WALKERS - BEGIN */
#ifdef s_k_walk

      for (j=0 ; j<k_walk ; j++)
	{
	  vertex = it;
	  long pr = -1;
	  visited[vertex]++;

	  for (i=1 ; i<(length/k_walk) ; i++)
	    {
	      if (deg[vertex] == 1)
		{
		  visited[edge[vertex][0]]++;
		  pr = vertex;
		  vertex = edge[vertex][0];
		}
	      else
		{
		  long nr;
		  do {
		    nr = RAN(&deg[vertex]);
		  } while (edge[vertex][nr] == pr);
		  
		  visited[edge[vertex][nr]]++;
		  pr = vertex;
		  vertex = edge[vertex][nr];
		}
	    }
	}

#endif
/* K RANDOM WALKERS - END */

/* FLOODING, WALKING - BEGIN */
#ifdef s_fwf

      //FLOODING

      first = 0;
      last = 0;
      prev[last] = -1;
      queue[last++] = vertex;
      visited[vertex]++;
      long pr;

      while (last-first < fwf)
        {
          pr = prev[first];
          vertex = queue[first++];

          struct node * tmp = e[vertex];
          while ((tmp) && (last-first < fwf))
            {
              if (tmp->v != pr)
                {
                  visited[tmp->v]++;
                  prev[last] = vertex;
                  queue[last++] = tmp->v;
                }
              tmp=tmp->next;
            }
        }

      nr_queries += 1.0*(last-1)/iters;               // number of quering for flooding

      // WALKING

      long vertex_ind=first;

      for (i=last ; i<length ; i++)
        {
	  pr = prev[vertex_ind];
	  vertex = queue[vertex_ind];

          if (deg[vertex] == 1)
            {
              prev[vertex_ind] = vertex;
              queue[vertex_ind] = edge[vertex][0];
	      visited[edge[vertex][0]]++;
            }
          else
            {
              long nr;
              do {
                nr = RAN(&deg[vertex]);
              } while (edge[vertex][nr] == pr);

              visited[edge[vertex][nr]]++;
              prev[vertex_ind] = vertex;
              queue[vertex_ind] = edge[vertex][nr];
            }

	  vertex_ind++;
          if (vertex_ind == last)
            vertex_ind = first;
        }

#endif
/* FLOODING, WALKING - END */

      for (i=0 ; i<n ; i++)
	vector[visited[i]]++;

      for (i=0 ; i<=n ; i++)
	if (vector[i])
	  {
	    final[i] += 1.0*vector[i]/iters;
	    if (i > max) max = i;
	  }
    }

#ifdef s_fwf

  fprintf(out, "R%dF%dW%d ", fwf, ((int)nr_queries), length-((int)nr_queries));
  
#endif

#ifdef s_k_walk

  fprintf(out, "W%d ", k_walk);

#endif

  for (i=0 ; i<=max ; i++)
    fprintf(out, "%f ", final[i]);
  fprintf(out, "\n");
  
  /* memory deallocation */
  for (i=0 ; i<n ; i++)
    {
      free(edge[i]);
      while (e[i] != NULL)
	{
	  struct node * tmp_node = e[i];
	  e[i] = e[i]->next;
	  free(tmp_node);
	}
    }
  free(e);
  free(edge);
  free(deg);
  free(visited);
  free(queue);
  free(vector);
  free(final);
  free(prev);

  fclose(out);

  return 0;
}