#include "alignment.h"
#include "typeConvert.h"
//#include "aaTools.h"

// default constructor
Alignment::Alignment() {
  if(DB) printf("entering default Alignment constructor\n");

  name = new char[30];
  strcpy(name, "default");

  nSequences = 0;
  nProteins = 0;
  maximumSequenceLength = 0;
  sequences = 0;

  if(DB) printf("leaving default Alignment constructor\n");
}

// constructor, assigns Alignment name to char* s
Alignment::Alignment(char *s) {
  if(DB) printf("entering Alignment char* constructor\n");

  if(strlen(s)>30) {
    printf("error in Alignment::Alignment(char*):  Alignment::name too long\n");
    return;
  }
  
  name = new char[30];
  strcpy(name, s);

  nSequences = 0;
  nProteins = 0;
  maximumSequenceLength = 0;
  sequences = 0;

  if(DB) printf("leaving Alignment char* constructor\n");
}


// Destructor
Alignment::~Alignment() {
  if(DB) printf("entering Alignment destructor\n");
  if(name!=0)
    delete[] name;
  if(sequences!=0)
    delete[] sequences;
  if(DB) printf("leaving Alignment destructor\n");
}


// copy constructor
Alignment::Alignment(const Alignment& m) {
  if(DB) printf("entering Alignment copy constructor\n");
  if(m.name==0) name=0;
  else {
    name = new char[strlen(m.name)+1];
    strcpy(name, m.name);
  }
  nSequences = m.nSequences;
  nProteins = m.nProteins;
  maximumSequenceLength = m.maximumSequenceLength;
  if(m.sequences==0) sequences=0;
  else {
    sequences = new Sequence[nSequences];
    for(int i=0; i<nSequences; i++)
      sequences[i] = m.sequences[i];
  }
  if(DB) printf("leaving Alignment copy constructor\n");
}


// Assignment
Alignment& Alignment::operator= (const Alignment& m) {
  if(DB) printf("entering Alignment operator=\n");
  if(this != &m) {
    if(m.name==0) name=0;
    else {
      if(name!=0) delete[] name;
      name = new char[strlen(m.name)+1];
      strcpy(name, m.name);
    }
    //int oldNSequences = nSequences; // to make sure we delete memory correctly
    //int oldMaxSeqLen = maximumSequenceLength;
    nSequences = m.nSequences;
    nProteins = m.nProteins;
    maximumSequenceLength = m.maximumSequenceLength;
    if(m.sequences==0) sequences=0;
    else {
      if(sequences!=0) 
        delete[] sequences;
      sequences = new Sequence[nSequences];
      for(int i=0; i<nSequences; i++)
        sequences[i] = m.sequences[i];
    }
  }
  if(DB) printf("leaving Alignment operator=\n");
  return *this;
}


/* outdated, do not use; still coded for a void::void func
Alignment Alignment::aaToSS(Alignment alignment) {
  FILE* fp;
  ssSeq = new char*[nSequences];
  for(int i = 0; i < nSequences; i++)
    ssSeq[i] = new char[sequences[i].length+1];
  char *temp1 = new char[MSL];
  char *temp2 = new char[MSL];
  int n, old, k;

  for(int i=0; i<nSequences; i++) {
    char *temp = new char[1200];
    strcpy(temp, "");
    strcat(temp, DSSP_PATH);
    strcat(temp, " -na ");
    strcat(temp, PDB_PATH);
    strcat(temp, sequences[i].name);
    strcat(temp, ".pdb > tempF");
    system(temp);
    delete[] temp;

    fp = fopen("tempF", "r");

    do {
      fgets(temp1, 200, fp);
    } while(temp1[2] != '#');

    k = 0;
    for(int j=0; j<sequences[i].length; j++) {
      if(sequences[i].residues[j]=='-')
        ssSeq[i][j] = '-';
      else {
        ssSeq[i][j] = temp1[16];
        if(sequences[i].residues[j] == ' ')
          ssSeq[i][j] = 'C';
        fgets(temp1, 200, fp);
      }
    }
//    printf("\n%s\n", ssSeq[i]);

  fclose(fp);

  }

  system("rm -f tempF");

  delete[] temp1; 
  delete[] temp2; 
}
*/


