/*****************************************************************************/
/* one-two.c								     */
/*									     */
/* two-locus-penetrance matrix derived from single-locus penetrance vectors  */
/*									     */
/* Manuel Mattheisen							     */
/* Institute for Medical Biometry, 					     */
/* Informatics, and Epidemiology, 				             */
/* University of Bonn                                                        */
/* 									     */
/* mail to: mmattheisen@uni-bonn.de					     */
/* or phone: 0049-228-287 5564						     */
/*									     */
/* Copyright (c) Manuel Mattheisen 02/2004, revised 1/2006, 		     */
/*	according to Strauch et al., Hum Hered. 2003;56(4):200-11            */
/*									     */
/* Number of formulae mentioned in this program refer to this publication    */
/*                                                                           */
/* Input:  two penetrance vectors with disease allele frequencies            */
/*                                                                           */
/* Output: file with two-locus matrix and disease allele                     */
/*         frequencies                                                       */
/*                                                                           */
/* NOTE: GENEHUNTER-TWOLOCUS is a program that takes imprinting              */
/*	 into account. 							     */
/*								 	     */
/*	 Penetrance vectors, manually entered, have to contain 		     */
/*	 four penetrances! 						     */
/*									     */
/*****************************************************************************/

/*
                            NO WARRANTY

THIS SOFTWARE IS PROVIDED 'AS IS' AND ANY EXPRESS OR IMPLIED WARRANTIES,
INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  THE ENTIRE RISK AS TO
THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.  SHOULD THE PROGRAM
PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, REPAIR, OR
CORRECTION.  IN NO EVENT SHALL THE COPYRIGHT HOLDER BE LIABLE FOR ANY DIRECT,
INDIRECT, INCIDENTAL, GENERAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; DATA BEING RENDERED INACCURATE; FAILURE OF
THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS; BUSINESS INTERRUPTION; OR
OTHER LOSSES) SUSTAINED BY YOU OR THIRD PARTIES HOWEVER CAUSED AND ON ANY
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OR
INABILITY TO USE THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
DAMAGES.

*/




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

#define MAX_NAME_SIZE 100	/* maximum length of pathname for inputfiles */
#define LINES 4  
#define COLUMNS 4

