/* simmss.c */
/* T J Finney */
/* use as you please */

#include <stdio.h>
#include <stdlib.h>         /* Define the rand(), srand() and calloc() functions */
#include <time.h>           /* Define the time() function */

#define	COPY		10		/* Inverse probability of a copy (years) */
#define	ERROR		20		/* Inverse probability of an error (variants copied) */
#define	CORRECT		40		/* Inverse probability of a correction (variants copied) */
#define	DESTROY		50		/* Inverse probability of MS destruction (years) */
#define	AGEMAX		100		/* Maximum MS age (years) */
#define	YEARMIN		70		/* Minimum copying year */
#define	YEARMAX		300		/* Maximum copying year */
#define	MSNUM		10		/* Maximum number of MSS to be chosen */
#define	MSTOT		100		/* Maximum total number of MSS */
#define	CELLNUM		20		/* Number of variation cells */

unsigned long int seed;

int random_number(int max)
{
	/* Generate a random number in the range 1..max */

	/* Change seed */
	seed = seed * 9973 + 4999;		/* Pick a seed */
    srand(seed);					/* Seed number generator */
    return rand() % max + 1;		/* Generate the number */
}

int make_copy(int date, int np, int nd, int *pp, int *dp, int *vp)
{
	/* Copy a parent manuscript to a daughter manuscript */
	
	int i;

	*dp = 1;							/* Is daughter MS alive? 1 = yes, 0 = no */
	*(dp + 1) = date;					/* Year copy made */
	*(dp + 2) = random_number(AGEMAX);	/* MS lifespan */
	*(dp + 3) = np;						/* Parent MS number */
	*(dp + 4) = nd;						/* Daughter MS number */

	for (i = 5; i < CELLNUM; ++i)
	{
		if (random_number(ERROR) == 1)
		{
			*(dp + i) = random_number(*(vp + i) + 1);	/* New variant */
		}
		else
		{
			*(dp + i) = *(pp + i);		/* Copy from parent to daughter */
		}
		if (*(dp + i) > *(vp + i))
		{
			++*(vp + i);				/* Increment appropriate variant cell */
		}
		if (random_number(CORRECT) == 1)
		{
			*(dp + i) = random_number(*(vp + i))/2 + 1;	/* Change to an existing variant */
		}

	}
}

void print_ms(int *msp)
{
	/* Print a MS */
	
	int i;

	for (i = 0; i < CELLNUM; ++i)
	{
		printf("%3d ", *(msp + i));	/* Print MS cell i */
	}
	printf("\n");
}

void file_ms(FILE *filep, int *msp)
{
	/* Print MS to a file */
	
	int i;

	for (i = 0; i < CELLNUM; ++i)
	{
		fprintf(filep, "%3d ", *(msp + i));	/* Print MS cell i to a file */
	}
	fprintf(filep, "\n");
}

main()
{
	/* Declare */
	
	int *ms_array;			/* Variant register and MS array */
	int *ptr;				/* Pointer for MS array */
	int *lptr;				/* Pointer to latest MS */
	int year;				/* Year counter */
	int latest_number;		/* Current total number of mss */
	int start_number;		/* Initial number of mss for each year */
	int i;					/* MS cell counter */
	int n, x;				/* MS number counters */
	int m;					/* Living MS number counter */
	FILE *fp;

	/* Allocate memory for ms_array */
	
	ms_array = calloc(MSTOT + 1, CELLNUM * sizeof(int));
	
	/* Initialise variant register */
	
	for (ptr = ms_array; ptr < ms_array + CELLNUM; ++ptr)
	{
		*ptr = 1;
	}
	
	/* Initialise random number generator seed */
	
	seed = time(NULL);

	/* Make MS 1 */

	ptr = ms_array + CELLNUM;			/* Point to first cell of MS 1 */
	*ptr = 1;							/* Is MS 1 extant? 1 = yes */
	*(ptr + 1) = YEARMIN;				/* Year copy made */
	*(ptr + 2) = random_number(AGEMAX);	/* MS lifespan */
	*(ptr + 3) = 1;						/* Parent MS number (Original) */
	*(ptr + 4) = 1;						/* Daughter MS number (Original) */
	for (i = 5 ; i < CELLNUM; ++i)
	{
		*(ptr + i) = 1;					/* Set rest of cells to 1 */
	}
	latest_number = 1;					/* Number of MSS is now 1 */

	/* Open history file */

	fp = fopen("History", "w");	/* Open file */

	/* Execute loop until copying ends */
	year = YEARMIN;
	while (year <= YEARMAX && latest_number <= MSTOT)
	{
		/* Terminate appropriate MSS */
		for (n = 1; n <= latest_number; ++n)		/* Look at existing MSS */
		{
			ptr = (ms_array + n * CELLNUM);		/* Points to cell 1 of MS n */
			if (*ptr == 1 && random_number(DESTROY) == 1)
			{
				*ptr = 0;
				fprintf(fp, "MS %d killed in %d\n", n, year);
			}
			if (*ptr == 1 && (*(ptr + 1) + *(ptr + 2) <= year))
			{
				*ptr = 0;
				fprintf(fp, "MS %d died in %d\n", n, year);
			}
		}
		
		/* Execute copying loop for existing MSS */
		start_number = latest_number;
		for (n = 1; n <= start_number; ++n)	/* Look at existing MSS */
		{
			ptr = (ms_array + n * CELLNUM);					/* Points to cell 1 of MS n */
			lptr = (ms_array + latest_number * CELLNUM);	/* Points to cell 1 of latest MS */

			/* Select appropriate manuscripts and call copy function */
			if (*ptr == 1 && (random_number(COPY) == 1) && (latest_number <= MSTOT))
			{
				make_copy(year, n, latest_number, ptr, lptr, ms_array);
				file_ms(fp, lptr);
				++latest_number;
			}
		}
		++ year;
	}

	fclose(fp);

	/* Count number of living MSS */
	m = 0;
	for (n = 1; n < latest_number; ++n)
	{
		ptr = (ms_array + n * CELLNUM);				/* Points to cell 1 of MS n */
		if (*ptr == 1)
		{
			++m;
		}
	}

	
	/* Open extant file */
	
	fp = fopen("Extant", "w");	/* Open file */
	
	/* Print year and MSS numbers reached */
	
	printf("  At the year: %d.    Generated MSS: %d.    Living MSS: %d.\n\n",
	year - 1, latest_number - 1, m);
	fprintf(fp, "  At the year: %d.    Generated MSS: %d.    Living MSS: %d.\n\n",
	year - 1, latest_number - 1, m);

	/* Print variant register */
	print_ms(ms_array);
	file_ms(fp, ms_array);

	/* Select the appropriate number of MSS and write to a file */
	if (m > MSNUM)
	{
		/* Choose appropriate number of MSS */
		m = 0;
		while (m < MSNUM)
		{
			x = random_number(latest_number);
			ptr = (ms_array + x * CELLNUM);			/* Points to cell 1 of MS x */
			if (*ptr == 1)
			{
				print_ms(ptr);
				file_ms(fp, ptr);
				*ptr = 0;		/* Kill the MS so it will not be chosen again */
				++m;
			}
		}
	}
	else
	{
		for (n = 1; n <= latest_number; ++n)
		{
			ptr = (ms_array + n * CELLNUM);			/* Points to cell 1 of MS n */
			if (*ptr == 1)
			{
				print_ms(ptr);
				file_ms(fp, ptr);
			}
		}
	fclose(fp);
	free(ms_array);	/* Free allocated memory */
	}
}