/*
 * This file is Appendix C of the paper:
 *   Populations on Fragmented Landscapes with Spatially Structured
 *   Heterogeneities: Landscape Generation and Local Dispersal
 *     by David Hiebeler
 *     Ecology (2000), vol. XX, pp. YY-ZZ.
 *
 *
 *
 * File: 2x1gen.c
 *   By: Dave Hiebeler
 *       hiebeler@cam.cornell.edu
 *       Original version: September 1996
 *       Revised for publication: December 1998
 *
 * This program generates a heterogeneous landscape on a rectangular lattice,
 * which has a given probability distribution of 2x1 blocks for the
 * two habitat types.
 *
 * By default, a text picture of the landscape is output to the file
 * "2x1gen.out", with the characters `.' and `X' representing habitat
 * types 0 and 1, respectively.
 *
 * If you use the "-pgm" argument, a Portable Graymap image will be output
 * instead.  There is a large collection of free utilities written by
 * Jef Poskanzer which manipulates files of this format and converts them to
 * many other image formats (such as GIF); these utilities come for free with
 * many distributions of Linux, or you can download them from:
 *   http://www.acme.com/software/pbmplus/
 *
 * If you run the program with no command-line arguments, information
 * about other command-line options will be displayed.
 *
 * EXAMPLE:
 * A simple example of how to use the program is: to generate a 100x100
 * landscape with p0=0.5 and q00=0.8 (see the paper for a description of
 * these parameters), and store the output in the default file "2x1gen.out",
 * you would run:
 *   2x1gen 100 100 0.5 0.8
 * To do the same thing, but generate a PGM file called "land.pgm", you
 * would run:
 *   2x1gen -pgm -out land.pgm 100 100 0.5 0.8
 *
 * Unless you use the "-q" option for "quiet mode", the program displays
 * some diagnostics about the parameters you specified, and how close
 * the landscape came to achieving the desired structure.
 *
 * On my Linux system, I compile this program as follows:
 *   gcc -O3 -o 2x1gen 2x1gen.c
 * (The "-O3" flag tells the compiler to optimize the code's performance
 * as much as possible.)
 *
 *
 * Portability concerns:
 * This program was developed on a machine running Linux.  It should
 * be very portable, with two possible exceptions:
 *
 * (1) The main thing that may be a problem is the "my_srandom()"
 *     function, which will probably need to be changed for non-Unix systems.
 *     Simply change it to call whatever the appropriate initialization routine
 *     is for your random number generator.  Note that my function generates
 *     a seed based on the current time, and then calls the system routine
 *     "srandom()" with that seed.  I also call "srandom()" from "main()",
 *     if the user supplies a seed, so you'll need to change that as well.
 *
 * (2) The random number generator itself may be different if you are not
 *     using Unix.  The random number generator I use is called "random()",
 *     which on my computer (running Linux) returns a random 32-bit integer.
 *     I use the "random()" routine in my "randFrac()" function, and also
 *     in the "gen2x1()" function to choose random coordinates of a site
 *     to try flipping.
 */


#include <stdio.h>
#include <math.h>
#include <sys/time.h>  /* for seeding the random number generator */

#define RAND01_PRECISION 100000 /* granularity for generating random numbers */
#define EPSILON 1.0e-5 /* used to test whether or not a number is zero */
#define DEFAULT_OUTFNAME "2x1gen.out" /* default output filename */
#define DEFAULT_MAXITERS 1000000 /* default maximum number of iterations */
#define DEFAULT_MAXERR 0.001 /* default tolerance between desired and measured
			      * probabilities */
#define DEFAULT_CHAR0 '.'
#define DEFAULT_CHAR1 'X'


typedef unsigned char byte;  /* just for convenience */


/*
 * Print an error message and exit
 */
void
quit(char *s)
{
    fprintf(stderr, "%s\n", s);
    exit(1);
}


/*
 * Print a usage message and exit.
 */
