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


/*
 USAGE: 
	 i) compile with: gcc -lm bootstrap.c -o bootstrap  
	 ii) run with: ./bootstrap N p n_runs print_images 
 */



/* first define some random number generators*/
float randomf(void) {
	float x = (float)rand()/(float)(RAND_MAX);
	return x;
}

int randomi(int N) {
	int x = rand()%N;
	return x;
}



/* system initializer */
int initialize(int N, float p, int* grid, int* active_x, int* active_y) {

	int n_active =0;
	int n_alive =p*N*N;
	int i,j;
	while (n_active<n_alive){
		i=randomi(N);
		j=randomi(N);
		while(*(grid+j*N+i)){
			i=randomi(N);
			j=randomi(N);
		}
		*(grid+j*N+i)=1;
		*(active_x+n_active)=i;
		*(active_y+n_active)=j;
		n_active++;
	}

	return n_active;
}

/* collect the coordinates of nearest neighbours */
int get_nbrs(int N, int* nbrs, int x, int y ){
	int n_nbrs=0;
	*(nbrs+n_nbrs)=x-1;
	if (*(nbrs+n_nbrs)<0) *(nbrs+n_nbrs)+=N;
	*(nbrs+n_nbrs+4)=y;
	n_nbrs++;
	*(nbrs+n_nbrs)=(x+1)%N;
	*(nbrs+n_nbrs+4)=y;
	n_nbrs++;
	*(nbrs+n_nbrs)=x;
	*(nbrs+n_nbrs+4)=y-1;
	if (*(nbrs+n_nbrs+4)<0) *(nbrs+n_nbrs+4)+=N;
	n_nbrs++;
	*(nbrs+n_nbrs)=x;
	*(nbrs+n_nbrs+4)=(y+1)%N;
	n_nbrs++;

	return n_nbrs;

}


/* check how many active NNs */
int check_nbrs(int N, int* grid, int x, int y) {
	int* nbrs = calloc(4*2, sizeof(int));
	int n_nbrs = get_nbrs(N, nbrs, x, y);

	int active=0;
	for (int i=0; i<n_nbrs; i++) active+= *(grid+*(nbrs+i)+*(nbrs+i+4)*N);


	free(nbrs);
	return active;
}



/* print the state of systemto stdout - only for testing purposes */
void print_grid(int N, int* grid) {
	for (int j=0; j<N; j++){
                for (int i=0; i<N; i++){
			if (*(grid+i+j*N))
				//printf("%d ", *(grid+i+j*N));
				printf("o ");
			else
				printf(". ");
			}
		printf("\n");
	}
	printf("########################################\n");
}


/* helper function to check if coordinates in array */
int not_in(int* list_x, int* list_y, int x, int y, int n){
	for (int i=0; i<n; i++){
		if (*(list_x+i)==x && *(list_y+i)==y )
			return 0;
	}
	return 1;
}


/* main routine */
int percolation(int N, float p, int print_images) {
        int* grid = calloc(N*N, sizeof(int)); // system state is saved here
        int* active_x = calloc(N*N, sizeof(int)); // x- and y-coordinates of active sites
        int* active_y = calloc(N*N, sizeof(int)); // 
        int* changes_x = calloc(N*N, sizeof(int)); // x- and y- coordinates of the collected sites to be activated
        int* changes_y = calloc(N*N, sizeof(int));
	int* nbrs = calloc(4*2, sizeof(int)); // nearest neighbour coordinates

	int n_active = initialize(N, p, grid, active_x, active_y);
	int step = 0, changes, n_nbrs;

	/* start the evolution of the system */	
	while (n_active>0) {
		if (print_images) print_grid(N, grid);

		/*collect the sites that become active at the current step and save them for the next step  */
		changes = 0;
		for (int i=0; i<n_active;i++ ){
			n_nbrs = get_nbrs(N, nbrs, *(active_x+i), *(active_y+i));
			for (int j=0; j<n_nbrs; j++) {
				if (*(grid+*(nbrs+j)+*(nbrs+j+4)*N)==0  && not_in(changes_x, changes_y, *(nbrs+j),*(nbrs+j+4), changes)  &&  check_nbrs(N,grid,*(nbrs+j),*(nbrs+j+4))>=2 ) {
					*(changes_x+changes) = *(nbrs+j);
					*(changes_y+changes) = *(nbrs+j+4);
					changes++;
				}
			}

		}
		/* then activate the collected sites */
		for (int i=0; i<changes; i++) *(grid+*(changes_x+i)+*(changes_y+i)*N )=1;
		/*here we keep track of the activated sites on the previous step to check only their neighbours on the next step */
		n_active=changes;
		memcpy(active_x, changes_x, sizeof(int)*changes);		
		memcpy(active_y, changes_y, sizeof(int)*changes);		
	}


	int all_active =0;
		
	for (int j=0; j<N; j++){
                for (int i=0; i<N; i++){
			all_active+=*(grid+i+j*N);
		}
	}




	free(grid);
	free(active_x);
	free(active_y);
	free(changes_x);
	free(changes_y);
	free(nbrs);

	/* return 1 if all sites active, 0 otherwise */
	if (all_active==N*N) return 1;
	return 0;
}




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

        srand((unsigned int)time(NULL));


	/* first read the command line parameters */
        if (argc<5) return 0;
        int const N = atoi(argv[1]);
        float const p = atof(argv[2]);
        int const n_runs = atoi(argv[3]);
        int const print_images = atoi(argv[4]);





	int all_infected =0;

	printf("Bootstrap percolation: %d runs of %dx%d with initial p=%f\n", n_runs, N, N, p);
	for (int i=0; i<n_runs; i++) all_infected+= percolation(N, p, print_images);

	printf("    systems entirely infected %d/%d\n", all_infected, n_runs);



	return 0;



}