/*
 ============================================================================
 Name        : furtszam.c
 Author      : Stippinger Marcell
 Version     : 1.0
 Copyright   : BME-TTK FI SzamSzim gyak segedanyag
 Description : Simple cluster counting in C, Ansi-style
 ============================================================================
 */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>

/* Spin value type */
typedef unsigned char byte;
/* Spins */
#define UP 0
#define DOWN 1
#define SEEN 255
/* Lattice dimension */
#define DIMENSION 2

/* Integer power implementation */
int power(int base, int exponent)
{
	int n, product=1;
	for (n=exponent; n>0; n--) { product*=base; }
	return product;
}


/* Initialize a random sample with given occupation probability */
void init_sample(gsl_rng *r, byte *sample, long volume, double probability)
{
	long i; /* index throught the elements of the sample */

	for (i=0; i<volume; i++) {
		sample[i] = (gsl_rng_uniform(r) < probability) ? UP : DOWN;
	}
}


/* Breadth-first search on the lattice */
long breadth_first_search(byte *sample, long volume, long *step, long *queue, long start)
{
	int j;					/* counter over the possible directions */
	long *queue_head,		/* the head of the queue structure		*/
		 *queue_tail,		/* the tail of the queue structure		*/
		 neighbour;			/* index of a neighbouring element		*/
	byte mycolor;			/* the color of the starting cell		*/

	queue_head = queue_tail = queue;	/* initialize pointers		*/
	*queue_head = start;	/* put the starting cell in the queue	*/
	mycolor = sample[start];/* and memorize its color				*/
	
	sample[start] = SEEN;	/* mark the starting cella as visited	*/

	while (queue_head<=queue_tail)		/* while the queue is not empty move to the next cell in the cluster */
	{
		for (j=0; j<2*DIMENSION; j++)	/* look into direction j	*/
		{
			/* The cells have only one composite index, it contains
			 * the coordiantes of all the DIMENSION directions.
			 * The step[] array identifies one of the orthogonal
			 * unit displacements along the base directions.
			 * We use helical BC-s, so we only have to check that we
			 * remain in the bulk sample (i.e. there is no index overrun
			 */
			neighbour = *queue_head + step[j];
			if (neighbour<0) { neighbour+=volume; }
			if (neighbour>=volume) { neighbour-=volume; }

			/* If the neighbour of the same spin, it belongs to the
			 * cluster so we have to remember to treat it
			 */
			if (sample[neighbour]==mycolor)
			{
				*(++queue_tail) = neighbour;	/* put the neighbour cell in the queue	*/
				sample[neighbour] = SEEN;		/* and mark it as visited				*/
			}
		}
		/* Move the head pointer of the queue to the next element in the queue
		 * please note, that we are only changing the variable "queue_head", not "queue"
		 */
		queue_head++;
	}

	return (queue_tail-queue)+1;	/* number of elements put into the queue */
}

/* Count the clusters of UP spins */
long cluster_count(byte* sample, long volume, long *step, long* queue)
{
	long i,					/* index throught the elements of the sample	*/
		 count=0,			/* cluster counter								*/
		 min_size=volume+1,	/* minimum cluster size found					*/
		 max_size=-1;		/* maximum cluster size found					*/
	long cluster_size;		/* the size of the cluster						*/

	/* Search all the lattice for UP spins */
	for (i=0; i<volume; i++) {
		if (sample[i]==UP){
			count++;
			cluster_size = breadth_first_search(sample, volume, step, queue, i);
			if (cluster_size < min_size) { min_size=cluster_size; }
			if (cluster_size > max_size) { max_size=cluster_size; }
		}
	}

	/* Do output */
	printf("%d\t%d\t%d\n", count, min_size, max_size);
	return count;
}

int main(int argc, char* argv[], char* envp[])
{
	int j;					/* step indexing variable						*/
	byte *sample;			/* lattice using helical boundary conditions	*/
	long *queue;			/* queue maximally the same size as the lattice */
	long sidelength = 1000;	/* sidelength of the lattice					*/
	long volume = power(sidelength,DIMENSION);	/* volume of the lattice	*/
	long step[2*DIMENSION];	/* possible steps (composite numbers)			*/
	double p;				/* occupation probability						*/
	const double dp = 0.010;/* delta p										*/

	/* GSL RNG */
	const gsl_rng_type * T;
	gsl_rng * r;
	gsl_rng_env_setup();
	T = gsl_rng_mt19937;
	r = gsl_rng_alloc (T);

	/* Create composite indices for all possible step directions */
	for (j=0; j<DIMENSION; j++) {
		step[2*j] = power(sidelength,j);
		step[2*j+1] = -power(sidelength,j);
	}

	/* Allocate the needed space */
	sample = calloc(volume, sizeof(byte));
	queue = calloc(volume, sizeof(long));

	/* Sweep the occupation probability */
	printf("#p\tcount\tmin_size\tmax_size\n");
	for (p=0; p<=1; p+=dp) {
		init_sample(r, sample, volume, p);
		printf("%lf\t",p);
		cluster_count(sample, volume, step, queue);
	}

	/* Free the memory */
	free(sample);
	free(queue);
	gsl_rng_free (r);

	return 0;
}