void
printusage(char *s)
{
    fprintf(stderr, "Usage: %s [-byte val0 val1 | -char char0 char1 | -pgm]\n", s);
    fprintf(stderr, "            [-out fname] [-q] [-iter maxiters maxerr] [-s seed]\n");
    fprintf(stderr, "            xsize ysize p0 q00\n", s);
    fprintf(stderr, "   -byte val0 val1: numerically specify the 2 values to output\n");
    fprintf(stderr, "   -char char0 char1: use characters to specify the 2 values to output\n");
    fprintf(stderr, "                      (default: the characters `.' and `X')\n");
    fprintf(stderr, "   -s seed: specify random number seed\n");
    fprintf(stderr, "   -out fname: specify where to store the results (use `-' for stdout)\n");
    fprintf(stderr, "               (default: `2x1gen.out')\n");
    fprintf(stderr, "   -q: quiet mode\n");
    fprintf(stderr, "   -iter maxiters maxerr: specify maximum number of iterations, and max error\n");
    fprintf(stderr, "              (default: 1000000 .001)\n");
    fprintf(stderr, "  xsize and ysize specify the size of the landscape\n");
    fprintf(stderr, "  p0 and q00 are the landscape parameters as described in the paper\n");
    exit(2);
}


/*
 * Seed the random number generator, based on the current time in
 * milliseconds.  This will probably need to be rewritten for non-Unix
 * systems.
 */
void
my_srandom()
{
    struct timeval tp;
    struct timezone *tzp;
    FILE *fp;

    tzp = NULL;
    gettimeofday(&tp, tzp);
    srandom((int) tp.tv_usec);
    if ((fp=fopen(".2x1gen.srandom", "w"))==NULL)
	return;
    fprintf(fp, "%d\n", tp.tv_usec);
    fclose(fp);
}


/*
 * Return 1 a specified fraction of the time, else return 0.
 */
int
randFrac(double frac)
{
    int i;
    double d;

    /* On the system I use, the random() function returns a random
     * 32-bit integer.  If you are using a computer which has a
     * function that returns a random floating-point number between
     * 0 and 1, then you can just use that function to set the value
     * of the variable "d" directly, which is basically what I achieve
     * below.
     */
    i = random() % RAND01_PRECISION;
    d = ((double) i) / (double)RAND01_PRECISION;
    if (d < frac)
	return 1;
    else
	return 0;
}


/*
 * Check to see whether a number is close enough to zero that we really
 * should call it zero.
 */
int
fis_zero(double x)
{
    if (fabs(x) > EPSILON)
	return 0;
    else
	return 1;
}


/*
 * Measure the 2x1 block counts in the lattice.
 */
void
measureBlockCounts(int counts[4],  /* array for returning the counts */
		   byte *lattice,  /* pointer to the landscape array */
		   int xsize, int ysize)  /* size of the landscape array */
{
    int i, x, y, x1, y1;

    /* First initialize the counts to zero */
    for (i=0; i < 4; i++)
	counts[i] = 0;
    
    for (y=0; y < ysize; y++)
	for (x=0; x < xsize; x++) {
	    /* north */
	    x1 = x;
	    y1 = (y-1+ysize)%ysize;
	    counts[lattice[y*xsize+x]*2 + lattice[y1*xsize+x1]]++;

	    /* south */
	    x1 = x;
	    y1 = (y+1)%ysize;
	    counts[lattice[y*xsize+x]*2 + lattice[y1*xsize+x1]]++;

	    /* east */
	    x1 = (x+1)%xsize;
	    y1 = y;
	    counts[lattice[y*xsize+x]*2 + lattice[y1*xsize+x1]]++;

	    /* west */
	    x1 = (x-1+xsize)%xsize;
	    y1 = y;
	    counts[lattice[y*xsize+x]*2 + lattice[y1*xsize+x1]]++;
	}
}


/*
 * Compute how the block counts would be affected if we were
 * to flip the state of the cell at position (x,y)
 */
