#include "includes.h"
#include "grid.h"
#include "vector.h"

#include "declare.h"

/* This file contains source code for functions called in function max()
   which estimates ancestral haplotype */

void makeinitcandanc(igrid &candanc, const int &C, 
		     const int &LE, const int &RE, const int &ctrall, 
		     const igrid &allelelist, int &numcand)
{
  //This function generates a first list of candidate ancestral haplotypes;
  //the haplotypes consist of two adjacent markers and locus C. All 
  //possibilities are listed.

  int i,j;//make variables for loops;
  int numCm1=0;//find out how many alleles at locus LE;
  int numCp1=0;//find out how many alleles at locus RE;
  
  //find out how many alleles at LE;
  for (i=0;i<allelelist.dimx();i++)
    {
      if (i<allelelist.dimx()-1 && allelelist(i+1,LE)==0) {numCm1=i+1;break;}
    }
  if (numCm1==0) numCm1=allelelist.dimx();
 
  //find out how many alleles at C+1;
  for (i=0;i<allelelist.dimx();i++)
    {
      if (i<allelelist.dimx()-1 && allelelist(i+1,RE)==0) {numCp1=i+1;break;}
    }
  if (numCp1==0) numCp1=allelelist.dimx();

  //number of rows used in matrix;
  numcand=numCm1*numCp1;
  if (numcand > candanc.dimx()) 
    {
      cout << "ERROR: Need more than " << candanc.dimx() << 
	" rows in candanc matrix" << endl;
      exit(1);
    }
  
  
  //make center column;
  for (i=0;i<numCm1*numCp1;i++) candanc(i,C)=ctrall;
        
  //make column at locus LE;
  for (i=0;i<numCm1;i++)
    {
      for (j=0;j<numCp1;j++)
	{
	  candanc(i*numCp1+j,LE)=allelelist(i,LE);
	}
    }
  
  //make column at locus RE;
  for (i=0;i<numCm1;i++)
    {
      for (j=0;j<numCp1;j++)
	{
	  candanc(i*numCp1+j,RE)=allelelist(j,RE);
	}
    }
  
}//makeinitcandanc

void makecandanc(igrid &candanc, const igrid &allelelist, int &dirjump, 
		 int &numcand, int &numkeep, const int &C, int &LE, int &RE,
		 const int &ctrall)
{
  //Given a list of n candidate ancestral haplotypes, this function adds
  //each allele at an adjacent marker to this list, yielding n*m candidates,
  //where m is the number of alleles at the added locus

  int i,j,k;//looping indices;
  int numCq1=0; int Cq1;

  if (dirjump==-1) { if (LE-1==C) Cq1=LE-2; else Cq1=LE-1; }
  else { if (RE+1==C) Cq1=RE+2; else Cq1=RE+1; }

  //find out how many times to copy core of haps to keep, i.e.,
      //find out how many alleles at Cq1;
  for (i=0;i<allelelist.dimx();i++)
    {
      if (i<allelelist.dimx()-1 && allelelist(i+1,Cq1)==0) {numCq1=i+1;break;}
    }
  if (numCq1==0) numCq1=allelelist.dimx();

  //find out number of rows used in candres;
  numcand=numkeep*numCq1;
  if (numcand > candanc.dimx()) 
    {
      cout << "ERROR: Need more than " << candanc.dimx() << 
	" rows in candanc matrix" << endl;
      exit(1);
    }

  //copy core numCq1-1 times;
  for (i=1;i<=numCq1-1;i++)//copy index
    {
      for (j=LE;j<=RE;j++)//locus index
	{
	  for (k=0;k<numkeep;k++)//core row
	    {
	      candanc(i*numkeep+k,j)=candanc(k,j);
	    }
	}
    }
  
  //make the jump column;
  for (i=0;i<numCq1;i++)
    {
      for (j=0;j<numkeep;j++)
	{
	  candanc(j+i*numkeep,Cq1)=allelelist(i,Cq1);
	}
    }


  //fill in C if C is not between LE and RE (1/12/00)
  if (C<LE || C>RE)
    {
      for (j=0;j<numcand;j++)
	{
	  candanc(j,C)=ctrall;
	}
    }

}//makecandanc

void sortcandancbylik(igrid &candanc, dgrid &candres, 
		      const int &numcand, const int &loclik)  
{
  //this function orders the candidate ancestral haplotypes and their
  //corresponding array of results by decreasing loglikelihood

  //note that LE and RE have already been adjusted;

  int row, row1;

  //sort;
  for (row=0; row<numcand-1; row++)
    {
      for (row1=row+1; row1<numcand; row1++)
	{
	  if(candres(row1,loclik) > candres(row,loclik))
	    {
	      rowswitchint(candanc,row,row1);
	      rowswitchdou(candres,row,row1);
	    }
	}
    }
}//sortcandancbylik

