#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#define F 0
#define T 1
#define MAXLOCI 12
#define MAXALL 15
#define MAXCOM (MAXALL*(MAXALL+1))/2
#define MAXFAC 10
#define MAXLINE 1000


  typedef struct locusvalues {
    int nalleles, format;
    double freq[MAXALL];
    union {
      struct {
	int nclass;
	double pene[MAXCOM][MAXCOM];
      } aff;
      struct {
	int numfact;
	int binfactors[MAXALL][MAXFAC];
      } binary;
    } type;
  } locusvalues;

locusvalues locuslist[MAXLOCI];

int table[MAXLOCI][MAXALL][MAXALL];


int getline(s,lim,file)
     char s[];
     int lim;
     FILE *file;
{
  int c,i;
  for (i=0; i< lim-1 && (c=fgetc(file))!=EOF && c!= '\n'; ++i)
    s[i] = c;
  
  if (c== '\n')
    {
      s[i]=c;
      ++i;
    }
  s[i] = '\0';
  return i;
}

void putline(s,lim,file)
     char s[];
     int lim;
     FILE *file;
{
  int c,i;
  
  for (i=0; i<lim; i++)
    {
      c = s[i];
      fputc(c,file);
    }
}      


  
void binand(a,b,c,length)
     int a[MAXFAC], b[MAXFAC], c[MAXFAC], length;
{
  int i;

  for (i=0; i<length; i++)
    {
      if ((a[i]==b[i])&&(a[i]==1))
	c[i] = 1;
      else 
	c[i] = 0;
    }
}

void binaff (bin,aff,liab,locus,length)
     int bin[],*aff,*liab,locus,length;
{
  int i,j, where,allone,k,result[MAXFAC],a[2],count,inversehap[MAXALL][MAXALL];
  
  count = 0;
  for (i=0; i<length; i++)
    if (bin[i]==1)
      count += 1;
    
  if (count==0)
    {
      *aff = 0;
      *liab = 1;
    }
  else
    {
      where = 0;
      for (j=0; j<locuslist[locus].nalleles; j++)
	{
	  binand(bin,locuslist[locus].type.binary.binfactors[j],result,length);
	      
	  count = 0;
	  allone =0;
	  for (i=0; i<length; i++)
	    if (result[i]==locuslist[locus].type.binary.binfactors[j][i])
	      count +=1;
	    else if (result[i]==1)
	      allone +=1;
	      
	  if ((count == length)||(allone == length))
	    {
	      a[where] = j;
	      where += 1;
	    }
	}
	    
      if (where ==1)
	a[1] = a[0];

      k=0;
      for (i=0; i<locuslist[locus].nalleles; i++)
	for (j=i; j<locuslist[locus].nalleles; j++)
	  {
	    inversehap[i][j] = k+1;         /*note inversehap[0][0] indicates AA*/
	    inversehap[j][i] = k+1;         /*and inversehap[0][1] = AB*/
	    k+=1;
	  }
      *aff=2;
      *liab=inversehap[a[0]][a[1]];
    }
}

void numaff(all1, all2,aff,liab,locus)
     int all1,all2,*aff,*liab,locus;
{
  
  if ((all1==0)&(all2==0))
    {
      *aff = 0;
      *liab = 1;
    }
  else
    {
      *liab = table[locus][all1-1][all2-1];
      *aff = 2;
    }
}
 