int main(void)          /* main program */
	{	
	FILE *fp_in1, *fp_in2, *fp_out;
	char intro[4], in_outf[4], in_outf2[3], input[200], input1[200], input2[200], model2[3], again[4];
	int model;
	int flag, flag1, flag2, flag3;
	int adjust = 0;
	int i, j;
	char outfile[MAX_NAME_SIZE];	
		/* name of outfile */
	double pen_vec1[4];	
		/* penetrance vector for disease locus 1 (dl1) */ 
	double maf1, maf2;	
		/* disease allele frequencies for dl1 and dl2 */
	double wta1, wta2; 		
		/* wild type allele frequences for dl1 and dl2 */
	double paf_1, paf_2;
		/* population attributable fractions */
	double pen_vec2[4]; 		
		/* penetrance vector for disease locus 2 (dl2) */ 
	double pen_matrix[LINES][COLUMNS];	
		/* penetrance matrix for two-locus */
	double k_1, k_2, k_m;		
		/* prevalences  */
	double matrix_max;		
		/* highest value in two-locus matrix */
	double pen_vec_c1[4];		
		/* control of penetrance vector 1 */
	double pen_vec_c2[4];		
		/* control of penetrance vector 2 */

	system("/usr/bin/clear");
	printf("\n\n'One-Two' is a program which enables you to easily derive a"); 
	printf("\ntwo-locus model from two single-locus models."); 

	printf("\n\nAs a 'first-time-user' be sure to read the introduction!");
	printf("\nDo you want to read the short introduction?\n");

	do {
		flag = 0;
		printf("Type (y)es or (n)o: ");
		fgets(input, sizeof(input), stdin);
		sscanf(input, "%3s", intro);
		if ((strcmp(intro, "yes") == 0) || (strcmp(intro,"YES") == 0) || 
			(strcmp(intro, "y") == 0) || (strcmp(intro,"Y") == 0)) {
			flag = 1;
	                system("/usr/bin/clear");
			printf("This program enables you to easily derive a two-locus model"); 
			printf("\nfrom two single-locus models. Therefore the single-locus"); 
			printf("\nmodels (penetrance vectors and disease allele frequencies)");
			printf("\nhave to be entered manually.");
			printf("\n\nPenetrance vectors have to contain four penetrances, i.e.,");
			printf("\nthe imprinting formulation of disease models is assumed here.");
			printf("\nPenetrances should be given in the order f+/+ fm/+ f+/m fm/m,"); 
			printf("\nseparated by spaces, with '+' denoting the wild-type and 'm'");
			printf("\ndenoting the disease allele, and the paternally inherited");
			printf("\nallele listed first. In case of no imprinting, please list");
			printf("\nthe heterozygote penetrance twice, i.e., the penetrance vector");
			printf("\nshould be f+/+ fHet fHet fm/m."),
			printf("\n\nAfter choosing the formula you want to apply (for further");
			printf("\ninformation: Strauch et al., Hum Hered. 2003;56(4):200-11),");
			printf("\na screen output of the derived two-locus penetrance matrix is given."); 
			printf("\nAn adaption of the penetrance matrix is necessary in two cases.");
			printf("\nIf the highest two-locus penetrance in this matrix is > 1,");
			printf("\nall penetrances will be divided by that value so that all");
			printf("\nelements of the matrix are <= 1. If values in the matrix");
			printf("\nare < 0, these will be set to 0.");
			printf("\n\n              (Hit return to continue)\n");
			fgets(input, sizeof(input), stdin);
			system("/usr/bin/clear");

			printf("Applying an averaged-model formula in the beginning will");
			printf("\nlead to an on-screen comparison of the single-locus models");
			printf("\nyou entered in contrast to the averaged single-locus models");
			printf("\nderived from the two-locus model. A deviation will occur");
			printf("\nif either the single-locus prevalences for both single-locus");
			printf("\nmodels differ, or if an adaption of the two-locus penetrance");
			printf("\nis necessary.");
			printf("\n\nIn the end you will be requested to enter the name of your");
			printf("\npenetrance file which will be used by GENEHUNTER-TWOLOCUS.");
			printf("\nDisease allele frequencies and the two-locus penetrance");
			printf("\nmatrix will be stored to this file.");
			printf("\n\n\nIf you have further questions, please feel free to contact me:");
			printf("\nManuel Mattheisen, IMBIE Bonn, mmattheisen@uni-bonn.de");
			printf("\n\n              (Hit return to continue)\n");
			fgets(input, sizeof(input), stdin);
			system("/usr/bin/clear");
		} else if ((strcmp(intro, "no") == 0) || (strcmp(intro,"NO") == 0) || 
			(strcmp(intro, "n") == 0) || (strcmp(intro,"N") == 0)) {
			flag = 1;
			system("/usr/bin/clear");
		} 
	} while (flag !=1); 


	printf("\nNOTE: Penetrance vectors have to contain four penetrances!\n");
	
	printf("\nPlease enter single-locus model for first disease locus. ");
	printf("\ndisease allele frequency: ");
	fgets(input, sizeof(input), stdin);
	sscanf(input, "%lf", &maf1);
	do {
		pen_vec1[3] = (-999.9999);
		printf("penetrance vector: ");
		fgets(input, sizeof(input), stdin);
		sscanf(input, "%lf %lf %lf %lf", 
			&pen_vec1[0], &pen_vec1[1], &pen_vec1[2], &pen_vec1[3]);
	} while (pen_vec1[3] == -999.9999);			
	
	printf("\nPlease enter single-locus model for the second disease locus.");
	printf("\ndisease allele frequency: ");
	fgets(input, sizeof(input), stdin);
	sscanf(input, "%lf", &maf2);
	do {
		pen_vec2[3] = (-999.9999);
		printf("penetrance vector: ");
		fgets(input, sizeof(input), stdin);
		sscanf(input, "%lf %lf %lf %lf", 
			&pen_vec2[0], &pen_vec2[1], &pen_vec2[2], &pen_vec2[3]);
	} while (pen_vec2[3] == -999.9999);
	 
	
	printf("\npenetrance vector 1 is           : %lf %lf %lf %lf ", 
		pen_vec1[0], pen_vec1[1], pen_vec1[2], pen_vec1[3]); 
	printf("\ndisease allele frequency 1 is    : %lf\n", maf1);
	wta1 = (1-maf1);
	printf("wild type allele frequency 1 is  : %lf\n", wta1);
	k_1 =((pen_vec1[0]*wta1*wta1)+((pen_vec1[1]+pen_vec1[2])*maf1*wta1)+(pen_vec1[3]*maf1*maf1));
	printf("prevalence k(1)                  : %14.12lf\n", k_1);          
	paf_1=((k_1-pen_vec1[0])/k_1);
	printf("population attrib. fraction (1)  : %lf\n", paf_1);          
	printf("\npenetrance vector 2 is           : %lf %lf %lf %lf ", 
		pen_vec2[0], pen_vec2[1], pen_vec2[2], pen_vec2[3]);
	printf("\ndisease allele frequency 2 is    : %lf\n", maf2);
	wta2 = (1-maf2);
	printf("wild type allele frequency 2 is  : %lf\n", wta2);
	k_2 =((pen_vec2[0]*wta2*wta2)+((pen_vec2[1]+pen_vec2[2])*maf2*wta2)+(pen_vec2[3]*maf2*maf2));
	printf("prevalence k(2)                  : %14.12lf\n", k_2);	 
	paf_2=((k_2-pen_vec2[0])/k_2);
	printf("population attrib. fraction (2)  : %lf\n", paf_2);          
	k_m =((k_1 + k_2)/2);
	if (k_1 != k_2) {
		printf("\n\nNOTE: The prevalences K(1) and K(2) are different!");
		printf("\n      In the following the average value K(a) will"); 
		printf("\n      be used in the formulas for the averaged models.");
	}	
	do {
		do {
			printf("\n\nCreate a (h)eterogeneity, (m)ultiplicative,"); 
			printf("\nor an (a)dditive two-locus model?: ");
			model = getchar();
			getchar();
			if (model == 'h') {					/* heterogeneity models */
                                printf("\nNote: If you have obtained the single-locus penetrances by");
                                printf("\n      a MOD-score analysis, such as with GENEHUNTER-MODSCORE,");
				printf("\n      the averaged-model formula makes more sense than the");
                                printf("\n      component-mechanism formula!");
				printf("\n\nDo you wish to apply the component-mechanism formula (2)");
	         		printf("\nor the formula (9) which derives the two-locus model from the");
				printf("\naveraged single-locus models?");
				do {
					printf("\n\nPlease enter the number of the formula you want to apply: ");
					fgets (input, sizeof(input), stdin);
					sscanf(input, "%2s", model2);
					if (strcmp(model2, "9") == 0 ) {
						for (i=0 ; i<LINES ; i++) 
							for (j=0 ; j<COLUMNS ; j++) 
								pen_matrix[i][j] = (1-((1-pen_vec1[i])*(1-pen_vec2[j]))/(1-k_m));
					} else if (strcmp(model2, "2") == 0) {
						for (i=0 ; i<LINES ; i++) 
							for (j=0 ; j<COLUMNS ; j++) 
								pen_matrix[i][j] =
									 ((pen_vec1[i]+pen_vec2[j])-(pen_vec1[i]*pen_vec2[j]));
					} 
				} while ((strcmp(model2, "2") != 0) && (strcmp(model2, "9") != 0));
			} else if (model == 'm') {		 				/* multiplicative models */
                                printf("\nNote: If you have obtained the single-locus penetrances by");
                                printf("\n      a MOD-score analysis, such as with GENEHUNTER-MODSCORE,");
				printf("\n      the averaged-model formula makes more sense than the");
                                printf("\n      component-mechanism formula!");
				printf("\n\nDo you wish to apply the component-mechanism formula (11)");
	         		printf("\nor the formula (14) deriving the two-locus model from the");
				printf("\naveraged single-locus models?");
				do {
					printf("\n\nPlease enter the number of the formula you want to apply: ");
					fgets (input, sizeof(input), stdin);
					sscanf(input, "%2s", model2);
					if (strcmp(model2, "14") == 0 ) {
						for (i=0 ; i<LINES ; i++) 
							for (j=0 ; j<COLUMNS ; j++) 
								pen_matrix[i][j] = (pen_vec1[i]*pen_vec2[j]/k_m);
					} else if (strcmp(model2, "11") == 0) {
						for (i=0 ; i<LINES ; i++) 
							for (j=0 ; j<COLUMNS ; j++) 
								pen_matrix[i][j] = (pen_vec1[i]*pen_vec2[j]);
					} 
				} while ((strcmp(model2, "14") != 0) && (strcmp(model2, "11") != 0));
			} else if (model == 'a') {					/* additive models */
                                printf("\nNote: If you have obtained the single-locus penetrances by");
                                printf("\n      a MOD-score analysis, such as with GENEHUNTER-MODSCORE,");
				printf("\n      the averaged-model formula makes more sense than the");
                                printf("\n      component-mechanism formula!");
				printf("\n\nDo you wish to apply the component-mechanism formula (16)");
	         		printf("\nor the formula (19), deriving the two-locus model from the");
				printf("\naveraged single-locus models?");
				do {
					printf("\n\nPlease enter the number of the formula you want to apply: ");
					fgets (input, sizeof(input), stdin);
					sscanf(input, "%2s", model2);
					if (strcmp(model2, "19") == 0 ) {
						for (i=0 ; i<LINES ; i++) 
							for (j=0 ; j<COLUMNS ; j++) 
								pen_matrix[i][j] = (pen_vec1[i]+pen_vec2[j]-k_m);
					} else if (strcmp(model2, "16") == 0) {
						for (i=0 ; i<LINES ; i++) 
							for (j=0 ; j<COLUMNS ; j++) 
								pen_matrix[i][j] = (pen_vec1[i]+pen_vec2[j]);
					} 
				} while ((strcmp(model2, "16") != 0) && (strcmp(model2, "19") != 0));
			}
		} while ((model != 'm') && (model != 'h') && (model != 'a'));
		printf("\nYour two-locus matrix:\n");                /* screen ouput of matrix, not adjusted */
		printf("  %14.12lf %14.12lf %14.12lf %14.12lf\n", 
			pen_matrix[0][0], pen_matrix[0][1], pen_matrix[0][2], pen_matrix[0][3]);	
		printf("  %14.12lf %14.12lf %14.12lf %14.12lf\n", 
			pen_matrix[1][0], pen_matrix[1][1], pen_matrix[1][2], pen_matrix[1][3]);	
		printf("  %14.12lf %14.12lf %14.12lf %14.12lf\n", 
			pen_matrix[2][0], pen_matrix[2][1], pen_matrix[2][2], pen_matrix[2][3]);	
		printf("  %14.12lf %14.12lf %14.12lf %14.12lf\n", 
			pen_matrix[3][0], pen_matrix[3][1], pen_matrix[3][2], pen_matrix[3][3]);	
		matrix_max = pen_matrix[0][0]; 
		for (i=0 ; i<LINES ; i++)				/* Search for highest value in matrix */ 
			for (j=0 ; j<COLUMNS ; j++) { 
				if (matrix_max < pen_matrix[i][j]) 
					matrix_max = pen_matrix[i][j];
			}		
		for (i=0 ; i<LINES ; i++) 				/* Adjusting produced matrix */
			for (j=0 ; j<COLUMNS ; j++) { 
				if (pen_matrix[i][j] < 0) { 		/* values less than 0 are fixed 0 */
					pen_matrix[i][j] = 0;
					adjust = 1;			
				} else if (matrix_max > 1) {	
									/* highest value in matrix > 1: */ 
					pen_matrix[i][j] = (pen_matrix[i][j]/matrix_max);	
									/* matrix scaled by this value */
					adjust = 1;
				}
			}
		if ((model == 'a') && (adjust == 1)) {
			printf("\nWARNING: You choose an additive component-mechanism formula and");
			printf("\n         the highest penetrance in your two-locus matrix is >1!");
		}
		if (adjust == 1) {             /* again screen-output if matrix adjusted */
			printf("\n\nNOTE: Your matrix was adjusted:\n");
			printf("  %14.12lf %14.12lf %14.12lf %14.12lf\n", 
				pen_matrix[0][0], pen_matrix[0][1], pen_matrix[0][2], pen_matrix[0][3]);	
			printf("  %14.12lf %14.12lf %14.12lf %14.12lf\n", 
				pen_matrix[1][0], pen_matrix[1][1], pen_matrix[1][2], pen_matrix[1][3]);	
			printf("  %14.12lf %14.12lf %14.12lf %14.12lf\n", 
				pen_matrix[2][0], pen_matrix[2][1], pen_matrix[2][2], pen_matrix[2][3]);	
			printf("  %14.12lf %14.12lf %14.12lf %14.12lf\n", 
				pen_matrix[3][0], pen_matrix[3][1], pen_matrix[3][2], pen_matrix[3][3]);	
		}
		for (i=0 ; i<4 ; i++)        				  	/* Controls for pen_vec1[i] */ 
			pen_vec_c1[i] = ((pen_matrix[i][0]*wta1*wta1)+
					 (pen_matrix[i][1]*wta1*maf1)+
					 (pen_matrix[i][2]*wta1*maf1)+
					 (pen_matrix[i][3]*maf1*maf1));
		for (i=0 ; i<4 ; i++) 						/* Controls for pen_vec2[j] */
			pen_vec_c2[i] = ((pen_matrix[0][i]*wta2*wta2)+
					 (pen_matrix[1][i]*wta2*maf2)+
					 (pen_matrix[2][i]*wta2*maf2)+
					 (pen_matrix[3][i]*maf2*maf2));
		printf("\n\nEntered penetrances for dis.loc1: %lf %lf %lf %lf",
				pen_vec1[0], pen_vec1[1], pen_vec1[2], pen_vec1[3]); 
		printf("\nInferred penetrances dis.loc1   : %lf %lf %lf %lf", 
				pen_vec_c1[0], pen_vec_c1[1], pen_vec_c1[2], pen_vec_c1[3]);
		printf("\n\nEntered penetrances for dis.loc2: %lf %lf %lf %lf",
				pen_vec2[0], pen_vec2[1], pen_vec2[2], pen_vec2[3]); 
		printf("\nInferred penetrances dis.loc2   : %lf %lf %lf %lf", 
				pen_vec_c2[0], pen_vec_c2[1], pen_vec_c2[2], pen_vec_c2[3]);
		printf("\n\n\nPlease enter name of outfile: ");
		do { 
			scanf("%s", outfile);
			getchar();
			strcpy(in_outf2, "yes");
	 		if ((fp_out = fopen (outfile, "r")) != NULL) {
				strcpy(in_outf2, "no");
				printf("\nfile %s already exists, overwrite?\n", outfile);
				printf("Type 'yes' or 'no': ");
				scanf("%s", in_outf2);
				getchar();
				if ((strcmp(in_outf2,"no") == 0) || (strcmp(in_outf2,"n") == 0)) {
					printf("\nPlease enter name of outfile: ");
				} else if ((strcmp(in_outf2, "yes") == 0) || (strcmp(in_outf2, "y") == 0)) {
					strcpy(in_outf2, "yes");
				}
			}
		} while ((strcmp(in_outf2, "yes") != 0));
		fp_out = fopen (outfile, "w");                              /* file-output */
		fprintf(fp_out, "%lf %lf           << DISEASE ALLELE FREQUENCY DISEASE LOCUS 1\n", wta1, maf1);
		fprintf(fp_out, "%lf %lf           << DISEASE ALLELE FREQUENCY DISEASE LOCUS 2\n", wta2, maf2);
		fprintf(fp_out, "1                             << NUMBER OF LIABILITY CLASSES\n");
		fprintf(fp_out, "%14.12lf %14.12lf %14.12lf %14.12lf \n", 
			pen_matrix[0][0], pen_matrix[0][1], pen_matrix[0][2], pen_matrix[0][3]);
		fprintf(fp_out, "%14.12lf %14.12lf %14.12lf %14.12lf \n", 
			pen_matrix[1][0], pen_matrix[1][1], pen_matrix[1][2], pen_matrix[1][3]);
		fprintf(fp_out, "%14.12lf %14.12lf %14.12lf %14.12lf \n", 
			pen_matrix[2][0], pen_matrix[2][1], pen_matrix[2][2], pen_matrix[2][3]);
		fprintf(fp_out, "%14.12lf %14.12lf %14.12lf %14.12lf \n", 
			pen_matrix[3][0], pen_matrix[3][1], pen_matrix[3][2], pen_matrix[3][3]);
		fclose(fp_out);

		printf("\n\n\nDo you wish to derive another two-locus matrix on the");	
		printf("\nbasis of the same entered single-locus penetrances?");
		printf("\nType (y)es or (n)o: ");
		fgets(input, sizeof(input), stdin);
		sscanf(input, "%3s", again);
		if ((strcmp(again, "yes") == 0) || (strcmp(again,"YES") == 0) || 
			(strcmp(again, "y") == 0) || (strcmp(again,"Y") == 0)) {
			flag3 = 1;
			system("/usr/bin/clear");
			printf("\npenetrance vector 1 is           : %lf %lf %lf %lf ", 
				pen_vec1[0], pen_vec1[1], pen_vec1[2], pen_vec1[3]); 
			printf("\ndisease allele frequency 1 is    : %lf\n", maf1);
			printf("wild type allele frequency 1 is  : %lf\n", wta1);
			printf("prevalence k(1)                  : %14.12lf\n", k_1);          
			printf("\npenetrance vector 2 is           : %lf %lf %lf %lf ", 
				pen_vec2[0], pen_vec2[1], pen_vec2[2], pen_vec2[3]);
			printf("\ndisease allele frequency 2 is    : %lf\n", maf2);
			printf("wild type allele frequency 2 is  : %lf\n", wta2);
			printf("prevalence k(2)                  : %14.12lf\n", k_2);	 
			if (k_1 != k_2) {
				printf("\n\nNOTE: The prevalences K(1) and K(2) are different!");
				printf("\n      In the following the average value K(a) will"); 
				printf("\n      be used in the formulas for the averaged models.");
			}	
		} else if ((strcmp(again, "no") == 0) || (strcmp(again,"NO") == 0) || 
			(strcmp(again, "n") == 0) || (strcmp(again,"N") == 0)) {
			flag3 = 0;
		}	
	} while (flag3 != 0);
	
	printf("\n\n\n  bye...\n\n\n");	
	return(0);
}
