/* Generator for sparse geometric graphs as described in:

S. Poljak and F. Rendl: Solving the MAX-CUT Problem using Eigenvalues.
Report No. 91735-OR, Forschungsinstitute fur Diskrete Mathematik,
University of Bonn, Germany, 1991.

  "...we take the integer lattice {1,...,x_max} x {1,...,y_max} and
select each gridpoint with a certain probability p as a vertex of a graph.
We introduce an edge of weight one between two points, whenever their
Euclidean distance is at most a certain threshold delta_max. "

How to run the program: a.out x_max y_max p delta_max filename

Output format: briefly specified in the function write_graph()

written June 1994 by Marcus Peinado
*/

#include <stdio.h>

#define TRUE    1
#define FALSE   0
#define AND     &&
#define OR      ||
#define NOT     !
#define MOD     %
#define BOOL	char
#define randomNat(n)	(random() % (n))

#define MAXLONG         0x7fffffff
extern long random();

#define MAX_NR_VERTICES		20000
#define MAX_NR_VERTICESdiv8	1250

typedef struct {
  int v1, v2;
  double wt;
} F_EDGE;      /* the form in which the edges are stored in the file */


BOOL bitmap[ MAX_NR_VERTICES ][ MAX_NR_VERTICESdiv8 ];
int nr_vert, deg[ MAX_NR_VERTICES ], v_nr[ MAX_NR_VERTICES ];
double p, epsilon;
int xdim, ydim;

char masks[ 8 ] = { 0x01, 0x02, 0x04, 0x08, 0x10, 0x20, 0x40, 0x80 };

double randomFloat()
{
        return ( (double)random() / (double)MAXLONG );
}

void set_edge( register int i, register int j, BOOL x )
{
	register int byte, bit, mask;

	bit  = 7 - (j & 0x00000007);
	byte = j >> 3;

	mask = masks[bit];
	if ( x == 1 ) bitmap[i][byte] |= mask;
		else  bitmap[i][byte] &= ~mask;
}


BOOL get_edge( register int i, register int j )
/* if {i,j} is an edge then TRUE otherwise FALSE */
{
        register int byte, bit, mask;

        bit  = 7-(j & 0x00000007);
        byte = j >> 3;

        mask = masks[bit];
        return( (bitmap[i][byte] & mask)==mask );
}


int round( int p, int q )
{
	if ( (p MOD q) == 0 ) return( p / q );
	else return( 1 + p/q );
}


void write_graph( char *fn )
/* The output file is a sequence of edges in the form given by F_EDGE.
   Every edge appears twice: (u,w) and (w,u). Clearly, it would be suffcient
   to store each edge only once. */
{
        int i,j, nr_edges=0;
	F_EDGE e[2];
	FILE *fp;

	fp = fopen( fn,"w" );
	if ( fp == NULL ) { printf("can't open file\n"); exit(10); }
	e[0].wt = 1.0; e[1].wt = 1.0;
        for ( i=0; i<nr_vert; i++ )
	  for ( j=0; j<i; j++ )
	      if ( get_edge(i,j) )
		{
		  nr_edges++;
		  e[0].v1 = v_nr[i]; e[0].v2 = v_nr[j];
		  e[1].v2 = v_nr[i]; e[1].v1 = v_nr[j];
		  fwrite( (void *)e, sizeof( F_EDGE ),2,fp );
		}
	fclose( fp );
	printf("nr_edges = %d\n",nr_edges );
}


void eliminate_disconn()
/* purpose: eliminate all vertices which are isolated, i.e. have no edge
            incident on them. These vertices are irrelevant for any CUT.
            Only vertices with degree>=1 get a vertex number  */
            
{
  int i, v=0;
  
  for ( i=0; i<nr_vert; i++ )
    if ( deg[i] > 0 )
      {
	v_nr[ i ] = v;
	v++;
      }
    else v_nr[i] = -1;
  printf("nr_vert = %d   ",v );
}


void empty_graph()
{
	int i, j;

	for ( i=0; i<nr_vert; i++ )
	  {
	    deg[i] = 0;
	    for ( j=0; j<nr_vert; j++ )
	      set_edge(i,j,FALSE);
	  }
}


double distance( int i, int j )
{
  int xi, yi, xj, yj, dx, dy;

  xi = i % xdim; yi = i / xdim;
  xj = j % xdim; yj = j / xdim;
  dx = xj - xi;
  dy = yj - yi;

  return( dx * dx + dy * dy );
}


void geo_graph()
{
  int i, j, no;
  BOOL active[ nr_vert ];
  
  for ( i=0; i<nr_vert; i++ )
    if ( randomFloat() < p ) active[i] = TRUE;
    else active[i] = FALSE;

  for ( i=0; i<nr_vert; i++ )
    {
      set_edge(i,i,FALSE );
      if ( active[i] )
	for ( j=0; j<i; j++ )
	  if ( active[j] AND (distance(i,j)<=epsilon) )
	    {
	      set_edge( i,j,TRUE );
	      set_edge( j,i,TRUE );
	      deg[i]++; deg[j]++;
	    }
	  else
	    {
	      set_edge( i,j,FALSE );
	      set_edge( j,i,FALSE );
	    }
    }
}


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

	srandom( time( (char *)0 ) ); 

	if ( argc != 6 ) {printf("usage: %s xdim ydim p eps file\n",argv[0]);
			  exit(10); }

	sscanf( argv[1],"%d",&xdim );
	sscanf( argv[2],"%d",&ydim );
	sscanf( argv[3],"%lg",&p );
	sscanf( argv[4],"%lg",&epsilon );
	epsilon = epsilon * epsilon;
	nr_vert = xdim * ydim;

	geo_graph();
	eliminate_disconn();
	write_graph( argv[5] );
}