// ----------------------------------------------------------------------
//
// File:        Samples.c
//
// Description: 
//
// Date:        January 29, 2004
//
// Author(s):   Giri Narasimhan, Florida International University 
//
// ----------------------------------------------------------------------

#include "Samples.h"

/*
struct individual
{
  float t1b1, t1b2, t1g1, t1g2, t2b1, t2b2, t2g1, t2g2;
  int isHom; // 1 for homozygous and 0 for heterozygous
};

const int MAXSAMPLES = 100;

class Samples{
      
  individual HoS[MAXSAMPLES];
  individual HeS[MAXSAMPLES];
  int SampleCount, HoSCount, HeSCount;
*/

int Samples::getSampleCount ()
{ return SampleCount;
}

int Samples::getHoSCount ()
{ return HoSCount;
}

int Samples::getHeSCount ()
{ return HeSCount;
}

float Samples::getThreshold ()
{ return TRACETHRESHOLD;
}

void Samples::setThreshold (float f)
{ TRACETHRESHOLD = f;
}

void Samples::ReadOneSample (ifstream& fIn)
{
  individual ind;
  /*
  for (int i =0; i<NAMELENGTH; i++)
	  cin.get(ind.name[i]);
  cin.width(6);
  cin >> ind.name;
  */
  fIn >> ind.t1b1 >> ind.t1b2 >> ind.t1g1 >> ind.t1g2 
      >> ind.t2b1 >> ind.t2b2 >> ind.t2g1 >> ind.t2g2;
  // cout << "Hello2\n";
  fIn.get(ind.name, NAMELENGTH, '\r');
  // cout << "Hello3\n";
  if (DEBUG) 
    {
      cout  << ind.name
	    << ": "  << ind.t1b1 << ','  << ind.t1b2 
	    << ','  << ind.t1g1 << ','  << ind.t1g2
	    << ','  << ind.t2b1 << ','  << ind.t2b2 
	    << ','  << ind.t2g1 << ','  << ind.t2g2;
      if (ind.isHom) cout << "  Homozygous" << endl;
      else cout << endl;
    }
  ind.isHom = testHomozygosity(ind);
  if (ind.isHom) 
    {
      HoS[HoSCount++] = ind;
      cout << "Homozygous Sample: " << HoSCount << ' ' << ind.name << endl;
    }
  else 
    {
      HeS[HeSCount++] = ind;
      cout << "Heterozygous Sample: " << HeSCount << ' ' << ind.name << endl;
    }
  SampleCount++; 
  if (DEBUG) 
    {
      cout  << ind.name
	    << ": "  << ind.t1b1 << ','  << ind.t1b2 
	    << ','  << ind.t1g1 << ','  << ind.t1g2
	    << ','  << ind.t2b1 << ','  << ind.t2b2 
	    << ','  << ind.t2g1 << ','  << ind.t2g2;
      if (ind.isHom) cout << "  Homozygous" << endl;
      else cout << endl;
    }
}

void Samples::ReadAllSamples (ifstream& fIn)
{
  cout << "========================================================\n";
  if (DEBUG) cout << "Reading all samples" << endl;
  char temp[500];
  fIn.getline(temp, 500, '\n'); // read and discard header line
  char c;
  while ((c = fIn.get()) != EOF)
    {
      // cout << "Hello1\n";
      ReadOneSample(fIn);
    }
  cout << "Finished reading " << SampleCount << " samples ("
       << HoSCount << " homozygous and " << HeSCount 
       << " heterozygous samples)" << endl;
}

void Samples::ReadOneHSample (ifstream& fIn)
{
  individual ind;
  /*
  for (int i =0; i<NAMELENGTH; i++)
	  cin.get(ind.name[i]);
  cin.width(6);
  cin >> ind.name;
  */
  fIn >> ind.t1b1 >> ind.t1g1 >> ind.t2b1 >> ind.t2g1;
  // fIn.get(ind.name,NAMELENGTH, '\n');
  // if (DEBUG) cout << ind.name;
  ind.t1b2 = ind.t1b1; ind.t1g2 = ind.t1g1;
  ind.t2b2 = ind.t2b1; ind.t2g2 = ind.t2g1;
  ind.isHom = true;
  HoS[HoSCount++] = ind;
  SampleCount++; 
  if (DEBUG) 
    {
      cout  << ind.t1b1 << "; "  << ind.t1g1 << "; "  
	    << ind.t2b1 << "; "  << ind.t2g1 ;
      if (ind.isHom) cout << "  Homozygous" << endl;
      else cout << endl;
    }
}

void Samples::ReadAllHSamples (ifstream& fIn)
{
  cout << "========================================================\n";
  if (DEBUG) cout << "Reading all samples" << endl;
  char temp[500];
  fIn.getline(temp, 500, '\n'); // read and discard header line
  char c;
  while ((c = fIn.get()) != EOF)
    {
      ReadOneHSample(fIn);
    }
  cout << "Finished reading " << HoSCount << " homozygous samples." << endl;
  if (DEBUG) 
    {
      cout  << HoS[0].t1b1 << "; "  << HoS[0].t1g1 << "; "  
	    << HoS[0].t2b1 << "; "  << HoS[0].t2g1 << endl;
    }
}

bool Samples::testHomozygosity (individual ind)
{
  // cout << "In testHo" << endl;
  if ((ind.t1b1==ind.t1b2) && (ind.t1g1==ind.t1g2) &&
      (ind.t2b1==ind.t2b2) && (ind.t2g1==ind.t2g2)) 
    return true;
  else return false;
}

