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

/*

CUE:  Calculating Undetected Errors
Joshua Tauberer, Derek Gordon, Jurg Ott

To find the proportion of errors that go unnoticed in n-tuples
of genomic data of nuclear families.

Note that allele combinations are stored as codes.  For example:
 	1=1/1, 2=1/2, and 3=2/2 for a two-allele locus.

Please refer to the CUE website for more information:
http://linkage.rockefeller.edu/derek/cue

*/


#define MIN(a,b) ((a>b) ? b : a)
#define MAX(a,b) ((a>b) ? a : b)
#define MINMAXIDX(a,b) [(a>b) ? b : a][((a>b) ? a : b)]

/* STRUCTS and GLOBALS and THE LIKE */
enum Consistency { consistent, inconsistent };

#define MAX_PEOPLE 15 // The max. number of people in an ntuple
#define MAX_COMBINATIONS 110 // The max. number of allele combinations
#define MAX_ALLELES 16  // The max. number of alleles in the polymorphism

struct NTuple {
       char people[MAX_PEOPLE]; // Stores allele combiation code
       double probability;
       int likenessCount;
       double PrN0[MAX_PEOPLE*2+1];  // Buffers Pr(N0|M,i)
};

// An AlleleCrossing stores the (coded) possible outcomes for a cross of two
// allele combinations (those that are consistent) and their probabilities.
struct HWProbability {
       int coeficient;
       int powers[MAX_ALLELES];
};
struct AlleleCrossing {
       int outcomes[MAX_COMBINATIONS+1]; // This stores 4* the probability.
       struct HWProbability probability;
};

double alpha = .5;
int n = 3;
double alleleProbabilities[MAX_ALLELES] = {0.2,0.8};

// This will store the NTuples- dynamically allocated.
struct NTuple *ntuples;
// The number of ntuples in the array that ntuple will point to.
long ntupleCount = 0;

// The number of alleles in the polymorphism.
int alleleCount = 2;
int alleleCombinationCount = 3;
char AlleleCombinations[MAX_COMBINATIONS+1][2];
struct AlleleCrossing AlleleCrossings[MAX_COMBINATIONS+1][MAX_COMBINATIONS+1];

int debug = 1;

/* PROTOTYPES */

void DefineAlleleCombinations(int debug);
double Combination(int a, int b);
double Factorial(int a);
double PrIErrors(int i, int pedSize);
enum Consistency IsConsistent(struct NTuple ntuple);
double power(double x, int y);
void DefineNTuples(int debug);
double ProbabilityOfConsistency(char *alleles, int i, int r, int s);
struct NTuple AllelesToNTuple(char *alleles, struct NTuple *ntuple);
char *NTupleToAlleles(struct NTuple ntuple, char *alleles);
void PrintAlleles(char *alleles);
void PrintNTuple(struct NTuple ntuple);
void PrintNTupleAsAlleles(struct NTuple ntuple);

