#include <stdlib.h>
#include <stdio.h>
 
 
int i, j, n;
int matrix[100][5] = {0};                              //Define a matrix with a maximum of 100 rows and 6 columns
int AA=0, AB=0,BB=0;                                  // Define the genotype counter     
int repeats = 0;
int norepeats;
float p, q;
float rand1;
float pp, freq;
void SetSeed();
 
int main () {
 
//FILE *output_file = fopen("output.txt", "w");
 
	printf("Enter the frequency of allele A: ");
	scanf("%f", &p);
	printf("\n");
	pp=p*10000;
	freq=0.5*10000;
	printf("Frequency of allele A (p) is: %.3f\n", p);       // States the frequency of the A allele
	q = 1 - p;                                               // The value of q is 1 - (the probability entered for p)
 
	printf("Frequency of allele B (q) is: %.3f\n", q);       // States the frequency of the B allele
	printf("\n");
	printf("How many individuals in the genealogy?");
	scanf("%d", &n);
 
	printf("\n");
	printf("Please enter the desired number of repeats: ");
	scanf("%d", &norepeats);                                 // Desired number of repeats of the program stored in 'norepeats'
	printf("\n");
 
	SetSeed(); 												 // Seed file for RNG set here		
	const float RANGE = 10000;                              // Random number generated is below 1
	rand1 = (float) rand() * RANGE / RAND_MAX;
 
	for(j=1;j<norepeats;j++){
 
		FILE *input_txt;        
 
		input_txt = fopen("input.txt","r");               
 
		for(i = 0; i < n; i++){
 
			fscanf(input_txt, "%d %d %d %d %d", &matrix[i][0], &matrix[i][1], &matrix[i][2], &matrix[i][3], &matrix[i][4]);
			if (feof(input_txt)) break;
		}
 
 
		for(i=0;i<n;i++){
			if (matrix[i][1] == 0 && matrix[i][2] == 0){             // If the individual is a founder:
 
				rand1 = (float) rand() * RANGE / RAND_MAX;            // Random number 1 generated 
				if (rand1 <= pp){matrix[i][3] = 0;}                // Place a 0 allele down in column 4 if random number 1 is below the set p value    
				else matrix[i][3] = 1;                    // Place a 1 allele down in column 4 if random number 1 is above the set p value
 
				rand1 = (float) rand() * RANGE / RAND_MAX;            // Random number 2 generated
				if (rand1 <= pp){matrix[i][4] = 0;}                                  // Place a 0 allele down in column 5 if random number 2 is below the set p value
				else matrix[i][4] = 1;                                         // Place a 1 allele down in column 5 if random number 2 is above the set p value
 
			}//end of if
 
		else{
			int dad = matrix[i][1];                               // 1 is "dad" column
			int mum = matrix[i][2];                               // 2 is "mum" column
 
			int dadRowIndex = dad - 1;                            // The row of the father is 'dad - 1' as the first row is set at zero
			int mumRowIndex = mum -1;                             // The row of the mother is 'mum - 1' as the first row is set at zero
 
 
			rand1 = (float) rand() * RANGE / RAND_MAX;            // Random number 3 generated
			if (rand1 <= freq){matrix[i][3] = matrix[dadRowIndex][3];}  // If random number 3 is less than a half, receive allele stored in column 4 from dad
			else{matrix[i][3] = matrix[dadRowIndex][4];}             // If random number 3 is more than a half, receive allele stored in column 5 from dad
 
 
			rand1 = (float) rand() * RANGE / RAND_MAX;            // Random number 4 generated
			if (rand1 <= freq){ matrix[i][4] = matrix[mumRowIndex][3];}        // If random number 4 is less than a half, receive allele stored in column 4 from mum
			else{matrix[i][4] = matrix[mumRowIndex][4];}            // If random number 4 is more than a half, receive allele stored in column 5 from mum
 
			}//end of else
 
            if (matrix[i][3] == 0 && matrix[i][4] == 0)   // Counting the number of times AA appear in the 1000 simulations.
     	AA++;
     if (matrix[i][3] == 0 && matrix[i][4] == 1)   // Counting the number of times AB appear in the 1000 simulations.
     	AB++;
     if (matrix[i][3] == 1 && matrix[i][4] == 0)
     AB++; 
     if (matrix[i][3] == 1 && matrix[i][4] == 1)   // Counting the number of times BB appear in the 1000 simulations.
     	BB++;
 
		}//end of for
 
 
 
	}//end of for
 
 
	printf("the matrix is:\n");
	for(i=0;i<n;i++){
		for(j=0;j<5;j++){
		printf("%d\t", matrix[i][j]);
		}
	printf("\n");
	}
 
	printf("\n");
	printf("Cumulative genotype totals:\n");        
	printf("AA: %d\n", AA);                      // Print total number of AA individuals
	printf("AB: %d\n", AB);                      // Print total number of heterozygotes
	printf("BB: %d\n", BB);                      // Print total number of BB individuals
	printf("\n");
	printf("\n");
 
 
 
			}//end of main
 
	void SetSeed() 
	{
                                                           // Determining the seed file
     FILE *seed_file;
     long int seed_value;
 
     seed_file = fopen("seed.txt","r");                   // Open seed file
     fscanf(seed_file,"%d",&seed_value);                  // Exctract seed from file
     fclose(seed_file);                                   // Close seed file
     srand(seed_value);                                   // Randomise seed value
     seed_file = fopen("seed.txt","w");                   // Create temporary output
     fprintf(seed_file,"%d",rand());                      // Print random number to output
     fclose(seed_file);                                   // Close temporary seed output file
	}