#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<ctype.h>
#include "subfunction1.c"
int RX,RY,RZ;
int nsample, nmarker, **gen;
 /*
   nsample: number of sampled individuals
   marker:  number of markers
   gen   :  matrix (i,j)th element is the genotype of 
            ith individual at jth marker
  */
  int compHG(int *hap, int *sgen, int nmarker, int *hapcom);
  /* test if haplotype hap is compatible with genotype sgen. 
     if hap is compatible with sgen, than find the haplotype hapcom 
	 such that hap and hapcom together is sgen
  */
   void compatible( int *sgen, int nmarker, int *numall);
    /* find all haplotypes pairs compatible with genotype sgen
       and write all these haplotype pairs in a file
    */

main()
{
	int i,j,nn,*numall;
	FILE *fp;
	nsample=200;
	nmarker=5;
	gen=imatrix(1,nsample,1,nmarker);
	numall=ivector(1,nmarker);
	for(i=1;i<=nmarker;i++)numall[i]=2;
    /* read genotype data from file hapdata.txt to matrix gen */

	fp=fopen("hapdata.txt","r");
	if(fp==NULL)
	{
		printf("can't open file hapdata.txt for reading\n");
		exit(1);
	}
	
	for(i=1;i<=nsample;i++)
	{
		fscanf(fp,"%d",&nn);
		for(j=1;j<=nmarker;j++)fscanf(fp,"%d",&gen[i][j]);
	}
    compatible(gen[35],nmarker,numall);
	free_imatrix(gen,1,nsample,1,nmarker);
	free_ivector(numall,1,nmarker);

}

//=========================================================  
  int compHG(int *hap, int *sgen, int nmarker, int *hapcom)
  {
	  int i, indicator;
	  indicator=1;
	  for(i=1;i<=nmarker;i++)
	  {
		  if(((hap[i]==0)&&(sgen[i]==2))||((hap[i]==1)&&(sgen[i]==0)))
		  {
			  indicator=0;
			  break;
		  }
	  }
	  if(indicator) //same as if(indicator==1)
	  {
		  for(i=1;i<=nmarker;i++)hapcom[i]=sgen[i]-hap[i];
	  }
	  return indicator;
  }

//===========================================================
  void compatible( int *sgen, int nmarker, int *numall)
  {
	  int i,j,nn,sum, *hap, *hapcom,numhap;
	  FILE *fp;

	  hap=ivector(1,nmarker);
	  hapcom=ivector(1,nmarker);
	  numhap=ipower(2,nmarker);
	  if((fp=fopen("haplotype.txt","w"))==NULL)
	  {
		  printf("can't open file haplotype.txt for writing\n");
		  exit(1);
	  }
	  fprintf(fp,"Genotype\n");
	  for(i=1;i<=nmarker;i++)
	  {
		  if(sgen[i]==0)fprintf(fp,"0/0 ");
		  else if(sgen[i]==1)fprintf(fp,"0/1 ");
		  else fprintf(fp,"1/1 ");
	  }
	  fprintf(fp,"\n\n");
	  sum=0;
	  for(i=1;i<=numhap;i++)
	  {
		  Arrayindex(i-1,hap,nmarker,numall);
		  if(compHG(hap,sgen,nmarker,hapcom))
		  {
			  nn=INdex(hapcom,nmarker,numall);
			  if(i-1<=nn)
			  {
			    sum++;
			    fprintf(fp,"\nHaplotype Pair %d\n",sum);
				for(j=1;j<=nmarker;j++)fprintf(fp,"%d ",hap[j]);
				fprintf(fp,"\n");
				for(j=1;j<=nmarker;j++)fprintf(fp,"%d ",hapcom[j]);
				fprintf(fp,"\n");
			  }
		  }
	  }
	  fclose(fp);
	  free_ivector(hap,1,nmarker);
	  free_ivector(hapcom,1,nmarker);
  }