int main(int argc, char *argv[]) {
	char *displayMode;
    int displayHelp = 0;

	int i, pedSize = 0;
    double sum = 0;
    double beta = 0;
    double PrN0MiPrMSum;
    double p, PrI = 0;
    double apparentErrorRate = 0, alphaIncrement = .00005;
    double pedErrorRate = 0;
    int ntuple;
    char alleles[MAX_PEOPLE*2];
    char *bp;
    time_t now;
    long statusCtr, statusMax;
    int readArg = 1;

    time(&now);

    printf("\nCalculating Undetected Errors\n");

	if (argc > 6) {
		displayMode = argv[readArg++];
		n = (int)strtol(argv[readArg++], &bp, 10);
        pedSize = n;
		if (n > MAX_PEOPLE) {
			printf("You may only specify up to %i members in the family.\n", MAX_PEOPLE);
	        displayHelp = 1;
        } else if (n < 3) {
			printf("There must be at least three members in a family (2 parents, at least 1 child).\n");
	        displayHelp = 1;
        }
        if (strcmp(argv[readArg], "true") == 0) {
        	alpha = strtod(argv[readArg+1], &bp);
            apparentErrorRate = 0.0;
            if (alpha == 0.0) {
		    	printf("You cannot have an error rate of 0.\n\n");
		        displayHelp = 1;
            }
        } else if (strcmp(argv[readArg], "observed") == 0) {
			apparentErrorRate = strtod(argv[readArg+1], &bp);
			alpha = 0.0;
            if (apparentErrorRate == 0.0) {
		    	printf("You cannot have an error rate of 0.\n\n");
		        displayHelp = 1;
            }
	    } else {
	    	printf("You did not specify a correct error type: 'true' or 'observed'.\n\n");
	        displayHelp = 1;
	    }
        readArg+=2;
		alleleCount = (int)strtol(argv[readArg++], &bp, 10);
		if (alleleCount > MAX_ALLELES-1) {
			printf("You may only specify up to %i alleles.  (You specified %i.)\n\n", MAX_ALLELES-1, alleleCount);
			displayHelp = 1;
		} else if (alleleCount <= 1) {
			printf("You did not provide a proper number of alleles.  Valid values are between 2 and %i.\n\n", MAX_ALLELES);
			displayHelp = 1;
        }
		if (argc < alleleCount + readArg - 1) {
			printf("You did not provide enough arguments.\n\n");
            displayHelp = 1;
		} else {
			sum = 0;
		    for (i=0; i<alleleCount - 1; i++) {
				alleleProbabilities[i] = atof(argv[readArg++]);
				sum += alleleProbabilities[i];
		    }
		    alleleProbabilities[alleleCount-1] = 1 - sum;
        }
	} else {
		displayHelp = 1;
	}

	if (displayHelp == 1) {
		printf("
Provide values for the following parameters on the command line (in order):
  the display mode ('web', 'terse', or 'verbose')
  n        the number of members in the family/sibs genotyped
  'true' if the next value is the true error rate or 'observed' if it
    is the observerd error rate.
  error    the error rate, true or observed based on the previous argument
  A        the number of alleles at the locus
  F1 F2... the frequency of each allele (you may omit the last)
");
		exit(1);
	}

	printf("
  number of members in the family, n = %d
  %s = %f
  the number of alleles at the locus, A = %d
  the frequency of each allele:
     ",
		n,
		(alpha!=0.0) ? "the actual rate of errors at alleles, alpha" : "the observed rate of errors of pedigrees, alpha_o",
		(alpha!=0.0) ? alpha : apparentErrorRate,
		alleleCount);
    for (i=0; i<alleleCount - 1; i++) { printf("F%d=%f, ", i+1, alleleProbabilities[i]); }
    printf("F%d=%f\n", alleleCount, alleleProbabilities[alleleCount-1]);

    // Allocate dynamic memory for the ntuples
	ntuples = (struct NTuple*)malloc(sizeof(struct NTuple) * 1);
	if (ntuples == NULL) { printf("\nFATAL: Cannot initialize enough memory for even one ntuple.\n"); exit(1); }

    // If not in web mode, display program start time.
    //if (strcmp(argv[1], "web")) printf("Program Start: %s", asctime(localtime(&now)));

    // Setup the lookup tables.
    DefineAlleleCombinations(!strcmp(argv[1],"verbose"));
    DefineNTuples(strcmp(argv[1], "web"));

    // Display and prepare status messages.
    if (strcmp(argv[1], "web")) printf("Enumerating errors and computing probabilities...  ");
    statusCtr = 0;
    statusMax = ((2*pedSize)*(ntupleCount) + 1);

    for (i=1; i<=2*pedSize; i++) {
		PrN0MiPrMSum = 0;
        if (alpha != 0.0) PrI = PrIErrors(i, pedSize);

        // If we can, skip this round if PrIErrors is negligible.
        //if (alpha!=0.0 && apparentErrorRate==0 && PrI<.00001) continue;

        // Do status/formatting text.
		if (!strcmp(argv[1], "verbose")) printf("\n");

        // Start doing magic on each ntuple.
		for (ntuple = 0; ntuple < ntupleCount; ntuple++) {

			// To begin the recursion, r (the initial recursive level)
			// must be zero.  s (the initial errorLocation) must be zero.
			NTupleToAlleles(ntuples[ntuple], alleles);
			p = ProbabilityOfConsistency(alleles, i, 0, 0);
			ntuples[ntuple].PrN0[i] = p;

			PrN0MiPrMSum += p * ntuples[ntuple].probability;

			statusCtr++;

			if (!strcmp(argv[1], "verbose")  ) {
				printf("%2ld%% ", 100*statusCtr/statusMax);
				//PrintAlleles(alleles);
				PrintNTuple(ntuples[ntuple]);
				printf(" Pr(M) * Pr(N0)=%.7f = %.7f\n", p, p * ntuples[ntuple].probability);
			} else if (strcmp(argv[1], "web") )  {
				if (statusCtr % (statusMax/5) == 0) printf("%2ld%% ", 1 + 100*statusCtr/statusMax);
			}
		}

        // If we're calculating forwards, add up beta.
		if (alpha != 0.0) {
        	beta += PrN0MiPrMSum * PrI;
        }
	}

	printf("\n");

    // Display initial analysis (going forward: beta in terms of alpha)
	if (alpha != 0) {
    	pedErrorRate = (1.0-power(1.0-alpha, 2*pedSize));
		printf("
Results:
  error detection rate, 1-beta = %.6f
  proportion of pedigrees with errors = %.6f
  observed error rate, alpha_o = %.6f
  ",
  1.0 - beta,
  pedErrorRate,
  pedErrorRate*(1.0-beta)
  		);
	}


    // Go backwards (alpha in terms of beta).
	if (apparentErrorRate != 0.0)  {
      	printf("Finding alpha...\n");
		p = 0;
		for (alpha = alphaIncrement; alpha < 1; alpha += alphaIncrement) {
			sum = 0;
			for (i=1; i<=2*pedSize; i++) {
				PrN0MiPrMSum = 0;
				for (ntuple = 0; ntuple < ntupleCount; ntuple++) {
				    PrN0MiPrMSum += ntuples[ntuple].PrN0[i] * ntuples[ntuple].probability;
				}
				sum += PrN0MiPrMSum * PrIErrors(i, pedSize);
		    }
            pedErrorRate = (1.0-power(1.0-alpha, 2*pedSize));
		    if ((pedErrorRate*(1.0-sum)) >= apparentErrorRate) {
				printf("
Results:
  true error rate, alpha = %.4f
  error detection rate, 1-beta = %.4f
  proportion of pedigrees with errors = %.4f
",
				alpha,
				1.0 - sum,
				pedErrorRate);

                break;
			}
			p = sum;
		}
	}

    // Free up the dynamic memory to store ntuples.
    free(ntuples);

    // Print citation/help information.
    printf("
When using these results, please cite the following:
	Gordon D, Heath SC, Ott J (1999) True pedigree errors more frequent
	than apparent errors for single nucleotide polymorphisms. Hum
	Hered 49:65-70
	Gordon D, Leal SM, Heath SC, Ott J (2000) An analytic solution to
	single nucleotide polymorphism error-detection rates in nuclear
	families: implications for study design. In Pacific Symposium on
	Biocomputing 5 (Altman, Dunker, Hunter, Lauderdale, Klein, eds.).
	Honolulu, Hawaii.
Questions should be directed to tauberer@for.net.
");

	// If not in web mode, display finish time.
    /*if (strcmp(argv[1], "web")) {
	    time(&now);
    	printf("\nProgram End: %s\n", asctime(localtime(&now)));
    }*/

    // Done!
    return 0;
}

int AllelesToCombination(int a, int b) {
	int c;
	for (c=1; c<alleleCombinationCount; c++) {
		if (AlleleCombinations[c][0] == MIN(a,b) && AlleleCombinations[c][1] == MAX(a,b)) return c;
    }
    printf("\nFATAL: AllelesToCombination called with invalid alleles passed (%i,%i).", a , b);
    exit(1);
    return 0;
}

// Initializes the AlleleCombinations 2-dim array.
void DefineAlleleCombinations(int debug) {
    int a,b,c;

    if (debug) printf("Generating allele combination lookup table.");

	// Find all of the combinations of the alleles.  (index begins at 1)
    alleleCombinationCount = 1;
    for (a=1; a<=alleleCount; a++) {
    for (b=a; b<=alleleCount; b++) {
		AlleleCombinations[alleleCombinationCount][0] = a;
		AlleleCombinations[alleleCombinationCount][1] = b;
        alleleCombinationCount++;
		if (alleleCombinationCount > MAX_COMBINATIONS) {
			printf("\nFATAL: This requires too many combinations of alleles.  The max. is %d.\n", MAX_COMBINATIONS);
			exit(1);
		}
	}
	}
	if (debug) printf("%i entries.\n", alleleCombinationCount);

    // For each crossing of the two alleles...
    for (a=1; a<alleleCombinationCount; a++) {
    for (b=a; b<alleleCombinationCount; b++) {
		// Reset outcomes...
		for (c = 1; c < alleleCombinationCount; c++)
	    	AlleleCrossings[a][b].outcomes[c] = 0;

		// Do a virtual punnet square marking outcomes as consistent
		// and incrementing a counter for each occurence...
		AlleleCrossings[a][b].outcomes[AllelesToCombination(AlleleCombinations[a][0], AlleleCombinations[b][0])]++;
		AlleleCrossings[a][b].outcomes[AllelesToCombination(AlleleCombinations[a][0], AlleleCombinations[b][1])]++;
		AlleleCrossings[a][b].outcomes[AllelesToCombination(AlleleCombinations[a][1], AlleleCombinations[b][0])]++;
		AlleleCrossings[a][b].outcomes[AllelesToCombination(AlleleCombinations[a][1], AlleleCombinations[b][1])]++;

		// ...but don't divide each by four (There are four outcomes in a punnet square.)
		// because we want to keep it integral.  We'll divide by four later.

		// Reset the prob. coeficient and for each combination of two
		// non-equal alleles, multiply it by two.  (This takes into account
		// the fact that a/b = b/a.)
		AlleleCrossings[a][b].probability.coeficient = 1;
		if (AlleleCombinations[a][0] != AlleleCombinations[a][1]) AlleleCrossings[a][b].probability.coeficient *= 2;
		if (AlleleCombinations[b][0] != AlleleCombinations[b][1]) AlleleCrossings[a][b].probability.coeficient *= 2;

		// Bear in mind that alleles begin with 1, but their indices begins
		// with 0.
		for (c=0; c<alleleCount; c++) {
			AlleleCrossings[a][b].probability.powers[c] = 0;
			if (AlleleCombinations[a][0] == (c+1)) AlleleCrossings[a][b].probability.powers[c]++;
			if (AlleleCombinations[a][1] == (c+1)) AlleleCrossings[a][b].probability.powers[c]++;
			if (AlleleCombinations[b][0] == (c+1)) AlleleCrossings[a][b].probability.powers[c]++;
			if (AlleleCombinations[b][1] == (c+1)) AlleleCrossings[a][b].probability.powers[c]++;
		}
	}
	}
}

/* aCb can be calculated as a(a-1)...(b+1) (a>b+1): (floats are used to avoid conversions later) */
double Combination(int a, int b) {
	double d;
	if (a < b) {
		printf("\nFATAL: Combination called a (%i) < b (%i)\n", a , b);
		exit(1);
	} else if (a == b || b == 0) {
		return 1;
	} else if ((a-1 == b) || b == 1) {
		return a;
	}
	d = Factorial(a) / (Factorial(b) * Factorial(a-b));
	return d;
}

double Factorial(int a) {
	double b = 1;
	int c;
	if (a <= 0) {
		printf("\nFATAL: Factorial called with a (%i) <= 0.\n", a);
		exit(1);
	}
	if (a > 60) printf("\nNON-FATAL: Factorial called with a (%i) > 60.  I sure wouldn't bet on results being accurate.", a);
	for (c=1; c<=a; c++) b *= c;
	return b;
}

// PrIErrors = Pr(i errors in tuple) is now:
double PrIErrors(int i, int pedSize) {
	double numerator = Combination(2*pedSize, i) * power(alpha, i) * power(1.0-alpha, (2*pedSize) - i);
    double denominator = 1.0 - power(1.0-alpha, 2*pedSize);
    if (denominator == 0) {
		printf("\nFATAL: PrIErrors(n=%i,alpha=%g,i=%i) called returning undefined as denominator is zero.\n", n, alpha, i);
		exit(1);
	}
    return (numerator/denominator);
}

double power(double x, int y) {
	int i;
	double result = 1.0;
	if (y < 0) {
		printf("\nFATAL: pow2 called with y < 0.\n");
		exit(1);
	}
	for (i=0; i<y; i++) result *= x;
	return result;
}

void DefineNTuples(int debug) {
	int counting[MAX_PEOPLE+1];
	int a,b,ct,f;
	struct AlleleCrossing c;
	int needMoreSpace = 1;
	int combinationCount[MAX_COMBINATIONS];
	int combinationCount2[MAX_COMBINATIONS];
    unsigned long statusCounter = 0, statusMax;

    if (debug) printf("Generating list of ntuples...  ");

    for (a=0; a<=n; a++) counting[a] = 1;

    statusMax = power(alleleCombinationCount-1, n);

	while (counting[0] == 1) {
		if (debug) {
	        statusCounter++;
	        if (statusCounter % (statusMax/5) == 0) { printf("%li%% ", (long int)((double)statusCounter/(double)statusMax*100.0)); }
        }

		// Here we have a unique array- turn it into a ntuple
		for (a=n; a>0; a--)
			ntuples[ntupleCount].people[a-1] = counting[a];

		ntuples[ntupleCount].likenessCount = 1;

		// If this is a consistent ntuple...
		if (IsConsistent(ntuples[ntupleCount]) == consistent) {
			// See if there's a previos ntuple that is this one in a
			// different order.
            f = 0;
			if (ntupleCount < 7500) {
				for (ct=1; ct<alleleCombinationCount; ct++)
					combinationCount[ct] = 0;
				for (ct=2; ct<n; ct++)
					combinationCount[(int)ntuples[ntupleCount].people[ct]]++;
				for (b=0; b<ntupleCount; b++) {
					if (MIN(ntuples[b].people[0],ntuples[b].people[1])!=MIN(ntuples[ntupleCount].people[0],ntuples[ntupleCount].people[1]) || MAX(ntuples[b].people[0],ntuples[b].people[1])!=MAX(ntuples[ntupleCount].people[0],ntuples[ntupleCount].people[1])) continue;
					for (ct=1; ct<alleleCombinationCount; ct++)
						combinationCount2[ct] = 0;
					for (ct=2; ct<n; ct++)
						combinationCount2[(int)ntuples[b].people[ct]]++;
					f = 0;
					for (ct=1; ct<alleleCombinationCount; ct++)
						if (combinationCount[ct] != combinationCount2[ct]) { f = 1; break; }
					if (f == 0) {
						ntuples[b].likenessCount++;
						f = 1;
						break;
					} else f = 0;
				}
			} else {
				f = 0;
			}

			if (f == 0) {
				// compute the probability of the ntuple
				c = AlleleCrossings[MIN(ntuples[ntupleCount].people[0], ntuples[ntupleCount].people[1])][MAX(ntuples[ntupleCount].people[0], ntuples[ntupleCount].people[1])];
				ntuples[ntupleCount].probability = (double) c.probability.coeficient;
				for (a=0; a<alleleCount; a++) { // p^n * q^m ...
					if (c.probability.powers[a] != 0) ntuples[ntupleCount].probability *= power(alleleProbabilities[a], c.probability.powers[a]);
				}
				for (a=2; a<n; a++) { // * prob for each child (divide by four because we didn't initially)
					ntuples[ntupleCount].probability *= 0.25 * c.outcomes[(int)ntuples[ntupleCount].people[a]];
	            }
    			// and increment the counter so it's saved in the array.
				ntupleCount++;

                // This dynamically allocates array space, and at this point
                // it needs more.
				if (ntupleCount == needMoreSpace) {
					if ((ntuples=realloc(ntuples, 3*sizeof(struct NTuple)*(ntupleCount+1))) == NULL) {
						printf("\nFATAL: Cannot increase NTuple allocation.  Stopped at %li.\n", ntupleCount-1);
					    exit(1);
					}
					needMoreSpace = ntupleCount*3;
				}
			}
		}

		// add one to the "units" digit
		counting[n]++;

		// Increment digits
		for (a=n; a>0; a--) {
			if (counting[a] == (alleleCombinationCount)) { // should never be greater than this
				counting[a] = 1;
				counting[a-1]++;
			}
		}
	}

	// Multiply the likenessCount by the calculated probs.
    for (ct=0; ct<ntupleCount; ct++)
		ntuples[ct].probability *= ntuples[ct].likenessCount;

    if (debug) printf("\n");
}

// Tests if ntuple is consistent by comparing children to parents.
enum Consistency IsConsistent(struct NTuple ntuple) {
	int a,x,y;
    x = MIN(ntuple.people[0], ntuple.people[1]);
    y = MAX(ntuple.people[0], ntuple.people[1]);
    for (a=2; a<n; a++) {
		if (x < 1 || y < 1 || ntuple.people[a] < 1) {
	    	printf("\nERROR: strange thing found in checking consistency.\n");
            PrintNTuple(ntuple);
	    	exit(1);
		}
		if (AlleleCrossings[x][y].outcomes[(int)ntuple.people[a]] == 0.0) return inconsistent;
	}
    return consistent;
}


// Returns the probability that ntuple with i errors introduced
// will be consistent.  This function is recursive, crazily too.
// To simplify matters, the "return globals" below are used
// to retain information during the recurision.
long PermutationCount = 0;
long ConsistentCount = 0;
double ProbabilityOfConsistency(char *alleles, int i, int r, int s) {
	char newAlleles[n*2];
	int errorLocation;
	struct NTuple ntuple;
	int errorUpDown;

	if (r==0) {
		PermutationCount = 0;
		ConsistentCount = 0;
	}

	for (errorLocation = s; errorLocation<n*2 - i + 1 + r; errorLocation++) {
	for (errorUpDown = 0; errorUpDown < 2; errorUpDown++) {
		// set the working copy to the passed array
		memcpy(newAlleles, alleles, n*2);

		// Make the error at the location.
        if (errorUpDown == 0) {
        	if (newAlleles[errorLocation] == 1) {
        		newAlleles[errorLocation]++;
            } else if (newAlleles[errorLocation] == alleleCount) {
            	newAlleles[errorLocation]--;
            } else {
            	newAlleles[errorLocation]++;
            }
        } else {
         	if (newAlleles[errorLocation] == 1) {
				newAlleles[errorLocation]++;
			} else if (newAlleles[errorLocation] == alleleCount) {
				newAlleles[errorLocation]--;
            } else {
				newAlleles[errorLocation]--;
            }
        }

        // At this stage in the recursion all of the errors have been
        // introduced, so we can add to the running counts of information.
		if ((i == r+1)) {
	   		AllelesToNTuple(newAlleles, &ntuple);
			PermutationCount++;
			if (IsConsistent(ntuple) == consistent) ConsistentCount++;

        // But in this stage, we haven't introduced all of the errors,
        // so we have to recurse another level.
		} else {
           ProbabilityOfConsistency(newAlleles, i, r+1, errorLocation+1);
        }
	}
    }

    if (r == 0) {
		if (PermutationCount == 0) {
	    	printf("\nFATAL: ProbabilityOfConsistency called (i=%d) returning undefined as PermutationCount equals 0.\n", i);
	    	exit(1);
		} else {
			if (!debug && ConsistentCount < PermutationCount) {
				printf("i=%i <", i);
				//PrintAlleles(alleles);
                PrintNTuple(AllelesToNTuple(alleles, &ntuple));
				printf(">");
				printf(" %li/%li\n", ConsistentCount, PermutationCount);
            }
			return ((double)ConsistentCount)/((double)PermutationCount);
		}
	} else {
    	return 0;
	}
}

struct NTuple AllelesToNTuple(char *alleles, struct NTuple *ntuple) {
	int p, a, f;
	for (p=0; p<n; p++) {
		f = 0;
		for (a=1; a<alleleCombinationCount; a++) {
			if (MIN(alleles[2*p],alleles[2*p +1])==AlleleCombinations[a][0] && MAX(alleles[2*p],alleles[2*p +1])==AlleleCombinations[a][1]) {
				f = 1;
				ntuple->people[p] = a;
				break;
			}
		}
		if (f == 0) {
			printf("\nFATAL: AllelesToNTuple found improper allele combination: p=%i=>(%i, %i).", p, alleles[2*p], alleles[2*p +1]);
			exit(1);
		}
	}
	return *ntuple;
}

char *NTupleToAlleles(struct NTuple ntuple, char *alleles) {
	int p;
	for (p=0; p<n; p++) {
		alleles[p*2] = AlleleCombinations[(int)ntuple.people[p]][0];
		alleles[p*2 +1] = AlleleCombinations[(int)ntuple.people[p]][1];
    }
    return alleles;
}

void PrintAlleles(char *alleles) {
	int a;
    for (a=0; a<n-1; a++)
		printf("%i/%i,", alleles[2*a], alleles[2*a +1]);
	printf("%i/%i", alleles[2*(n-1)], alleles[2*(n-1) +1]);
}

void PrintNTuple(struct NTuple ntuple) {
	int a;
	for (a=0; a<n-1; a++)
		printf("%i,", ntuple.people[a]);
	printf("%i", ntuple.people[n-1]);
}

void PrintNTupleAsAlleles(struct NTuple ntuple) {
	char alleles[n*2];
	PrintAlleles(NTupleToAlleles(ntuple, alleles));
}