/*



              +==================================+
              | Random Graph Generation Software |
              +==================================+




       This source code is (c) Copyright 2001 by Christopher R. Palmer
       and J. Gregory Steffan.  It may be freely redistributed and
       included in other packages provided this copyright notice is
       included.




   Version 1.0
       January 12, 2001 - Initial public release

   ----------------------------------------------------------------------

   Using the generator:
   ====================

   This file is the entire distribution.  Compile it with your favorite
   C compiler and you can then run the program as follows:

      ./gen [-pgm fname | -pgm2 fname] logn m alpha beta gamma epsilon

    which will then randomly generate a graph with 2**logn nodes.
    There will be exactly m distinct edges in the graph (each edge
    is weighted).  The recursive probability parameters are

       alpha + beta + gamma + epsilon = 1

    For more information please see the paper that describes this
    graph generation algorithm.  A citation appears at the end of
    this document.

     The two optional arguments

          -pgm fname
          -pgm2 fname

     will create the named file in pgm format.  Using the first option
     will visualize the probabilities for the adjacency matrix (darker
     is less probable) and the second form will visualize the actual
     adjacency matrix created (darker is fewer links).  These options
     are mostly useful for debugging.


     Output format:
     ==============

     The output format is a simple ascii representation of the graph.
     The first line contains a single integer that is the number of
     nodes in the graph (2**logn).  Each subsequent line contains three
     integers:

        u v w

     which represents an edge from node u to node v with weight w.


     More Information / Etc:
     =======================

     The paper that describes this work is:

         Christopher R. Palmer and J. Gregory Steffan, Generating
         Network Topologies that Obey Power Laws, IEEE Globecom 2000,
         San Francisco, CA, November 2000.

         http://www.cs.cmu.edu/~crpalmer/pubs.html/globecom2000.ps

     please cite it when referring to this graph generation software.
     Further information, bug reports, etc. may be sent to:

         crpalmer@cs.cmu.edu


     ------------------------------------------------------------------

*/

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

#define SYMMETRIC 1

typedef struct hashEntS {
    int		w;
    int		u, v;
    struct hashEntS *next;
} hashEnt_t;

hashEnt_t **htab;
int        nhtab;

int n, nlinks;
double alpha, beta, Gamma, epsilon;

#define PGM_RAW		0
#define PGM_WIDTH	32		/* size of an atom in pgm file */
#define PGM_HEIGHT	32
#define PGM_MAX		255

#define mymax(a,b) ((a) > (b) ? (a) : (b))


inline int
hash(int uu, int vv)
{
    long long int u = uu, v = vv;

    return (u*n+v)%nhtab;
}

int
hashadd(int u, int v)
{
    hashEnt_t **h = &htab[hash(u,v)];

#if SYMMETRIC
    if (v < u) return hashadd(v, u);
#endif

    while (*h) {
	if ((*h)->u == u && (*h)->v == v) {
	    (*h)->w++;
	    return 0;
	}
	h = &(*h)->next;
    }

    if ((*h = malloc(sizeof(**h))) == NULL) {
	perror("malloc");
	exit(1);
    }

    (*h)->u = u;
    (*h)->v = v;
    (*h)->w = 1;
    (*h)->next = NULL;

    return 1;
}

void
gen(int n, int *u, int *v)
{
    double p;

    if (n <= 1) {
	*u = *v = 0;
	return;
    }

    p = random()/(RAND_MAX+1.0);
    gen(n/2, u, v);

    if (p < alpha) {
	/* nop because it's just + (0,0) */
    } else if (p < alpha+beta) {
	*v += n/2;
    } else if (p < alpha+beta+Gamma) {
	*u += n/2;
    } else {
	*u += n/2;
	*v += n/2;
    }
}

void
gensym(int n, int *u, int *v)
{
    double p;
    double T;
    double A;
    double F;
    double norm;

    if (n <= 1) {
	*u = *v = 0;
	return;
    }

    T = n*n-n;
    A = (n/2)*(n/2)-(n/2);
    F = (n/2)*(n/2);
    norm = 1/(T*((alpha+epsilon)/A + (Gamma+beta)/F));

    p = random()/(RAND_MAX+1.0);

    if (p < alpha*norm*(T/A)) {
	gensym(n/2, u, v);
	/* nop because it's just + (0,0) */
    } else if (p < (alpha+epsilon)*norm*(T/A)) {
	gensym(n/2, u, v);
	*u += n/2;
	*v += n/2;
    } else {
	gen(n/2, u, v);
	*v += n/2;
    }
}
    
/* genprob(A, sim, n, topr, topc, size, skip)
     A - matrix in which to store the result
     sim          is t/f does this need to be symmetric?
     n            is the number/prob of having edges in this matrix
     (topr, topc) is the top left corner of the submatrix
     size         is the width and height of the submatrix
     skip         is the number of elements in the entire matrix row
 */