float Samples::computeTrace (individual ind1, individual ind2)
{ 
  /* (DEBUG) cout << ind1.t1b1 << ind2.t1b1 << ind1.t1g1 << ind2.t1g1
     << ind1.t2b1 << ind2.t2b1 << ind1.t2g1 << ind2.t2g1 << endl; */
  float Trace;
  if ((ind1.isHom) && (ind2. isHom))
    Trace = 
      fabs(ind1.t1b1 - ind2.t1b1) + fabs(ind1.t1g1 - ind2.t1g1) +
      fabs(ind1.t2b1 - ind2.t2b1) + fabs(ind1.t2g1 - ind2.t2g1);
  else if (((ind1.isHom) && (!(ind2.isHom))) || ((!(ind1.isHom)) && (ind2.isHom)))
    Trace =
      MIN(fabs(ind1.t1b1 - ind2.t1b1), fabs(ind1.t1b2 - ind2.t1b2)) + 
      MIN(fabs(ind1.t1g1 - ind2.t1g1), fabs(ind1.t1g2 - ind2.t1g2)) +
      MIN(fabs(ind1.t2b1 - ind2.t2b1), fabs(ind1.t2b2 - ind2.t2b2)) + 
      MIN(fabs(ind1.t2g1 - ind2.t2g1), fabs(ind1.t2g2 - ind2.t2g2));
  else // Hetero-Hetero case: this is complicated. Need to think about it
    Trace = 100000;
  return Trace;
}

int Samples::HomoSameAllele (individual hoP1, individual hoP2)
{
  float Trace = computeTrace(hoP1, hoP2);
  if (Trace <= TRACETHRESHOLD)
    return true;
  else return false;
}

int Samples::HomoHeteroSameAllele (individual hoP1, individual heP2)
{
  float Trace = computeTrace(hoP1, heP2);
  if (Trace <= TRACETHRESHOLD)
    return true;
  else return false;
}

int Samples::HeteroHeteroCommonAllele (individual heP1, individual heP2)
{
  return true; // not implemented yet!
}

float Samples::computeThreshold (ofstream& fOut)
{
  if ((HoSCount == SampleCount) && (SampleCount > 0)) 
    {
      float MaxThreshold, Max2 = 0;
      int tempI, tempJ, tempI2, tempJ2;
      if (DEBUG)
	{
	  cout << endl 
	       << "========================================================" 
	       << endl;;
	  cout << "Computing Threshold" << endl;
	}
      fOut << "Sample 1\t Sample 2\t Trace" << endl;
      for (int i = 0; i<HoSCount; i++)
	for (int j = i+1; j<HoSCount; j++)
	  { 
	    float temp = computeTrace(HoS[i], HoS[j]);
	    fOut << i << '\t' << j << '\t' << temp << endl;
	    if (DEBUG) 
	      {
		cout << i << ' ' << j << ' ' << temp << "; ";
		if ((!(j == 0)) && (j%5 == 0)) cout << endl;
	      }
	    if (temp > MaxThreshold) 
	      {
		Max2 = MaxThreshold;
		tempI2 = tempI; tempJ2 = tempJ;
	        MaxThreshold = temp;
		tempI = i; tempJ = j;
	      }
	  }
      cout << "Threshold = " << MaxThreshold << " for indices "
	      << tempI << " and " << tempJ << endl;
      cout << "2nd highest Threshold = " << Max2 << " for indices "
	      << tempI2 << " and " << tempJ2 << endl;
      fOut << "Threshold = " << MaxThreshold << " for indices "
	      << tempI << " and " << tempJ << endl;
      fOut << "2nd highest Threshold = " << Max2 << " for indices "
	      << tempI2 << " and " << tempJ2 << endl;
      fOut.close();
      return MaxThreshold;
    }
  cout << "Cannot compute threshold." << endl;
  return -1;
}

void Samples::ProcessAllSamples ()
{
  cout << "========================================================\n";
  cout << "Processing all homozygous samples" << endl;
  bool Ho[HoSCount], He[HeSCount];
  for (int i=0; i<HoSCount; i++) Ho[i] = false;
  for (int i=0; i<HeSCount; i++) He[i] = false;
  for (int i = 0; i<HoSCount; i++)
    if (!Ho[i])
      {
	cout << i+1 << " (" << HoS[i].name 
	     << ") has same allele as: ";
	for (int j = i+1; j<HoSCount; j++)
	  {
	    if (computeTrace(HoS[i], HoS[j]) <= TRACETHRESHOLD)
	      {
		cout << j+1 << "; "; // " (" << HoS[j].name << "); ";
		Ho[j] = true;
	      }
	  }
      cout << endl;
    }
  cout << "Processing all heterozygous samples" << endl;
  for (int j = 0; j<HeSCount; j++)
    {
      cout << j+1 << " (" << HeS[j].name 
	   << ") has same allele as homozygous samples: ";
      for (int i = 0; i<HoSCount; i++)
	if (!Ho[i])
	  {
	    if (computeTrace(HoS[i], HeS[j]) <= TRACETHRESHOLD)
	      cout << i+1 << "; "; // << " (" << HoS[i].name << "); ";
	  } 
      cout << endl;
    }
}

void Samples::OutputResults ()
{
  cout << "========================================================\n";
	cout << "Results: " << endl;
}

void Samples::setDEBUG(bool b)
{
  DEBUG = b;
}