void Alignment::printFasta() {
  char *temp = new char[maximumSequenceLength+1];
  for(int j = 0; j < nSequences; j++) {
    printf(">%s\n", sequences[j].name);
    //printf(">%s\n", sequences[j].name, temp);
    for(int i = 0; i < maximumSequenceLength/60+1; i++) {
      strncpy(temp, sequences[j].residues+i*60,60);
      temp[60]='\0';
      if(strlen(temp)>0)
        printf("%s\n", temp);
    }
    if(strlen(temp)>0)
      printf("\n");
  }
  delete[] temp;
}


void Alignment::printMsa() {
  char *temp = new char[maximumSequenceLength+1];
  for(int i = 0; i < maximumSequenceLength/60+1; i++) {
    for(int j = 0; j < nSequences; j++) {
      strncpy(temp, sequences[j].residues+i*60,60);
      temp[60]='\0';
      if(strlen(temp)>0)
        printf("%s: %s\n", sequences[j].name, temp);
    }
    if(strlen(temp)>0)
      printf("\n");
  }
  delete[] temp;
}



void Alignment::printPwa() {
  if(DB) printf("entering Alignment::printPwa\n");
  char *temp= new char[maximumSequenceLength+1];

  for(int k = 0; k < nSequences/2; k++) {
    for(int i = 0; i < sequences[2*k].length/60+1; i++) {
      for(int j = 0; j < 2; j++) {
        strncpy(temp, sequences[2*k+j].residues+i*60,60);
        temp[60]='\0';
        if(strlen(temp)>0)
          printf("%s: %s\n", sequences[2*k+j].name, temp);
      }
      if(strlen(temp)>0)
        printf("\n");
    }
    printf("\n");
  }
  delete[] temp;
  if(DB) printf("leaving Alignment::printPwa\n");
}


// Reads in information from a PairWise or Multiple Sequence Alignment
//   Expects to find sequences written in FASTA format
void Alignment::readFromFile(char *s) {
  if(DB) printf("entering Alignment::readFromFile\n");

  // clear old data first
  nSequences = 0;
  nProteins = 0;
  maximumSequenceLength = 0;
  sequences = 0;

  //std::ifstream inFile;
  //inFile.open(s,std::ios::in);
  FILE *inFile = fopen(s, "r");

  int count=0;
  // int maxSeqLines=0;
  char *temp = new char[MSL+1];
  char *temp2 = new char[MSL+1];
  //while(inFile.getline(temp, MSL)) {
  while ( fgets(temp,MSL,inFile) ) {
    if (strstr(temp,"space")!=0) break;
    if (temp[0] == '>') {
      nSequences++;
      count=0; 
      printf("%d\n",nSequences);
    }
    else 
      count++; 
    //if(count>maxSeqLines)
    //  maxSeqLines = count;
  }
  fclose(inFile);
  //maximumSequenceLength = maxSeqLines*60;

  sequences = new Sequence[nSequences];

  //inFile.seekg(0, std::ios::beg);
  //inFile.clear();
  //inFile.seekg(0, std::ios::beg);

  //inFile.getline(temp,60+1);
  inFile = fopen(s, "r");
  fgets(temp,60+1,inFile);
  for (int i=0; i<nSequences; i++) {
    //sequences[i].setName(temp+1,7);
    for (int j=0; j<(int)strlen(temp); j++) {
      if (temp[j] == '\n') {
	temp[j] = '\0';
	break;
      }
    }
    sequences[i].setName(temp+1);
    //inFile.getline(temp,60+1);
    fgets(temp,60+1,inFile);
    do {
      for (int j=0; j<(int)strlen(temp); j++) {
	if (temp[j] == '\n') {
	  temp[j] = '\0';
	  break;
	}
      }
      strncat(temp2, temp, strlen(temp));
      //if(inFile.getline(temp,60+1)==0) break;
      fgets(temp,60+1,inFile);
      if (feof(inFile)>0) break;
    } while (temp[0] != '>');
    sequences[i].setResidues(temp2);
    strcpy(temp2,"\0");
  }

  fclose(inFile);

//  for(int i=0; i<nSequences; i++) {
//    for(int j=0; j<7; j++) {
//      sequences[i].name[j] = sequences[i].name[j+1];
//    }
//    sequences[i].name[6] = '\0';
//    sequences[i].name[7] = '\0';
//    if(sequences[i].name[5]=='(') sequences[i].name[5]='*';
//  }

  maximumSequenceLength=0;
  for(int i=0; i<nSequences; i++)
    if(sequences[i].length>maximumSequenceLength)
      maximumSequenceLength=sequences[i].length;

  // default for msa
  nProteins = nSequences;

  // check if Pwa, then nProteins is number of unique proteins
  if(nSequences>2)
    if(sequences[1].length!=sequences[2].length)
      nProteins = (int)floor(.1+.5*(1+sqrt((float)1+4*nSequences)));

  delete[] temp;
  delete[] temp2;

  if(DB) printf("leaving Alignment::readFromFile\n");
}

