#include "affinityScore.h"


// Computes the global affinity score for the proteins in msa
//   using the subsitution Matrix subMat
AffinityScore::AffinityScore(Alignment msa, Matrix subMat) {
  int i, j, k, A, B=0;
  // int n, totNumRes=0;
  length = msa.maximumSequenceLength;
  val = new float[msa.maximumSequenceLength];
  for(i=0;i<msa.maximumSequenceLength;i++) val[i]=0;
  subVal = NULL;
  subName = NULL;

  for(k=0; k<msa.nSequences; k++)
    for(j=0; j<k; j++)
      for(i=0; i<msa.maximumSequenceLength; i++) {
        A = subMat.residueCharToInt(msa.sequences[k].residues[i]);
        B = subMat.residueCharToInt(msa.sequences[j].residues[i]);
        if(A>-1 && B>-1) {
          val[i]+=subMat.cell[A][B];
        } else {
          #if ADD_GAP_PENALTY
            if(A>-1 || B>-1) val[i]+=.5*subMat.minimumCellValue;
          #endif
        }
      }

  // Normalize Scores
  for(j=0; j<i; j++)
    val[j] *= (float)1/(msa.nSequences*(msa.nSequences-1));

}

      
// Computes subMatrices for all proteins in the set
//   and creates modified pdb file for all
void AffinityScore::doAllSub(Alignment msa, Matrix subMat) {
  if(DB) printf("entering AffinityScore::doAllSub\n");
  for(int i=0; i<msa.nSequences; i++) {
    #if !DO_NC2
      getSubN(msa, subMat, msa.sequences[i].name);
    #endif
    getSub(msa, msa.sequences[i].name);
    printSub();
    modifyPDB();
  }
  if(DB) printf("leaving AffinityScore::doAllSub\n");
}


// Computes some stats for the affinity score
void AffinityScore::getStats() {
  float tot=0, sumDiff=0;

  for(int i=0; i<length; i++) {
    tot += val[i];
  }
  avg = tot/length;
  for(int i=0; i<length; i++) {
    sumDiff += (val[i]-avg)*(val[i]-avg);
  }
  stdDev = sqrt(sumDiff/length);

  printf("stdDev = %f\n", stdDev);

}


// Computes the affinity score for a particular protein (name)
void AffinityScore::getSub(Alignment msa, char* name) {
  subVal = new float[msa.maximumSequenceLength];
  int i;
  for(i=0; i<msa.nSequences; i++)
    if(strncmp(msa.sequences[i].name,name,5)==0)
      break;
  int k=0;
  for(int j=0; j<length; j++)
    if(msa.sequences[i].residues[j] != '-') {
      printf("%c", msa.sequences[i].residues[j]);
      subVal[k++] = val[j];
    }
  printf("\n\n");
  subLength = k;
}


// Computes the n-wise affinity score for a particular protein (name)
void AffinityScore::getSubN(Alignment msa, Matrix subMat, char* name) {
  if(DB) printf("entering AffinityScore::getSubN\n");
  int i, j, k, n, A, B=0;
  // int totNumRes=0;

  subVal = new float[msa.maximumSequenceLength];
  for(i=0;i<msa.maximumSequenceLength;i++) subVal[i]=0;
  if(subName!=NULL) delete[] subName;
  subName = new char[strlen(name)];
  strncpy(subName, name, 5);
  subName[5]='\0';
  for(i=0; i<msa.nSequences; i++)
    if(strncmp(msa.sequences[i].name,name,5)==0)
      break;
  n = i;
  //int count=0;
  for(k=0; k<msa.nSequences; k++) {
    if(k!=n) {
      for(i=0; i<msa.maximumSequenceLength; i++) {
//  for(j=0;j<msa.maximumSequenceLength;j++)
//    printf("%d %d %f\n",i,j,val[j]);
        A = subMat.residueCharToInt(msa.sequences[k].residues[i]);
        B = subMat.residueCharToInt(msa.sequences[n].residues[i]);
        if(A>-1 && B>-1) {
          val[i]+=subMat.cell[A][B];
        } else {
          if(A>-1 || B>-1) 
            if(ADD_GAP_PENALTY)
              val[i]+=.5*subMat.minimumCellValue;
        }
      }
    }
  }

  // Normalize Scores
  for(j=0; j<i; j++)
    val[j] *= (float)1/msa.nSequences;
  printf("msa.nSequences = %d\n",msa.nSequences);
  printf("msa.maximumSequenceLength = %d\n",msa.maximumSequenceLength);

  for(i=0,k=0; i<msa.maximumSequenceLength; i++) {
    if(subMat.residueCharToInt(msa.sequences[n].residues[i])>-1) {
//      printf("%c %d %f \n", msa.sequences[n].residues[i],k+1, val[i]);
//      printf("%d %f \n", k+1, val[i]);
      k++;
    }
  }
  if(DB) printf("leaving AffinityScore::getSubN\n");
}


// Copies the pdb file corresponding to AffinityScore::name and modifies its occupancy field
//   so that it is floor(100*AffinityScore) for each residue
void AffinityScore::modifyPDB()
{
  FILE * oPdbFile, *nPdbFile;
  char *temp = new char[MSL];
  char *temp1 = new char[MSL];
  char *temp2 = new char[MSL];
  int n, old;
  strcpy(temp, "cp ");
  strcat(temp, PDB_PATH); 
  strcat(temp, lower(subName)); 
  strcat(temp, "* "); 
  strcat(temp, lower(subName)); 
  strcat(temp, ".pdb"); 

  system(temp);
  strcat(temp, "2"); 
  system(temp);

  strcpy(temp1, lower(subName));
  strcat(temp1, ".pdb2");
  strcpy(temp2, lower(subName));
  strcat(temp2, ".pdb");

  oPdbFile = fopen(temp1, "r");
  nPdbFile = fopen(temp2, "w");

  old = -99;
  int count = 0;
  while(fgets(temp, 500, oPdbFile)!=NULL) {
    if(strncmp(temp, "ATOM", 4)!=0) {
      fputs(temp, nPdbFile);
    }
    else {
      strcpy(temp1,temp+23);
      temp1[4] = '\0';
      n = charToInt(temp1);
      if(n!=old) {
        old=n;
        count++;
      }
//      strcpy(temp1, intToString((int)floor(100*subVal[count-1]),5));
//      strcpy(temp1, intToString(normVal(subVal[count-1]),5));
      strcpy(temp1, intToString(normVal(val[count-1]),5));
      for(int i=0; i<5; i++)
        temp[55+i] = temp1[i];
      fputs(temp, nPdbFile);
    }
  }

  strcpy(temp2, "rm -f ");
  strcat(temp2, lower(subName));
  strcat(temp2, ".pdb2");
  system(temp2);

  delete[] temp;
  delete[] temp1;
  delete[] temp2;
}

  
// Computes the some stats for the affinity score
int AffinityScore::normVal(float f) {
  if(f<-3*stdDev)
    return (int)floor(100*(-3*stdDev + avg));
  if(f>3*stdDev)
    return (int)floor(100*(3*stdDev + avg));
  return (int)floor(100*f);
}


// Prints the global affinity score calculated in the constructor
void AffinityScore::print() {
  for(int i=0; i<length; i++) {
    printf("%d %f \n", i+1, val[i]);
  }
}


// Prints affinity score for a particular protein (AffinityScore::name)
void AffinityScore::printSub() {
  for(int i=0; i<subLength; i++) {
    printf("%d %f \n", i+1, subVal[i]);
  }
}