void
computeCountChanges(int changes[4], /* array to return the vector of changes */
		    byte *lattice,  /* landscape array */
		    int xsize, int ysize, /* size of landscape */
		    int x, int y)  /* which site we want to flip */
{
    int x1, y1, i;

    /* Clear out array just to be safe */
    for (i=0; i < 4; i++)
	changes[i] = 0;

    /* north */
    x1 = x;
    y1 = (y-1+ysize)%ysize;
    /* First decrease the appropriate count for the current configuration.
    * Note we always increase or decrease by 2, because every time we make
    * a flip, it will change the 2x1 block as seen from both ends, so to
    * speak, so we need to count the change twice. */
    changes[lattice[y*xsize+x]*2 + lattice[y1*xsize+x1]] -= 2;
    /* And then increase the count for the new configuration */
    changes[(1-lattice[y*xsize+x])*2 + lattice[y1*xsize+x1]] += 2;

    /* south */
    x1 = x;
    y1 = (y+1)%ysize;
    changes[lattice[y*xsize+x]*2 + lattice[y1*xsize+x1]] -= 2;
    changes[(1-lattice[y*xsize+x])*2 + lattice[y1*xsize+x1]] += 2;

    /* east */
    x1 = (x+1)%xsize;
    y1 = y;
    changes[lattice[y*xsize+x]*2 + lattice[y1*xsize+x1]] -= 2;
    changes[(1-lattice[y*xsize+x])*2 + lattice[y1*xsize+x1]] += 2;

    /* west */
    x1 = (x-1+xsize)%xsize;
    y1 = y;
    changes[lattice[y*xsize+x]*2 + lattice[y1*xsize+x1]] -= 2;
    changes[(1-lattice[y*xsize+x])*2 + lattice[y1*xsize+x1]] += 2;

    /* Now handle symmetry, lumping together the changes for [01] and [10]
    * blocks, since they are really the same by symmetry assumptions. */
    i = changes[1] + changes[2];
    if (i%2)
	quit("computeCountChanges(): Hey, 'i' was odd!!! That shouldn't be possible!");

    changes[1] = changes[2] = i/2;
}


/*
 * Given the parameters p0 and q00, compute the 2x1 block
 * probabilities, p00, p01, and p11.
 */
void
condProbsTo2x1Probs(double p0, double q00,  /* input parameters */
		    double *p00, double *p01, double *p11)  /* return parms */
{
    *p00 = p0 * q00;				/* p00 */
    *p01 = p0 - *p00;				/* p01 */
    *p11 = 1.0 - (*p00 + 2.0 * (*p01));		/* p11 */

    /* Sometimes we end up with a parameter which is 1.0e-10, i.e. zero for
     * all practical purposes.  Here we set such values to actually BE zero,
     * for simplicity */
    if (fis_zero(*p00))
	*p00 = 0.0;
    if (fis_zero(*p01))
	*p01 = 0.0;
    if (fis_zero(*p11))
	*p11 = 0.0;
    return;
}


/*
 * This is the heart of the program.  It generates the landscape with
 * spatially structured heterogeneities.
 */