/*
// Reads in information from a PairWise or Multiple Sequence Alignment
//   Expects to find sequences written in FASTA format
void Alignment::readFromFile(char *s) {
  if(DB) printf("entering Alignment::readFromFile\n");

  // clear old data first
  nSequences = 0;
  nProteins = 0;
  maximumSequenceLength = 0;
  sequences = 0;

  std::ifstream inFile;
  inFile.open(s,std::ios::in);

  int count=0, maxSeqLines=0;
  char *temp = new char[MSL+1];
  char *temp2 = new char[MSL+1];
  while(inFile.getline(temp, MSL)) {
    if(strstr(temp,"space")!=0) break;
    if(temp[0] == '>') {
      nSequences++;
      count=0; 
      //printf("%d\n",nSequences);
    }
    else 
      count++; 
    if(count>maxSeqLines)
      maxSeqLines = count;
  }
  maximumSequenceLength = maxSeqLines*60;

  sequences = new Sequence[nSequences];

  inFile.seekg(0, std::ios::beg);
  inFile.clear();
  inFile.seekg(0, std::ios::beg);

  inFile.getline(temp,60+1);
  for(int i=0; i<nSequences; i++) {
    //sequences[i].setName(temp+1,7);
    sequences[i].setName(temp+1);
    inFile.getline(temp,60+1);
    do {
      strncat(temp2, temp, strlen(temp));
      if(inFile.getline(temp,60+1)==0) break;
    } while(temp[0] != '>');
    sequences[i].setResidues(temp2);
    strcpy(temp2,"\0");
  }

//  for(int i=0; i<nSequences; i++) {
//    for(int j=0; j<7; j++) {
//      sequences[i].name[j] = sequences[i].name[j+1];
//    }
//    sequences[i].name[6] = '\0';
//    sequences[i].name[7] = '\0';
//    if(sequences[i].name[5]=='(') sequences[i].name[5]='*';
//  }

  maximumSequenceLength=0;
  for(int i=0; i<nSequences; i++)
    if(sequences[i].length>maximumSequenceLength)
      maximumSequenceLength=sequences[i].length;

  // default for msa
  nProteins = nSequences;

  // check if Pwa, then nProteins is number of unique proteins
  if(nSequences>2)
    if(sequences[1].length!=sequences[2].length)
      nProteins = (int)floor(.1+.5*(1+sqrt((float)1+4*nSequences)));

  delete[] temp;
  delete[] temp2;

  if(DB) printf("leaving Alignment::readFromFile\n");
}
*/