main()
{
  int i, j, k, l,m,n,count,o,tmpi;
  double tmpd;
  FILE  *ped, *newped, *dat, *newdat;

  /*outped*/
  int numind, numloci,numcom,indv,child3,mom,dad,child1,child2,sex,proband;
  int all1, all2,obsbin[MAXCOM],pedno,numfactors,liabclasses,locitype,aff;
  int liab,numalleles,lociorder[MAXLOCI];

  /*outdat*/
  int risk, sexlink,progcode,mutlocus,sexdiff,interfere,itlocus1,itlocus2,len,hapfreq;
  double mutmale,mutfeme,recombfrac;
  char line[MAXLINE];

  /*outdat*/
  dat = fopen("indat","r");
  newdat = fopen("outdat","w");

  fscanf(dat,"%d %d %d %d",&numloci,&risk,&sexlink,&progcode);
  getline(line,MAXLINE,dat);
  fprintf(newdat,"%d %d %d %d\n",numloci,risk,sexlink,progcode);
  fscanf(dat,"%d %lf %lf %d",&mutlocus, &mutmale, &mutfeme,&hapfreq);
  getline(line,MAXLINE,dat);
  fprintf(newdat, " %d %8.6lf %8.6lf %d\n ",mutlocus,mutmale,mutfeme,hapfreq);
  
  for (i=0; i<numloci; i++)
    {
      fscanf(dat,"%d",&tmpi);
      lociorder[i]=tmpi;
      fprintf(newdat,"%d ",lociorder[i]);
    }
  getline(line,MAXLINE,dat);
  fprintf(newdat,"\n\n");
  
  for (j=0; j<numloci; j++)
    {
      fscanf(dat,"%d %d",&locitype,&numalleles);
      getline(line,MAXLINE,dat);
      fprintf(newdat,"%d %d\n",1,numalleles);

      locuslist[j].nalleles = numalleles;
      locuslist[j].format = locitype;      
      
      for (k=0; k<numalleles; k++)
	{
	  fscanf(dat,"%lf",&tmpd);
	  locuslist[j].freq[k] = tmpd;
	  fprintf(newdat,"%6.4lf ",tmpd);
	}
      fprintf(newdat,"\n ");
      getline(line,MAXLINE,dat);
      
      numcom=((numalleles+1)*numalleles)/2;
	      
      switch (locitype)
	{
	case 1: 
	  {
	    fscanf(dat,"%d",&liabclasses);
	    fprintf(newdat,"%d\n ",liabclasses);
	    getline(line,MAXLINE,dat);

	    locuslist[j].type.aff.nclass=liabclasses;
	    
	    for (l=0; l<liabclasses; l++)
	      {
		for (n=0; n<numcom; n++)
		  {
		    fscanf(dat,"%lf",&tmpd);
		    locuslist[j].type.aff.pene[l][n]=tmpd;
		    fprintf(newdat,"%8.6lf ",tmpd);
		  }
		getline(line,MAXLINE,dat);
		fprintf(newdat,"\n ");
	      }
	    fprintf(newdat,"\n");
	  }
	  break;
	case 2:
	  {
	    fscanf(dat,"%d",&numfactors);
	    getline(line,MAXLINE,dat);
	    fprintf(newdat,"%d\n ",numcom);

	    locuslist[j].type.binary.numfact=numfactors;

	    for (l=0; l<numalleles; l++)
	      {
		for (n=0; n<numfactors; n++)
		  {
		    fscanf(dat,"%d",&tmpi);
		    locuslist[j].type.binary.binfactors[l][n]=tmpi;
		  }
		getline(line,MAXLINE,dat);
	      }

	    for (n=1; n<=numcom; n++)
	      {
		for (m=1; m<=numcom; m++)
		  {
		    if (n!=m)
		      fprintf(newdat,"%5.3lf ",0.0);
		    else
		      fprintf(newdat,"%5.3lf ",1.0);
		  }
		fprintf(newdat,"\n ");
	      }
	    fprintf(newdat,"\n");
	  }
	  break;
	case 3:
	  {
	    fprintf(newdat,"%d\n ",numcom);
	    for (n=1; n<=numcom; n++)
	      {
		for (m=1; m<=numcom; m++)
		  {
		    if (n!=m)
		      fprintf(newdat,"%5.3lf ",0.0);
		    else
		      fprintf(newdat,"%5.3lf ",1.0);
		  }
		fprintf(newdat,"\n ");
	      }
	    fprintf(newdat,"\n");
	  }
	  break;
	}
    }

  fscanf(dat,"%d %d",&sexdiff, &interfere);
  fprintf(newdat," %d %d ",sexdiff, interfere);

  while (( len = getline(line, MAXLINE, dat)) > 0)
    putline(line, len, newdat);

  fclose(dat);
  fclose(newdat);


  /*outped file creation */
  ped = fopen("inped","r");
  newped = fopen("outped","w");

  /*Create Table*/
  for (i=0; i<MAXLOCI; i++)
    for (k=0; k<MAXALL; k++)
      for (j=0; j<MAXALL; j++)
	table[i][j][k] = 0;

  for (i=0; i<numloci; i++)
    {
      count = 1;
      for (k=0; k<locuslist[i].nalleles; k++)
	for (j=k; j<locuslist[i].nalleles; j++)
	  {
	    table[i][k][j] = count;
	    table[i][j][k] = count;
	    count += 1;
	  }
    }
  
    
  while (fscanf(ped,"%d ",&pedno)!=EOF)
    {
      fscanf(ped, "%d  %d  %d  %d  %d  %d  %d %d",&indv,&mom,&dad,&child1,&child2,&child3,&sex,&proband);

      fprintf(newped, "%4d  %3d  %3d  %3d  %3d  %3d  %3d %d  %d  ",pedno,indv,mom,dad,child1,child2,child3,sex,proband);

      for (k=0; k<numloci; k++)
	{
	  switch (locuslist[k].format)
	    {
	    case 1:
	      {
		if (locuslist[k].type.aff.nclass==1)
		  {
		    fscanf(ped,"%d ",&aff);
		    fprintf(newped,"%2d  ",aff);
		  }
		else
		  {
		    fscanf(ped,"%d %d",&aff,&liab);
		    fprintf(newped,"%2d %2d  ",aff,liab);
		  }
	      }
	      break;
	    case 2:
	      {
		for (o=0; o<locuslist[k].type.binary.numfact; o++)
		  {
		    fscanf(ped,"%d ",&tmpi);
		    obsbin[o]=tmpi;
		  }

		binaff(obsbin,&aff,&liab,k,locuslist[k].type.binary.numfact);
		fprintf(newped,"%2d %2d  ",aff,liab);
	      }
	      break;
	    case 3:
	      {
		fscanf(ped," %d %d",&all1, &all2);
		numaff(all1,all2,&aff,&liab,k);
		fprintf(newped,"%2d %2d  ",aff,liab);
	      }
	      break;
	    }
	  
	}
      getline(line,MAXLINE,ped);
      fprintf(newped,"\n");
    }
      
  fclose(ped);
  fclose(newped);

}