void
genprob(double *prob, int sim, double n, int topr, int topc, int size, int skip)
{
    if (n < 0) return;
    if (n > 1) n = 1;

    if (sim && size == 2) {
	prob[topr*skip+topc+1] = n;
    } else if (size == 1) {
	prob[topr*skip+topc] = n;
    } else {
	if (sim) {
	    double T = size*size-size;
	    double A = (size/2)*(size/2)-(size/2);
	    double F = (size/2)*(size/2);
	    double norm = 1/(T*((alpha+epsilon)/A + (Gamma+beta)/F));

	    genprob(prob, 1, alpha*norm*(T/A)*n, topr, topc, size/2, skip);
	    genprob(prob, 1, epsilon*norm*(T/A)*n, topr+size/2, topc+size/2, size/2, skip);
	    genprob(prob, 0, (alpha+epsilon)*norm*(T/F)*n, topr, topc+size/2, size/2, skip);
	} else {
            genprob(prob, 0, alpha *n,  topr,        topc,        size/2, skip);
            genprob(prob, 0, beta*n,    topr,        topc+size/2, size/2, skip);
            genprob(prob, 0, Gamma*n,   topr+size/2, topc,        size/2, skip);
            genprob(prob, 0, epsilon*n, topr+size/2, topc+size/2, size/2, skip);
	}
    }
}

int
lincmp(const void *A, const void *B)
{
    const hashEnt_t *a = *(const hashEnt_t **) A;
    const hashEnt_t *b = *(const hashEnt_t **) B;

    if (a->u != b->u) return a->u - b->u;
    return a->v - b->v;
}

void
dumplinks(void)
{
    hashEnt_t *h;
    hashEnt_t **lin;
    int i, nlin;

    printf("%d\n", n);
    fflush(stdout);

    lin = malloc(sizeof(*lin)*nlinks);
    nlin = 0;

    for (i = 0; i < nhtab; i++) {
	for (h = htab[i]; h; h = h->next) {
	    lin[nlin++] = h;
	}
    }

    qsort(lin, nlin, sizeof(*lin), lincmp);

    for (i = 0; i < nlin; i++) {
	printf("%d %d %d\n", lin[i]->u, lin[i]->v, lin[i]->w);
    }
}

void
recordlinks(double *p, int n)
{
    hashEnt_t *h;
    int i;

    memset(p, 0, sizeof(*p)*n*n);

    for (i = 0; i < nhtab; i++) {
	for (h = htab[i]; h; h = h->next) {
	    p[h->u*n+h->v] = h->w;
	}
    }
}

void
dumppgm(char *fname, double *P, int n)
{
    FILE *f;
    int   r, c, i, j;
    double mx = 0;

    if ((f = fopen(fname, "w")) == NULL) {
	perror(fname);
	return;
    }

    for (r = 0; r < n; r++) for (c = 0; c < n; c++) mx = mymax(mx, P[r*n+c]);

    fprintf(stderr, "maximum %g\n", mx);

    fprintf(f, "P%c\n%d %d\n%d\n", PGM_RAW ? '5' : '2', n*PGM_WIDTH, n*PGM_HEIGHT, PGM_RAW ? PGM_MAX : (int) (mx*1000));

    for (r = 0; r < n; r++) {
	for (i = 0; i < PGM_HEIGHT; i++) {
	    for (c = 0; c < n; c++) {
		for (j = 0; j < PGM_WIDTH; j++) {
		    if (PGM_RAW) putc((unsigned char) (PGM_MAX - P[r*n+c]/mx*(PGM_MAX-1) + 1), f);
		    else fprintf(f, "%u ", (unsigned int) (mx*1000 - P[r*n+c]*1000));
		}
	    }
	    if (! PGM_RAW) fprintf(f, "\n");
	}
    }

    fclose(f);
}

int
main(int argc, char **argv)
{
    int genlinks;
    int u, v;
    double *prob;
    char *pgm_name = NULL;
    char *pgm2_name = NULL;

    for (;;) {
        if (argc >= 3 && strcmp(argv[1], "-pgm") == 0) {
	    pgm_name = argv[2];
	    argc -= 2;
	    argv += 2;
        } else if (argc >= 3 && strcmp(argv[1], "-pgm2") == 0) {
	    pgm2_name = argv[2];
	    argc -= 2;
	    argv += 2;
	} else break;
    }

    if (argc != 7 || (n = atoi(argv[1])) < 1 || (nlinks = atoi(argv[2])) < 1 ||
	(alpha = atof(argv[3])) < 0 || (beta = atof(argv[4])) < 0 ||
	(Gamma = atof(argv[5])) < 0 || (epsilon = atof(argv[6])) < 0) {
	fprintf(stderr, "usage: [-pgm fname | -pgm2 fname] logn nlinks alpha beta Gamma epsilon\n");
	exit(1);
    }

    srandom(time(NULL));

    if (fabs(alpha+beta+Gamma+epsilon-1) > 0.00001) {
	fprintf(stderr, "error: probabilities sum to %g\n", alpha+beta+Gamma+epsilon);
	exit(1);
    }

    n = 1 << n;
    genlinks = 0;
    nhtab = nlinks*2;

    htab = calloc(nhtab, sizeof(*htab));

    if (pgm_name || pgm2_name) {
        prob = calloc(n*n, sizeof(*prob));
    }

    if (pgm_name) {
        genprob(prob, 1, 1, 0, 0, n, n);
        dumppgm(pgm_name, prob, n);
    }

    fprintf(stderr, "generating %d nodes with %d links (avg degree = %g)\n",
		n, nlinks, (2.0*nlinks)/n);

    while (genlinks < nlinks) {
#if SYMMETRIC
	gensym(n, &u, &v);
#else
	gen(n, &u, &v);
#endif
	if (u == v) continue;
	if (hashadd(u, v)) genlinks++;
    }

    dumplinks();

    if (pgm2_name) {
	recordlinks(prob, n);
	dumppgm(pgm2_name, prob, n);
    }

    return 0;
}