// Read and parse pdb file
void Alignment::readCaCoordinates() {
  //FILE *fpError = fopen("pdbError","w");
  //fprintf(fpError,"");
  //fclose(fpError);
  for(int i=0;i<nSequences;i++) {
    sequences[i].structure.readFromPdbFile();
    if(sequences[i].lengthWithoutGaps!=sequences[i].structure.length) {
      printf("WARNING: sequence length (%d) != structure length (%d) for protein %s\n",
          sequences[i].lengthWithoutGaps,sequences[i].structure.length,
          sequences[i].name);
    }
  }
}




// Sets all sequences' isSecondary flags to true
void Alignment::setIsSecondary() {
  if(DB) printf("entering Alignment::setIsSecondary\n");
  for(int i=0; i<nSequences; i++)
    sequences[i].isSecondary=true;
  if(DB) printf("leaving Alignment::setIsSecondary\n");
}




// Initializes Sequence
/*Sequence::Sequence(char *s) {
  if(DB) printf("entering Sequence constructor \n");
  name = new char[strlen(s)+1];
  strcpy(name,s);

  length = 0;
  lengthWithoutGaps = 0;
  isSecondary = false;
  isGap = 0;
  isEndGap = 0;
  residues = 0;
  residuesWithoutGaps = 0;
  residuesIndexToResiduesWithoutGapsIndex = 0;
  residuesWithoutGapsIndexToResiduesIndex = 0;
  structure.setPdbFileName(s);
  if(DB) printf("leaving Sequence constructor \n");
}


// Default constructor
Sequence::Sequence() {
  if(DB) printf("entering Sequence constructor \n");
  name = 0;
  length = 0;
  lengthWithoutGaps = 0;
  isSecondary = false;
  isGap = 0;
  isEndGap = 0;
  residues = 0;
  residuesWithoutGaps = 0;
  residuesIndexToResiduesWithoutGapsIndex = 0;
  residuesWithoutGapsIndexToResiduesIndex = 0;
  if(DB) printf("leaving Sequence constructor \n");
}

// Destructor
Sequence::~Sequence() {
  if(DB) printf("entering Sequence destructor \n");
  if(isGap!=0) delete[] isGap;
  if(isEndGap!=0) delete[] isEndGap;
  if(name!=0) delete[] name;
  if(residues!=0) delete[] residues;
  if(residuesWithoutGaps!=0) delete[] residuesWithoutGaps;
  if(residuesIndexToResiduesWithoutGapsIndex!=0)
    delete[] residuesIndexToResiduesWithoutGapsIndex; 
  if(residuesWithoutGapsIndexToResiduesIndex!=0)
    delete[] residuesWithoutGapsIndexToResiduesIndex; 
  if(DB) printf("leaving Sequence destructor \n");
}

// copy constructor
Sequence::Sequence(const Sequence& sequence) {
  if(DB) printf("entering Sequence copy constructor \n");
  isSecondary = sequence.isSecondary;
  length = sequence.length;
  lengthWithoutGaps = sequence.lengthWithoutGaps;
  structure = sequence.structure;
  if(sequence.isGap==0) 
    isGap=0;
  else {
    isGap = new bool[length];
    for(int i=0;i<length;i++)
      isGap[i] = sequence.isGap[i];
  }
  if(sequence.isEndGap==0) 
    isEndGap=0;
  else {
    isEndGap = new bool[length];
    for(int i=0;i<length;i++)
      isEndGap[i] = sequence.isEndGap[i];
  }
  if(sequence.name==0) 
    name=0;
  else {
    name = new char[strlen(sequence.name)+1];
    strcpy(name,sequence.name);
  }
  if(sequence.residues==0) 
    residues=0;
  else {
    residues = new char[sequence.length+1];
    strcpy(residues,sequence.residues);
  }
  if(sequence.residuesWithoutGaps==0) 
    residuesWithoutGaps=0;
  else {
    residuesWithoutGaps = new char[sequence.lengthWithoutGaps+1];
    strcpy(residuesWithoutGaps,sequence.residuesWithoutGaps);
  }
  if(residuesIndexToResiduesWithoutGapsIndex==0)
    residuesIndexToResiduesWithoutGapsIndex=0;
  else {
    residuesIndexToResiduesWithoutGapsIndex = new int[length];
    for(int i=0;i<length;i++)
      residuesIndexToResiduesWithoutGapsIndex[i] = 
        sequence.residuesIndexToResiduesWithoutGapsIndex[i];
  }
  if(residuesWithoutGapsIndexToResiduesIndex==0)
    residuesWithoutGapsIndexToResiduesIndex=0;
  else {
    residuesWithoutGapsIndexToResiduesIndex = new int[lengthWithoutGaps];
    for(int i=0;i<lengthWithoutGaps;i++)
      residuesWithoutGapsIndexToResiduesIndex[i] = 
        sequence.residuesWithoutGapsIndexToResiduesIndex[i];
  }




  if(DB) printf("leaving Sequence copy constructor \n");
}

// Assignment
Sequence& Sequence::operator= (const Sequence& sequence) {
  if(DB) printf("entering Sequence::operator= \n");
  if(this != &sequence) {
    structure = sequence.structure;
    isSecondary = sequence.isSecondary;
    length = sequence.length;
    lengthWithoutGaps = sequence.lengthWithoutGaps;
    if(sequence.isGap==0) 
      isGap=0;
    else {
      isGap = new bool[length];
      for(int i=0;i<length;i++)
        isGap[i] = sequence.isGap[i];
    }
    if(sequence.isEndGap==0) 
      isEndGap=0;
    else {
      isEndGap = new bool[length];
      for(int i=0;i<length;i++)
        isEndGap[i] = sequence.isEndGap[i];
    }
    if(sequence.name==0) name=0;
    else {
      if(name!=0)
        delete[] name;
      name = new char[strlen(sequence.name)+1];
      strcpy(name,sequence.name);
    }
    if(sequence.residues==0) residues=0;
    else {
      if(residues!=0)
        delete[] residues;
      residues = new char[sequence.length+1];
      strcpy(residues,sequence.residues);
    }
    if(sequence.residuesWithoutGaps==0) residuesWithoutGaps=0;
    else {
      if(residuesWithoutGaps!=0) delete[] residuesWithoutGaps;
      residuesWithoutGaps = new char[sequence.lengthWithoutGaps+1];
      strcpy(residuesWithoutGaps,sequence.residuesWithoutGaps);
    }
    if(residuesIndexToResiduesWithoutGapsIndex==0)
      residuesIndexToResiduesWithoutGapsIndex=0;
    else {
      residuesIndexToResiduesWithoutGapsIndex = new int[length];
      for(int i=0;i<length;i++)
        residuesIndexToResiduesWithoutGapsIndex[i] = 
          sequence.residuesIndexToResiduesWithoutGapsIndex[i];
    }
    if(residuesWithoutGapsIndexToResiduesIndex==0)
      residuesWithoutGapsIndexToResiduesIndex=0;
    else {
      residuesWithoutGapsIndexToResiduesIndex = new int[lengthWithoutGaps];
      for(int i=0;i<lengthWithoutGaps;i++)
        residuesWithoutGapsIndexToResiduesIndex[i] = 
          sequence.residuesWithoutGapsIndexToResiduesIndex[i];
    }
  }
  if(DB) printf("leaving Sequence::operator= \n");
  return *this;
}


// Sets the name field
void Sequence::setName(char *s) {
  if(DB) printf("entering Sequence::setName\n");
  if(name!=0)
    delete[] name;
  name = new char[strlen(s)+1];
  strcpy(name,s);
  if(structure.pdbFileName==0) structure.setPdbFileName(s);
  if(DB) printf("leaving Sequence::setName\n");
}

// Sets the name field to the first n characters of s
void Sequence::setName(char *s, int n) {
  if(DB) printf("entering Sequence::setName\n");
  if(name!=0)
    delete[] name;
  name = new char[n+1];
  strncpy(name,s,n);
  name[n]='\0';
  if(structure.pdbFileName==0) structure.setPdbFileName(name);
  if(DB) printf("leaving Sequence::setName\n");
}

// Sets the residues field
void Sequence::setResidues(char *s) {
  if(DB) printf("entering Sequence::setResidues\n");
  if(residues!=0)
    delete[] residues;
  if(residuesWithoutGaps!=0)
    delete[] residuesWithoutGaps;

  residues = new char[strlen(s)+1];
  strcpy(residues,s);
  length = strlen(s);

  int count=0;
  for(int i=0;i<length;i++)
    if(residues[i]!='-') count++;
  lengthWithoutGaps = count;

  isGap = new bool[length];
  isEndGap = new bool[length];
  residuesWithoutGaps = new char[lengthWithoutGaps+1];
  residuesWithoutGapsIndexToResiduesIndex = new int[lengthWithoutGaps];
  residuesIndexToResiduesWithoutGapsIndex = new int[length];

  count=0;
  for(int i=0;i<length;i++) {
    if(residues[i]!='-') {
      isGap[i] = false;
      residuesWithoutGaps[count] = residues[i];
      residuesWithoutGapsIndexToResiduesIndex[count]=i;  
      residuesIndexToResiduesWithoutGapsIndex[i]=count;  
      count++;
    }
    else {
      isGap[i] = true;
      residuesIndexToResiduesWithoutGapsIndex[i]=-1;
    }
  }

  for(int i=0;         i<length; i++) isEndGap[i] = false;
  for(int i=0;         isGap[i]; i++) isEndGap[i] = true;
  for(int i=length-1;  isGap[i]; i--) isEndGap[i] = true;



//  residuesWithoutGaps[lengthWithoutGaps]='\0';
  
  if(DB) printf("leaving Sequence::setResidues\n");
}


// Returns the index in residuesWithoutGaps of the residue
// immediately preceding the gapped position pos
int Sequence::getLeft(int pos) {
  if(DB) printf("entering Sequence::getLeft\n");
  if(residues==0) return -1;
  int i;
  for(i=pos; i>-1 && isGap[i]==true; i--)
    ;
  if(DB) printf("leaving Sequence::getLeft\n");
  if(i==-1) return -1;
//  return residuesIndexToResiduesWithoutGapsIndex[i];
  return i;
}

// Returns the index in residuesWithoutGaps of the residue
// immediately following the gapped position pos
int Sequence::getRight(int pos) {
  if(DB) printf("entering Sequence::getRight\n");
  if(residues==0) return -1;
  int i;
  for(i=pos; i<length && isGap[i]==true; i++)
    ;
  if(DB) printf("leaving Sequence::getRight\n");
  if(i==length) return -1;
//  return residuesIndexToResiduesWithoutGapsIndex[i];
  return i;
}



// Initializes Structure
Structure::Structure(char *s, char *p) {
  if(DB) printf("entering Structure constructor \n");
  pdbFilePath = new char[strlen(p)+1];
  strcpy(pdbFilePath,p);
  pdbFileName = new char[strlen(s)+1];
  strcpy(pdbFileName,s);
  pdbChain = '_';  // default chain name equivalent to wildcard
  caCoordinates = 0;
  pdbNumbering = 0;
  if(DB) printf("leaving Structure constructor \n");
}

// Initializes Structure
Structure::Structure(char *s) {
  if(DB) printf("entering Structure constructor \n");
  pdbFilePath = new char[strlen(PDB_PATH)+1]; // assume default path
  strcpy(pdbFilePath,PDB_PATH);
  pdbFileName = new char[strlen(s)+1];
  strcpy(pdbFileName,s);
  pdbChain = '_';  // default chain name equivalent to wildcard
  length = 0;
  caCoordinates = 0;
  pdbNumbering = 0;
  if(DB) printf("leaving Structure constructor \n");
}

// Initializes Structure
Structure::Structure() {
  if(DB) printf("entering Structure constructor \n");
  pdbFilePath = new char[strlen(PDB_PATH)+1]; // assume default path
  strcpy(pdbFilePath,PDB_PATH);
  pdbFileName = 0; // set to 0 for now
  pdbChain = '_';  // default chain name equivalent to wildcard
  length = 0;
  caCoordinates = 0;
  pdbNumbering = 0;
  if(DB) printf("leaving Structure constructor \n");
}

// Destructor
Structure::~Structure() {
  if(DB) printf("entering Structure destructor \n");
  if(pdbFileName!=0) delete[] pdbFileName;
  if(pdbFilePath!=0) delete[] pdbFilePath;
  if(pdbNumbering!=0) delete[] pdbNumbering;
  if(caCoordinates!=0)
    for(int i=0; i<length; i++)
      delete[] caCoordinates[i];
  if(DB) printf("leaving Structure destructor \n");
}

// copy constructor
Structure::Structure(const Structure& structure) {
  if(DB) printf("entering Structure copy constructor \n");
  length = structure.length;
  pdbChain = structure.pdbChain;
  if(structure.pdbFileName==0)
    pdbFileName=0;
  else {
    pdbFileName = new char[strlen(structure.pdbFileName+1)];
    strcpy(pdbFileName,structure.pdbFileName);
  }
  if(structure.pdbFilePath==0)
    pdbFilePath=0;
  else {
    pdbFilePath = new char[strlen(structure.pdbFilePath+1)];
    strcpy(pdbFilePath,structure.pdbFilePath);
  }
  if(structure.caCoordinates==0) 
    caCoordinates=0;
  else {
    caCoordinates = new float*[length];
    for(int i=0; i<length; i++)
      caCoordinates[i] = new float[3];
    for(int i=0; i<length; i++)
      for(int j=0; j<3; j++)
        caCoordinates[i][j] = structure.caCoordinates[i][j];
  }
  if(structure.pdbNumbering==0)
    pdbNumbering=0;
  else {
    pdbNumbering = new int[length];
    for(int i=0; i<length; i++)
      pdbNumbering[i] = structure.pdbNumbering[i];
  }
  if(DB) printf("leaving Structure copy constructor \n");
}

// Assignment
Structure& Structure::operator= (const Structure& structure) {
  if(DB) printf("entering Structure::operator= \n");
  if(this != &structure) {
    length = structure.length;
    pdbChain = structure.pdbChain;
    if(structure.pdbFileName==0) pdbFileName=0;
    else {
      pdbFileName = new char[strlen(structure.pdbFileName)+1];
      strcpy(pdbFileName,structure.pdbFileName);
    }
    if(structure.pdbFilePath==0)
      pdbFilePath=0;
    else {
      pdbFilePath = new char[strlen(structure.pdbFilePath)+1];
      strcpy(pdbFilePath,structure.pdbFilePath);
    }
    if(structure.caCoordinates==0) 
      caCoordinates=0;
    else {
      caCoordinates = new float*[length];
      for(int i=0; i<length; i++)
        caCoordinates[i] = new float[3];
      for(int i=0; i<length; i++)
        for(int j=0; j<3; j++)
          caCoordinates[i][j] = structure.caCoordinates[i][j];
    }
    if(structure.pdbNumbering==0)
      pdbNumbering=0;
    else {
      pdbNumbering = new int[length];
      for(int i=0; i<length; i++)
        pdbNumbering[i] = structure.pdbNumbering[i];
    }
  }
  if(DB) printf("leaving Structure::operator= \n");
  return *this;
}



// Sets pdbFileName
void Structure::setPdbFileName(char *s) {
  if(DB) printf("entering Structure::setPdbFileName\n");
  if(pdbFileName!=0) delete[] pdbFileName;
  pdbFileName = new char[strlen(s)+1];
  strcpy(pdbFileName,s);
  if(DB) printf("leaving Structure::setPdbFileName\n");
}


// Sets pdbFilePath
void Structure::setPdbFilePath(char *s) {
  if(DB) printf("entering Structure::setPdbFilePath\n");
  if(pdbFilePath!=0) delete[] pdbFilePath;
  pdbFilePath = new char[strlen(s)+1];
  strcpy(pdbFilePath,s);
  if(DB) printf("leaving Structure::setPdbFilePath\n");
}


// Read and parse pdb file
void Structure::readFromPdbFile() {
  int i,j,num,numL=0; float r;
  // int k;
  char *temp = new char[1000];
  char *tempL = new char[1000];
  char *temp2 = new char[100];
  char *tName = new char[100];
  // Write warnings/errors to this file
  FILE *fpError = fopen("pdbError","a");

  float **tempCaCoordinates = new float*[MSL];
  for(i=0;i<MSL;i++)
    tempCaCoordinates[i] = new float[3];

  // Open file and get initial info
  strcpy(tName,pdbFilePath);
  strcat(tName,pdbFileName); 
  strcat(tName,".ent");
  FILE *fp = fopen(tName,"r");
  if(fp==0) {
    strcpy(tName,pdbFilePath);
    strcat(tName,pdbFileName); 
    strcat(tName,".pdb");
    fp = fopen(tName,"r");
  }
  if(fp==0) {
    printf("Error opening pdb file %s, with either .ent or .pdb extensions\n",tName); 
    return;}
  printf("tName = %s\n",tName);
  i=0;

  // Read through the file line by line
  while(fgets(temp,1000,fp)!=0) {

    // If it is an NMR structure, warn
    if(strstr(temp,"ENDMDL")!=0) {
      fprintf(fpError,"WARNING: In file %s, protein %s is an NMR structure\n",
          tName,pdbFileName);
      break;
    }

    // Only care about C-Alpha Atoms
    if(strstr(temp,"ATOM  ")!=0 && strstr(temp,"CA")!=0) {
      strncpy(temp2,temp+23,3); temp2[3]='\0'; num=charToInt(temp2);
      //        if(i==0) seqOffset[iSeq]=num;
      strncpy(temp2,temp+30,8); temp2[8]='\0';  //X
      tempCaCoordinates[i][0]=charToFloat(temp2);
      strncpy(temp2,temp+38,8); temp2[8]='\0';  //Y
      tempCaCoordinates[i][1]=charToFloat(temp2);
      strncpy(temp2,temp+46,8); temp2[8]='\0'; //Z
      tempCaCoordinates[i][2]=charToFloat(temp2);

      // If the numbering is off, warn
      if(i>0 && num-numL!=1) {
        r = sqrt( pow(tempCaCoordinates[i][0] - tempCaCoordinates[i-1][0],2)
                + pow(tempCaCoordinates[i][1] - tempCaCoordinates[i-1][1],2)
                + pow(tempCaCoordinates[i][2] - tempCaCoordinates[i-1][2],2) );
        fprintf(fpError,"WARNING: In file %s, CA neighbors %d and %d\n",tName,numL,num);
        if(r>3.000) {
          fprintf(fpError,"CA distance = %1.3f > 3.000, so assuming consecutive residues\n",r);
        } else {
          fprintf(fpError,"CA distance = %1.3f < 3.000, so assuming duplicate residues\n",r);
          i--;
        }
        fprintf(fpError,"%s%s",tempL,temp);
      }
      i++; numL=num; strcpy(tempL,temp);
    }
    length = i;
  }
  fclose(fp);

  // now reallocate memory correctly; wasteful - should be fixed
  caCoordinates = new float*[length];
  for(i=0;i<length;i++)
    caCoordinates[i] = new float[3];
  for(i=0;i<length;i++)
    for(j=0;j<3;j++)
      caCoordinates[i][j] = tempCaCoordinates[i][j];
  for(i=0;i<MSL;i++)
    delete[] tempCaCoordinates[i];

  delete[] temp;
  delete[] tempL;
  delete[] temp2;
  delete[] tName;

  return;
}*/