void rowswitchint(igrid &candanc, const int &row, const int &row1)
{
  //this function switches two rows in a grid of integers

  int j;

  ivector temprow(candanc.dimy());
  for (j=0;j<candanc.dimy();j++) temprow[j]=candanc(row,j);
  for (j=0;j<candanc.dimy();j++) candanc(row,j)=candanc(row1,j);
  for (j=0;j<candanc.dimy();j++) candanc(row1,j)=temprow(j);
}//rowswitchint

void rowswitchdou(dgrid &candres, const int &row, const int &row1)
{
  //this function switches two rows in a grid of doubles

  int j;

  dvector temprow(candres.dimy());
  for (j=0;j<candres.dimy();j++) temprow[j]=candres(row,j);
  for (j=0;j<candres.dimy();j++) candres(row,j)=candres(row1,j);
  for (j=0;j<candres.dimy();j++) candres(row1,j)=temprow(j);
}//rowswitchdou

int whichdir(const int &LE, const int &RE, const dvector &x)
{  
  //this function determines whether the next locus added to the 
  //ancestral haplotype will be to the right or left; the decision
  //is based on picking the next marker closest to the gene--this
  //marker presumably contains more info about the parameters than
  //any other remaining choice

  int dir;
  if (LE>0 && RE<x.dim()-1)
    {
      if (x[LE-1] < x[RE+1]) dir = -1;
      else dir = 1;
    }  
  else
    {
      if (LE<=0 && RE>=x.dim()-1) dir = 0;
      else
	{    
	  if (LE <= 0) dir = 1;
	  else dir = -1;
	}
    }
  return dir;
}//whichdir

int findnumkeep(const dgrid &candres, const int &dirjump, 
		const double &threshfix, const double &threshvar, 
		int &numcand)
{
  //when estimating the extent of LD, not mapping, the number of candidate
  //ancestral haplotypes designated to continue being "grown" depends on
  //the magnitude of cstar. This function counts the number of candidates
  //fed into makecandanc().

  //find column by which weve sorted;
  int pos; if (dirjump==-1) pos=6; else pos=7;

  int i; int dummy=1; int numkeep=0;
      
  //count number of haplotypes above the threshold;
  for (i=0;i<numcand;i++)
    if(candres(i,pos)<threshfix || candres(i,pos)<threshvar) 
      {
	numkeep=i;
	dummy=0;
	break;
      }
  if (dummy==1) numkeep=numcand;
  return numkeep;
}//findnumkeep

void cleantemps(igrid &candanc, dgrid &candres, const int &numcand, 
		const int &numkeep, const int &LE, const int &RE) 
{
  //this function deletes discarded candidates in candanc

  int i;
  int j;
  
  for (i=numkeep;i<numcand;i++)
    {
      for (j=LE;j<=RE;j++) candanc(i,j)=0;
      for (j=0;j<candres.dimy();j++) candres(i,j)=0;
    }
}//cleantemps

void correctlik(dgrid &candres, const int &numcand, const int &loclik)
{
  //in earlier versions of the software, a multiple of log(2) was added 
  //to the model and null-loglikelhoods for some candidate ancestral
  //haplotypes. We added multiples of log(2) to each to make them comparable.

  int i;
  double templikmod=0; double templiknull=0;
    
  //find maximum null loglik
  double max=candres(0,loclik+1);
  for (i=1;i<numcand; i++)
    {
      if (candres(i,loclik+1)>max) max=candres(i,loclik+1);
    }
  
  //equalize null lik's and add to model lik.
  for (i=0;i<numcand && candres(i,loclik+1) != 0; i++)
    {
      if (candres(i,loclik+1) < max - 0.001)
	{
	  templikmod=candres(i,loclik);
	  templiknull=candres(i,loclik+1);
	  while ( candres(i,loclik+1) < max - 0.0001 )
	    {
	      candres(i,loclik)=candres(i,loclik)+log(2);
	      candres(i,loclik+1)=candres(i,loclik+1)+log(2);
	    }
	  if (candres(i,loclik+1) < max - 0.0001 || 
	      candres(i,loclik+1) > max + 0.0001) 
	    cout << "WARNING: the null likelihood for the " << i << 
	      "th candidate ancestral haplotype (" << templiknull << 
	      ") does not differ from the max null lik (" << max <<
	      ") by a multiple of log(2)" << endl;
	}
    }
}//correctlik