void
gen2x1(byte *lattice,  /* landscape array */
       int xsize, int ysize,  /* size of landscape */
       double p0, double q00,  /* the main parameter values */
       int quiet,   /* whether or not to work quietly */
       int iterates, /*maximum number of times to loop through the algorithm*/
       double maxErr)  /* the maximum difference between measured and desired
			* probabilities that we will accept when stopping */
{
    int
	x, y, x1, y1, i, j,
	blockCounts[4],  /* used to count the numbers of occurances of the
			  * various possible 2x1 blocks on the landscape */
	tmpCounts[4],  /* a temporary array of counts */
	targetCounts[4],  /* the counts we would see if we had the landscape
			   * we are trying to get */
	countChanges[4],  /* for computing changes to counts */
	oldErr, newErr,  /* differences between desired and achieved COUNTS */
	timesFlipped,  /* for counting how many times we flipped sites */
	stop;  /* flag for stopping the algorithm */
    double p00, p01, p11,  /* the 2x1 block probabilities */
	p1,
	probErr, tmpProbErr; /* differences between desired and
			      * achieved PROBABILITIES */

    timesFlipped = 0;

    /* Compute the 2x1 block probabilities */
    condProbsTo2x1Probs(p0, q00, &p00, &p01, &p11);
    p1 = 1.0 - p0;

    /* First fill the lattice randomly, using the correct patch
     * occupancy probability */
    bzero(lattice, xsize*ysize);  /* Clear out the lattice */
    for (y=0; y < ysize; y++)
	for (x=0; x < xsize; x++)
	    if (randFrac(p1))
		lattice[y*xsize + x] = 1;

    /* Now measure the block counts from the lattice */
    measureBlockCounts(blockCounts, lattice, xsize, ysize);

    /* Compute target block counts from the block probabilities.
     * The reason for the factor of 4 is because when we count blocks,
     * we check 4 neighbors of each site.
     */
    targetCounts[0] = (int)(p00 * xsize*ysize*4);
    targetCounts[1] = targetCounts[2] = (int)(p01 * xsize*ysize*4);
    targetCounts[3] = (int)(p11 * xsize*ysize*4);

    if (quiet == 0) {  /* Print out some diagnostic info */
	printf("Target counts are ");
	for (j=0; j < 4; j++)
	    printf("%d  ", targetCounts[j]);
	printf("  (sum = %d)", targetCounts[0]+targetCounts[1]+
	       targetCounts[2]+targetCounts[3]);
	printf("\n");

	printf("Starting counts are ");
	for (j=0; j < 4; j++)
	    printf("%d  ", blockCounts[j]);
	printf("  (sum = %d)", blockCounts[0]+blockCounts[1]+
	       blockCounts[2]+blockCounts[3]);
	printf("\n");
    }


    oldErr = 0;
    /* Initialize the discrepancy between desired and achieved counts */
    for (i=0; i < 4; i++)
	oldErr += abs(targetCounts[i] - blockCounts[i]);

    /* Now comes the main loop.  We repeatedly pick a site at random, and
     * see if flipping its value (0 to 1, or 1 to 0) would move the block
     * counts (or equivalently, the 2x1 block probabilities) closer to the
     * desired values.  If so, flip the site.  Otherwise, try again. */
    stop = 0;
    for (i=0; (i < iterates) && (stop == 0); i++) {
	for (j=0; j < 4; j++)  /* copy the current counts into a temp. array */
	    tmpCounts[j] = blockCounts[j];

	/* Pick a random site.
	 * If your random number generator returns a number between 0 and 1,
	 * you will need to do something like this [say your random number
	 * generator is called rand01() ]:
	 *   x = (int)(rand01() * xsize);
	 *   y = (int)(rand01() * ysize);
	 * That assumes rand01() returns a number greater than or equal
	 * to zero, but strictly less than 1.
	 */
	x = random()%xsize;
	y = random()%ysize;

	/* now see how the counts would be affected if we were to
	 * change the state of this cell */
	computeCountChanges(countChanges, lattice, xsize, ysize, x, y);

	/* Merge the changes back into the count array */
	for (j=0; j < 4; j++)
	    tmpCounts[j] += countChanges[j];

	/* And compute what the new discrepancy between desired and achieved
	 * counts would be if we make this flip */
	newErr = 0;
	for (j=0; j < 4; j++)
	    newErr += abs(targetCounts[j] - tmpCounts[j]);

	/* if the new error would be smaller, flip the site's value */
	if (newErr <= oldErr) {
	    lattice[y*xsize+x] = 1-lattice[y*xsize+x];
	    oldErr = newErr;
	    /* Move the temporary copy of the counts into the regular array */
	    for (j=0; j < 4; j++)
		blockCounts[j] = tmpCounts[j];
	    timesFlipped++;  /* keep track of how many flips we make */
	}

	/* Convert the discrepancy in number of blocks into discrepancy
	 * in probabilities, by normalizing */
	probErr = ((double)oldErr)/(double)(xsize*ysize*4.0);

	/* And if the error is small enough, stop. */
	if (probErr < maxErr)
	    stop = 1;
    }

    if (quiet == 0) {
	printf("Achieved counts are ");
	for (j=0; j < 4; j++)
	    printf("%d  ", blockCounts[j]);
	printf("  (sum = %d)", blockCounts[0]+blockCounts[1]+
	       blockCounts[2]+blockCounts[3]);
	printf("\n");

	printf("Changed site's value a fraction %lg of the time, after %d iterations\n",
	       ((double)timesFlipped)/(double)i, i);
    }
}


/*
 * Main() routine.
 */
main(int argc, char **argv)
{
    int xsize, ysize,  /* size of landscape */
	iterates=DEFAULT_MAXITERS,  /* maximum number of iterations */
	quiet,  /* flag for controlling output */
	x, y, i, tmpInt,
	asciiFlag=1,  /* whether or not we should output the landscape in
		      * a human-readable format */
	pgmFlag=0,
	blockCounts[4],
	rseed, gotSeed=0;
    double p00, p01, p11,  /* the 2x1 block probabilities */
	p0, q00,  /* marginal and conditional probability parameters */
	dsum, actualProbs[4], actualp0, actualq00,
	maxErr=DEFAULT_MAXERR;  /* maximum discrepancy in probabilities that
				 * we will accept */
    char outFname[256],  /* name of output file */
	char0=DEFAULT_CHAR0, char1=DEFAULT_CHAR1;  /* output characters */
    FILE *fp;
    byte *lattice;  /* the main landscape array */

    quiet=0;
    sprintf(outFname, DEFAULT_OUTFNAME);


    /*
     * Parse command-line arguments
     */
    i = 1;
    while ((i < argc) && (argv[i][0] == '-')) {
	if (!strcmp(argv[i], "-byte")) {
	    if (i >= argc-2)
		printusage(argv[0]);
	    if (sscanf(argv[i+1], "%d", &tmpInt) != 1)
		printusage(argv[0]);
	    char0 = (char)tmpInt;
	    if (sscanf(argv[i+2], "%d", &tmpInt) != 1)
		printusage(argv[0]);
	    char1 = (char)tmpInt;
	    asciiFlag = 0;
	    i += 3;
	}
	else if (!strcmp(argv[i], "-char")) {
	    if (i >= argc-2)
		printusage(argv[0]);
	    if (strlen(argv[i+1]) != 1)
		printusage(argv[0]);
	    char0 = argv[i+1][0];
	    if (strlen(argv[i+2]) != 1)
		printusage(argv[0]);
	    char1 = argv[i+2][0];
	    i += 3;
	}
	else if (!strcmp(argv[i], "-pgm")) {
	    pgmFlag = 1;
	    asciiFlag = 0;
	    char0 = 0;
	    char1 = 1;
	    i++;
	}
	else if (!strcmp(argv[i], "-s")) {
	    gotSeed = 1;
	    if (i >= argc-1)
		printusage(argv[0]);
	    if (sscanf(argv[i+1], "%d", &rseed) != 1)
		printusage(argv[0]);
	    srandom(rseed);
	    i += 2;
	}
	else if (!strcmp(argv[i], "-out")) {
	    if (i >= argc-1)
		printusage(argv[0]);
	    strncpy(outFname, argv[i+1], 255);
	    i += 2;
	}
	else if (!strcmp(argv[i], "-q")) {
	    quiet = 1;
	    i++;
	}
	else if ((!strcmp(argv[i], "-iter")) || (!strcmp(argv[i], "-iters"))) {
	    if (i >= argc-2)
		printusage(argv[0]);
	    if (sscanf(argv[i+1], "%d", &iterates) != 1)
		printusage(argv[0]);
	    if (sscanf(argv[i+2], "%lg", &maxErr) != 1)
		printusage(argv[0]);
	    i += 3;
	}
	else
	    printusage(argv[0]);
    }
    if (gotSeed == 0)
	my_srandom();
    if (argc != i+4)
	printusage(argv[0]);
    if (sscanf(argv[i++], "%d", &xsize) != 1)
	printusage(argv[0]);
    if (sscanf(argv[i++], "%d", &ysize) != 1)
	printusage(argv[0]);
    if (sscanf(argv[i++], "%lg", &p0) != 1)
	printusage(argv[0]);
    if (sscanf(argv[i++], "%lg", &q00) != 1)
	printusage(argv[0]);

    /* Compute 2x1 block probabilities */
    condProbsTo2x1Probs(p0, q00, &p00, &p01, &p11);

    if (quiet == 0) {
	printf("p0 = %lg, q00 = %lg\n", p0, q00);
	printf("p00 = %lg, p01 = %lg, p11 = %lg\n", p00, p01, p11);
    }

    /* Check to make sure there aren't any problems with the parameters
     * entered (even problems that may come up due to round-off). */
    if ((p00 < 0.0) || (p01 < 0.0) || (p11 < 0.0) ||
	(p00 > 1.0) || (p01 > 1.0) || (p11 > 1.0))
	quit("Invalid values for p0 and q00");
    if ((p00 + 2.0*p01 - 1.0 > EPSILON) || (p11 + 2.0*p01 - 1.0 > EPSILON))
	quit("Probabilities were too large!");

    /* Allocate memory for the landscape array */
    if ((lattice=(byte *)malloc(xsize*ysize))==NULL)
	quit("Couldn't malloc lattice");

    /* Open output file */
    if (!strcmp(outFname, "-"))
	fp=stdout;
    else
	if ((fp=fopen(outFname, "w"))==NULL)
	    quit("Couldn't open 2x1gen.out");


    /* Now generate the landscape. */
    gen2x1(lattice, xsize, ysize, p0, q00, quiet, iterates, maxErr);

    /* And output it. */
    /* First I output a standard header, which my simulation program
     * will recognize, to make sure I use the right kind of landscape files
     * with the simulations... */
    if (asciiFlag) {
	fprintf(fp, "ascii landscape v0.2\n");
	fprintf(fp, "xsize=%d ysize=%d 0char=%c 1char=%c\n", xsize, ysize,
		char0, char1);
    }
    else if (pgmFlag) {
	fprintf(fp, "P5\n");  /* header for PGM file */
	fprintf(fp, "%d %d\n", xsize, ysize);  /* image size */
	fprintf(fp, "1\n");   /* maximum gray value in the image */
    }
    else {
	fprintf(fp, "landscape v0.2\n");
	fprintf(fp, "xsize=%d ysize=%d 0byte=%d 1byte=%d\n", xsize, ysize,
		(int)char0, (int)char1);
    }
    for (y=0; y < ysize; y++) {
	for (x=0; x < xsize; x++)
	    if (lattice[y*xsize+x] == 0)
		putc(char0, fp);
	    else
		putc(char1, fp);
	if (asciiFlag)
	    putc('\n', fp);
    }
    fclose(fp);

    /* Now compute some diagnostics */
    measureBlockCounts(blockCounts, lattice, xsize, ysize);

    /* Compute the 2x1 block probabilities in the landscape we generated */
    for (i=0; i < 4; i++)
	actualProbs[i] = ((double)blockCounts[i])/(xsize*ysize*4.0);
    actualp0 = actualProbs[0] + actualProbs[1];
    if (actualp0 == 0.0)
	actualq00 = 0.0;
    else
	actualq00 = actualProbs[0]/actualp0;  /*  q00 = p[00]/p0  */

    if (quiet == 0) {
	printf("=====Desired 2x1 block probabilities=====\n");
	printf("p[00]=%lg  p[01]=%lg  p[11]=%lg\n", p00, p01, p11);
	printf("=====Actual 2x1 block probabilities=====\n");
	printf("p[00]=%lg  p[01]=%lg  p[11]=%lg\n", actualProbs[0],
	       actualProbs[1], actualProbs[3]);
	printf("\n");
    }

    if (quiet == 0) {
	printf("====Desired marginal and conditional probabilities====\n");
	printf("p0 = %lg, q00 = %lg\n", p0, q00);
	printf("====Actual marginal and conditional probabilities====\n");
	printf("p0 = %lg, q00 = %lg\n", actualp0, actualq00);
    }

    dsum = fabs(actualProbs[0] - p00) + fabs(actualProbs[1] - p01) +
	fabs(actualProbs[2] - p01) + fabs(actualProbs[3] - p11);

    if (quiet == 0)
	printf("1-norm of difference in 2x1 block probs is %lg\n", dsum);
    exit(0);
}