src/common/prospect_variables.cc

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <assert.h>


#include <libxml/tree.h>
#include <libxml/parser.h>
#include <libxml/xpath.h>
#include <libxml/xpathInternals.h>

#include "open_prospect.h"
#include "sm_matrix.h"



void get_cores(TemplateStruct *template_data);
void get_interactions(TemplateStruct *template_data);


//---template data



xmlNodePtr get_first_child(xmlNodePtr parent, char *str) {
  xmlNodePtr curnode;
  
  curnode = parent->children;
  while (curnode != NULL) {
    if (!strcmp((char *)curnode->name, str)) {
      return curnode;
    }
    curnode = curnode->next;
  }
  return NULL;
}


char *get_first_child_content( xmlNodePtr parent, char *str) {
        xmlNodePtr node = get_first_child(parent, str);
        if ( node == NULL )
                return NULL;
        return (char *)xmlNodeGetContent( node );
}



TemplateStruct::TemplateStruct() : TargetStruct() {
        this->len = 0;
        this->cb_dist = NULL;
        this->template_name = NULL;
        this->residue = NULL;
        this->cores = NULL;
        this->core_count = 0;
        this->t_residue_info = NULL;
        this->twobody_interactions = NULL;
}


TemplateStruct::~TemplateStruct() {
        FreeMem();
}


void TemplateStruct::SetLen(int new_len) {
        if ( new_len != len ) {
                TargetStruct::SetLen( new_len );        
                t_residue_info = new TemplateResidue[ new_len ];        
        }
}


void TemplateStruct::FreeMem(void) {
        TargetStruct::FreeMem();
        if (this->template_name != NULL)
                free(this->template_name);
        if (this->cb_dist) {
                free(this->cb_dist[0]);
                free(this->cb_dist);
        }
        this->len = 0;
        this->cb_dist = NULL;
        this->template_name = NULL;
        this->residue = NULL;
}

void TemplateStruct::CopyScoreTable( pValType **scores ) {
        for (int i = 0; i < len; i++ ) {
                int rs = AA1toAANum( residue[i].RS );
                for ( int j = 0; j < 20; j++ ) {
                        t_residue_info[i].Score[j] = scores[ rs ][ j ];
                }               
        }       
}

void TemplateStruct::CopyScoreTable( pValType scores[21][21] ) {
        for (int i = 0; i < len; i++ ) {
                int rs = AA1toAANum( residue[i].RS );
                for ( int j = 0; j < 20; j++ ) {
                        t_residue_info[i].Score[j] = scores[ rs ][ j ];
                }               
        }       
}



int TemplateStruct::SetSequence(char *seq) {
        TargetStruct::SetSequence( seq );
        t_residue_info = (TemplateResidue *)malloc(sizeof(TemplateResidue) * len );
}

int TemplateStruct::LoadFile( char *template_path ) {
  xmlDocPtr doc;
  xmlNodePtr root, fssp_node, profile_node, pdbcode_node, pdbchain_node, resSeq, iCode, 
    template_loc, template_start, template_stop;
  int i, j;
  float N[3], Ca[3], C[3], O[3], Cb[3], Cm[3];
  char *read_buffer, *cur_line;
  pValType tmp_val;

 
//  fprintf(stderr, "Loading %s\n", template_path);
  doc = xmlParseFile(template_path);
  if (doc == NULL)
    return 1;

  root = xmlDocGetRootElement(doc);
  if (root == NULL)
    return 1;
  //Load in the residue information from the file

  
  this->threebody_count = 0;
  this->threebody_interactions = NULL;
  this->fourbody_count = 0;
  this->fourbody_interactions = NULL;
          
  
  fssp_node = get_first_child(root, "dsspData");
  read_buffer = (char *)xmlNodeGetContent(fssp_node);
  cur_line = read_buffer;

  TargetResidue *tmp_residue = new TargetResidue[5000];
  TemplateResidue *tmp_residue_info = new TemplateResidue[5000];
  ResidueAtoms *tmp_structure = new ResidueAtoms[5000];
  
  i=0;
  while (*cur_line != 0) {
    if (!strncmp(cur_line,"RES", 3)) {
                /*
      sscanf(cur_line ,"RES %d%c%c%d%c%c%c%d%f%f%f%f%f%f%f%f%f%f%f",
             &tmp_residue[i].F,
             &dummyc,
             &tmp_residue[i].RS,
             &tmp_residue[i].NUM,
             &dummyc,&dummyc,
             &tmp_residue[i].SS,
             &tmp_residue[i].ACC,
             &tmp_residue[i].xCa,
             &tmp_residue[i].yCa,
             &tmp_residue[i].zCa,
             &tmp_residue[i].xCb,
             &tmp_residue[i].yCb,
             &tmp_residue[i].zCb,
             &tmp_residue[i].xCm,
             &tmp_residue[i].yCm,
             &tmp_residue[i].zCm,
             &tmp_residue[i].phi,
             &tmp_residue[i].psi);
*/
                if (24 != sscanf(cur_line ,"RES %c %d %c %d %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f",
                           &tmp_residue[i].RS,
                           &tmp_residue_info[i].NUM,
                           &tmp_residue[i].SS,
                           &tmp_residue_info[i].ACC,
                                                 &N[0],  &N[1],  &N[2], 
                                                 &Ca[0],  &Ca[1],  &Ca[2], 
                                                 &C[0],  &C[1],  &C[2], 
                                                 &O[0],  &O[1],  &O[2], 
                                                 &Cb[0],  &Cb[1],  &Cb[2], 
                                                 &Cm[0],  &Cm[1],  &Cm[2],
                                                 &tmp_residue[i].phi,
                                                 &tmp_residue[i].psi) ) {
                        fprintf(stderr, "Template Formating error: \"%s\"", cur_line);
                        delete tmp_residue;
                        delete tmp_residue_info;
                        delete tmp_structure;
                        return 1;
                } else {
                        for (int j = 0; j < 3; j++) {
                                tmp_structure[i].atom[anN][j] = N[j];
                                tmp_structure[i].atom[anCa][j] = Ca[j];
                                tmp_structure[i].atom[anC][j] = C[j];
                                tmp_structure[i].atom[anO][j] = O[j];
                                tmp_structure[i].atom[anCb][j] = Cb[j];
                                tmp_structure[i].Cm[j] = Cm[j];                         
                        }
                        tmp_structure[i].atom_name[anN] = anN;
                        tmp_structure[i].atom_name[anCa] = anCa;
                        tmp_structure[i].atom_name[anC] = anC;
                        tmp_structure[i].atom_name[anO] = anO;
                        tmp_structure[i].atom_name[anCb] = anCb;

                        tmp_structure[i].bb_num[anN] = anN;
                        tmp_structure[i].bb_num[anCa] = anCa;
                        tmp_structure[i].bb_num[anC] = anC;
                        tmp_structure[i].bb_num[anO] = anO;
                        tmp_structure[i].bb_num[anCb] = anCb;
                        
                        //flag 1 if you want to  get the  CB from the template file
#if 0
                        if ( tmp_residue[i].RS != 'G' ) {
                                strcpy(tmp_residue[i].structure.name[anCb], "CB");
                                tmp_structure[i].bb_count = 4;
                                tmp_structure[i].atom_count = 5;
                        } else {
                                tmp_structur[i].bb_count = 4;
                                tmp_structure[i].atom_count = 4;
                        }
#else 
                        tmp_structure[i].atom_count = 4;
#endif
                        tmp_structure[i].chain = ' ';
                        tmp_structure[i].res_mod = ' ';
                        
                //      printf("read: %c %d %c %f %f %f\n", tmp_residue[i].RS, tmp_residue[i].NUM, tmp_residue[i].SS, 
                //                 tmp_residue[i].xCa, tmp_residue[i].yCa, tmp_residue[i].zCa);
                }
                tmp_residue[i].single_del = 0.0;
                tmp_residue[i].single_ins = 0.0;
                
                switch (tmp_residue[i].SS) {
                        case 'H': tmp_residue[i].SS='H'; break;
                        case 'G': tmp_residue[i].SS='H'; break;
                        case 'I': tmp_residue[i].SS='H'; break;
                        case 'E': tmp_residue[i].SS='E'; break;
                        case 'B': tmp_residue[i].SS='E'; break;
                        default : tmp_residue[i].SS='C'; break;
                }
                i++;
    }
          //advance pointer to next line or end of file
    while (cur_line[0] != 0 && cur_line[0] != '\n') {
      cur_line++;
    }
    if (cur_line[0] == '\n')
      cur_line++;
  }
  this->len = i;
  free(read_buffer);
  
  if (this->len == 0) {
          delete tmp_residue;
          delete tmp_residue_info;
          delete tmp_structure;
          xmlFreeDoc(doc); 
          return 1;
  }

#ifdef ALTIVEC
  this->residue=(TemplateResidue *)valloc(this->len * sizeof(TemplateResidue));  
  assert(this->residue);
#else
  this->residue =(TargetResidue *)malloc(this->len * sizeof(TargetResidue));
#endif  
  this->t_residue_info = new TemplateResidue[ this->len ];
  this->structure = new ResidueAtoms[ this->len ];
  
  for (i = 0; i < this->len; i++) {
          this->residue[i] = tmp_residue[i]; 
          this->t_residue_info[i] = tmp_residue_info[i];
          this->structure[i] = tmp_structure[i];
  }
  delete tmp_residue;
  delete tmp_residue_info;
  delete tmp_structure;
  
  profile_node = get_first_child(root, "frequency_matrix");
  read_buffer = (char *)xmlNodeGetContent(profile_node);
  cur_line = read_buffer;
  while (cur_line[0] == '\n') {
      cur_line++;
  }

  for (i=0; i < this->len; i++) {
    sscanf(cur_line, "%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f",
           &this->residue[i].FREQ[0],
           &this->residue[i].FREQ[1],
           &this->residue[i].FREQ[2],
           &this->residue[i].FREQ[3],
           &this->residue[i].FREQ[4],
           &this->residue[i].FREQ[5],
           &this->residue[i].FREQ[6],
           &this->residue[i].FREQ[7],
           &this->residue[i].FREQ[8],
           &this->residue[i].FREQ[9],
           &this->residue[i].FREQ[10],
           &this->residue[i].FREQ[11],
           &this->residue[i].FREQ[12],
           &this->residue[i].FREQ[13],
           &this->residue[i].FREQ[14],
           &this->residue[i].FREQ[15],
           &this->residue[i].FREQ[16],
           &this->residue[i].FREQ[17],
           &this->residue[i].FREQ[18],
           &this->residue[i].FREQ[19]);
          this->residue[i].FREQ[20] = 0;
    while (cur_line[0] != '\n') {
      cur_line++;
    }
    if (cur_line[0] == '\n')
      cur_line++;
    //if we hit the end of the profile info early
    if (cur_line[0] == 0 && i < this->len-1) 
      return 1;
  }
  free(read_buffer);


  //do the log score read in 
  
  
  profile_node = get_first_child(root, "score_matrix");
  read_buffer = (char *)xmlNodeGetContent(profile_node);
  cur_line = read_buffer;
  while (cur_line[0] == '\n') {
      cur_line++;
  }
  
  for (i=0; i < this->len; i++) {
          long score[20];
          sscanf(cur_line, "%ld%ld%ld%ld%ld%ld%ld%ld%ld%ld%ld%ld%ld%ld%ld%ld%ld%ld%ld%ld",
                         &score[0],
                         &score[1],
                         &score[2],
                         &score[3],
                         &score[4],
                         &score[5],
                         &score[6],
                         &score[7],
                         &score[8],
                         &score[9],
                         &score[10],
                         &score[11],
                         &score[12],
                         &score[13],
                         &score[14],
                         &score[15],
                         &score[16],
                         &score[17],
                         &score[18],
                         &score[19]);
          
          for (int j = 0; j < 20; j++) 
                  this->t_residue_info[i].Score[j] = score[j];
          this->t_residue_info[i].Score[20] = 0;
          while (cur_line[0] != '\n') {
                  cur_line++;
          }
          if (cur_line[0] == '\n')
                  cur_line++;
          //if we hit the end of the profile info early
          if (cur_line[0] == 0 && i < this->len-1) 
                  return 1;
  }
  free(read_buffer);
  
  
  
  char *single_del_text = get_first_child_content(root, "singledel");
  if (single_del_text) {
          ReadSingleDelStr( single_del_text );
          free(single_del_text);
  }
  char *single_ins_text = get_first_child_content(root, "singleins");
  if (single_ins_text) {
          ReadSingleInsStr( single_ins_text );
          free(single_ins_text);
  }
  
  
//-------------------------------------------------------
/*
  xmlNodePtr core_node;
  read_buffer = (char *)xmlNodeGetContent(core_node);
  cur_line = read_buffer;
  sscanf(cur_line, "%s%c%c%d",dummy,&dummyc,&dummyc,&this->core_count);
  printf("%d\n",this->core_count);
  free(read_buffer);
*/
//-------------------------------------------------------
  
  //get PDB code
  pdbcode_node = get_first_child(root, "pdbcode");
  read_buffer = (char *)xmlNodeGetContent(pdbcode_node);
  for (i = 0; i < 4; i++) {
    this->pdbcode[i] = read_buffer[i];
  }
  this->pdbcode[4] = 0;
  free(read_buffer);

  //pdb chain
  pdbchain_node = get_first_child(root, "pdbchain");
  read_buffer = (char *)xmlNodeGetContent(pdbchain_node);
  this->pdbchain = read_buffer[0];
  free(read_buffer);

  //template name
  this->template_name = (char *)xmlGetProp(root, (xmlChar *)"name");

  //template_loc
  template_loc = get_first_child(root, "template_loc");
  template_start = get_first_child( template_loc, "start");
  template_stop = get_first_child( template_loc, "stop");

  resSeq = get_first_child(template_start, "resSeq");
  iCode = get_first_child(template_start, "iCode");
  read_buffer = (char *)xmlNodeGetContent(resSeq);
  this->start_res = atoi( read_buffer );
  free(read_buffer);
  read_buffer = (char *)xmlNodeGetContent(iCode);
  if (read_buffer[0] == 0)
    read_buffer[0] = ' ';
  this->start_icode = read_buffer[0];
  free(read_buffer);

  resSeq = get_first_child(template_stop, "resSeq");
  iCode = get_first_child(template_stop, "iCode");
  read_buffer = (char *)xmlNodeGetContent(resSeq);
  this->stop_res = atoi( read_buffer );
  free(read_buffer);
  read_buffer = (char *)xmlNodeGetContent(iCode);
  if (read_buffer[0] == 0)
    read_buffer[0] = ' ';
  this->stop_icode = read_buffer[0];
  free(read_buffer);


  xmlFreeDoc(doc); 

  //  Determine template cores
  this->cores=(TemplateCoreDefinition *)malloc(1000*sizeof(TemplateCoreDefinition));
  get_interactions(this);
  get_cores(this);


  //Calculate CB distances between all residues

  this->cb_dist = (pValType **)malloc(sizeof(pValType *) * this->len);
  this->cb_dist[0] = (pValType *)malloc(sizeof(pValType) * this->len * this->len);
  for (i = 0; i < this->len; i++) {
    this->cb_dist[i] = this->cb_dist[0] + i * this->len;
  }


  for (i = 0; i < this->len; i++) {
          int i_cb = this->structure[i].bb_num[ anCb ];
          for (j = i; j < this->len; j++) {
                  int j_cb = this->structure[j].bb_num[ anCb ];
                  if ( i_cb != -1 && j_cb != -1 ) {
                          tmp_val = sqrt( pow(this->structure[i].atom[ i_cb ][0] - 
                                                                  this->structure[j].atom[ j_cb ][0], 2) +
                                                          pow(this->structure[i].atom[ i_cb ][1] - 
                                                                  this->structure[j].atom[ j_cb ][1], 2) +
                                                          pow(this->structure[i].atom[ i_cb ][2] - 
                                                                  this->structure[j].atom[ j_cb ][2], 2) );      
                  } else {
                          tmp_val = 10000;
                  }
                  this->cb_dist[i][j] = tmp_val;
                  this->cb_dist[j][i] = tmp_val;
          }
  }

  return 0;
}




//--------------------------------------------------------------------
void get_interactions(TemplateStruct *template_data) {
        int i,j,cnt;
        float x,y,z,distance;
        cnt = 0;
        
        for (i=0;i<template_data->len;i++) {
                for (j=i+MinDistTwoBody;j<template_data->len;j++) {                     
                        x = template_data->structure[i].atom[anCb][0] - template_data->structure[j].atom[anCb][0];
                        y = template_data->structure[i].atom[anCb][1] - template_data->structure[j].atom[anCb][1];
                        z = template_data->structure[i].atom[anCb][2] - template_data->structure[j].atom[anCb][2];
                        
                        distance = x*x + y*y + z*z;
                        
                        if (distance <= PairDistCutoff) {
                                cnt += 1;
                        }
                }
        }
        
        template_data->twobody_count = cnt;     
        template_data->twobody_interactions = (InteractionElement *)malloc(template_data->twobody_count*sizeof(InteractionElement));    
        cnt = 0;
        
        for (i=0;i<template_data->len;i++) {
                for (j=i+MinDistTwoBody;j<template_data->len;j++) {
                        
                        x = template_data->structure[i].atom[anCb][0] - template_data->structure[j].atom[anCb][0];
                        y = template_data->structure[i].atom[anCb][1] - template_data->structure[j].atom[anCb][1];
                        z = template_data->structure[i].atom[anCb][2] - template_data->structure[j].atom[anCb][2];                      
                        distance = x*x + y*y + z*z;                     
                        if (distance <= PairDistCutoff) {
                                template_data->twobody_interactions[cnt].member_count = 2;
                                template_data->twobody_interactions[cnt].residue_number[0] = i;
                                template_data->twobody_interactions[cnt].residue_number[1] = j;
                                cnt++;
                        }
                }
        }
}

//--------------------------------------------------------------------

//my $MinHelixLen = 5;   # old 4
//my $MinSheetLen = 2;   # old 3

#define cMinHelix  5
#define cMinStrand  3

//#define DEBUG 1

void get_cores(TemplateStruct *template_data) {
  char last_state='C';
  int cur_core=0;

  //first identify all possible cores
  for (int i=0; i < template_data->len; i++) {
    if (last_state != 'C' && last_state != template_data->residue[i].SS) {
      if ( cur_core > 0)
              template_data->cores[cur_core-1].len = i - template_data->cores[cur_core-1].head;
    }
    if ( (template_data->residue[i].SS == 'H' || template_data->residue[i].SS == 'E') && 
        last_state != template_data->residue[i].SS) {
      template_data->cores[cur_core].head=i;
      cur_core++;
    }
    last_state = template_data->residue[i].SS;
  }
  if (last_state != 'C') {
          template_data->cores[cur_core-1].len = template_data->len - template_data->cores[cur_core-1].head;
  }
  template_data->core_count=cur_core;

  
   
  //remove cores that are too short
  for (int i = 0; i < template_data->core_count; i++) {
#ifdef DEBUG
          fprintf(stderr, "Core %d (%c): %d\n", i, template_data->residue[ template_data->cores[i].head ].SS,  template_data->cores[i].len );
#endif
          if ( template_data->residue[ template_data->cores[i].head ].SS == 'H' && template_data->cores[i].len < cMinHelix ) {
                  for (int j = i; j < template_data->core_count-1; j++) {
                          template_data->cores[j] = template_data->cores[j+1];
                  }
                  template_data->core_count--;
                  i--;
          } else if ( template_data->residue[ template_data->cores[i].head ].SS == 'E' && template_data->cores[i].len < cMinStrand ) {
                  for (int j = i; j < template_data->core_count-1; j++) {
                          template_data->cores[j] = template_data->cores[j+1];
                  }
                  template_data->core_count--;
                  i--;
          }
  }
  
#if 1

  //remove cores with not enough contacts  
  int core_residue_count = 0, core_interaction_count = 0;
  for (int i = 0; i < template_data->core_count; i++) {
          core_residue_count += template_data->cores[i].len;
  }
  for (int i = 0; i < template_data->core_count; i++) {
          for (int j = 0; j < template_data->core_count; j++) {
                  if ( i != j) {
                          for (int k= 0; k < template_data->twobody_count; k++) {
                                  if ( template_data->twobody_interactions[k].residue_number[0] >= template_data->cores[i].head &&
                                           template_data->twobody_interactions[k].residue_number[0] <  template_data->cores[i].head + template_data->cores[i].len &&
                                           template_data->twobody_interactions[k].residue_number[1] >= template_data->cores[j].head &&
                                           template_data->twobody_interactions[k].residue_number[1] <  template_data->cores[j].head + template_data->cores[i].len ) {
                                          core_interaction_count++;                                       
                                  }
                          }
                  }
          }       
  }  
  float core_interaction_ratio = (float)core_interaction_count / (float) core_residue_count;  
#ifdef DEBUG 
        fprintf(stderr, "protein core interaction ratio: %f\n", core_interaction_ratio);
#endif

  for (int i = 0; i < template_data->core_count; i++) {
          long cur_interaction_count = 0;
#ifdef DEBUG
        fprintf(stderr, "core %d\n", i);
#endif
          for (int j = 0; j < template_data->core_count; j++) {
                  if ( i != j) {
                          for (int k= 0; k < template_data->twobody_count; k++) {
                                  if (( template_data->twobody_interactions[k].residue_number[0] >= template_data->cores[i].head &&
                                           template_data->twobody_interactions[k].residue_number[0] <  template_data->cores[i].head + template_data->cores[i].len &&
                                           template_data->twobody_interactions[k].residue_number[1] >= template_data->cores[j].head &&
                                           template_data->twobody_interactions[k].residue_number[1] <  template_data->cores[j].head + template_data->cores[i].len) ||
                                        ( template_data->twobody_interactions[k].residue_number[1] >= template_data->cores[i].head &&
                                           template_data->twobody_interactions[k].residue_number[1] <  template_data->cores[i].head + template_data->cores[i].len &&
                                           template_data->twobody_interactions[k].residue_number[0] >= template_data->cores[j].head &&
                                           template_data->twobody_interactions[k].residue_number[0] <  template_data->cores[j].head + template_data->cores[i].len)

 ) {
                                          cur_interaction_count++;      
#ifdef DEBUG
                                          fprintf(stderr, "\t%d -- %d\n", template_data->twobody_interactions[k].residue_number[0],
                                                                template_data->twobody_interactions[k].residue_number[1] );
#endif                            
                                  }
                          }
                  }
          }       
#ifdef DEBUG
        fprintf(stderr, "\tinteraction ratio %f\n", ( (float) cur_interaction_count / (float) template_data->cores[i].len ) );
#endif
          if ( 0.3 * core_interaction_ratio > ( (float) cur_interaction_count / (float) template_data->cores[i].len ) ) {
                  //removing the core
                  for (int j = i; j < template_data->core_count-1; j++) {
                          template_data->cores[j] = template_data->cores[j+1];
                  }
                  template_data->core_count--;          
                  i--;
          }       
  }
  

  //remove terminal residues that don't have enough contacts  
  for (int i = 0; i < template_data->core_count; i++) {
          long cur_interaction_count = 0;
          long residue_interaction_count[ template_data->cores[i].len ];
          for (int j = 0; j < template_data->cores[i].len; j++) {
                  residue_interaction_count[j] = 0;
          }
          for (int j = 0; j < template_data->core_count; j++) {
            if ( i != j) {
              for (int k= 0; k < template_data->twobody_count; k++) {
                if ( template_data->twobody_interactions[k].residue_number[0] >= template_data->cores[i].head &&
                     template_data->twobody_interactions[k].residue_number[0] <  template_data->cores[i].head + template_data->cores[i].len &&
                     template_data->twobody_interactions[k].residue_number[1] >= template_data->cores[j].head &&
                     template_data->twobody_interactions[k].residue_number[1] <  template_data->cores[j].head + template_data->cores[j].len 
                     ) {
                  cur_interaction_count++;
                  residue_interaction_count[ template_data->twobody_interactions[k].residue_number[0] - template_data->cores[i].head ]++;                                         
                } else if (template_data->twobody_interactions[k].residue_number[1] >= template_data->cores[i].head &&
                           template_data->twobody_interactions[k].residue_number[1] <  template_data->cores[i].head + template_data->cores[i].len &&
                           template_data->twobody_interactions[k].residue_number[0] >= template_data->cores[j].head &&
                           template_data->twobody_interactions[k].residue_number[0] <  template_data->cores[j].head + template_data->cores[j].len
                           ) {
                  cur_interaction_count++;
                  residue_interaction_count[ template_data->twobody_interactions[k].residue_number[1] - template_data->cores[i].head ]++;
                }
              }
            }
          }

#ifdef DEBUG
          
          for (int k = 0; k < template_data->cores[i].len; k++) {
            fprintf(stderr, "%d ", residue_interaction_count[k]);
          }
#endif
          char removed = 0;
          long total_removed = 0, new_head = 0, new_len = template_data->cores[i].len;  
          do {
                  removed = 0;
                  if ( residue_interaction_count[ new_head ] < residue_interaction_count[ new_len-1 ] ) {
                          if (  total_removed + residue_interaction_count[new_head] < 0.1 * (float) cur_interaction_count ) {
                                  total_removed += residue_interaction_count[new_head];
                                  new_head++;
                                  removed = 1;
                          }
                  } else {
                          if (  total_removed + residue_interaction_count[new_len-1] < 0.1 * (float) cur_interaction_count ) {
                                  total_removed += residue_interaction_count[new_len-1];
                                  new_len--;
                                  removed = 1;
                          }
                  }
          } while (removed);    
#ifdef DEBUG 
          fprintf(stderr, "Core %d (%c)  head %d (%d)  len: %d (%d )\n", i, template_data->residue[ template_data->cores[i].head ].SS,
                         template_data->cores[i].head + new_head, template_data->cores[i].head,
                         new_len - new_head,  template_data->cores[i].len );
#endif  
          template_data->cores[i].head += new_head;
          template_data->cores[i].len = new_len - new_head;  
          
  }
  
  //again, remove cores that are too short
  for (int i = 0; i < template_data->core_count; i++) {
          if ( template_data->residue[ template_data->cores[i].head ].SS == 'H' && template_data->cores[i].len < cMinHelix-1 ) {
                  for (int j = i; j < template_data->core_count-1; j++) {
                          template_data->cores[j] = template_data->cores[j+1];
                  }
                  template_data->core_count--;
                  i--;
          } else if ( template_data->residue[ template_data->cores[i].head ].SS == 'E' && template_data->cores[i].len < cMinStrand-1 ) {
                  for (int j = i; j < template_data->core_count-1; j++) {
                          template_data->cores[j] = template_data->cores[j+1];
                  }
                  template_data->core_count--;
                  i--;
          }
  }
  
#endif  

  
  
  
  //label residues by their core numbers
  for (int i = 0; i < template_data->len; i++) {
    template_data->t_residue_info[i].core_num = -1;
  }
  for (int i = 0; i < template_data->core_count; i++) {
          for (int j = 0; j < template_data->cores[i].len; j++) {
                  template_data->t_residue_info[ template_data->cores[i].head + j ].core_num = i;
          }
  }  
  
  //go through and label interactions by core #  
  for (int i = 0; i < template_data->core_count; i++) {
          for (int j = 0; j < template_data->twobody_count; j++) {
                  for (int k = 0; k < 2; k++) {
                          if ( template_data->twobody_interactions[j].residue_number[k] >= template_data->cores[i].head &&
                                   template_data->twobody_interactions[j].residue_number[k] < template_data->cores[i].head + template_data->cores[i].len ) {
                                  template_data->twobody_interactions[j].core_number[k] = i;
                          }
                  }
          }
  }
}

/*

void ProteinLinearAtoms::extract_residueatom_array( ResidueAtoms *r_atoms, int count) { 
        int counter = 0;
        for (int i= 0; i < count; i++) {
                counter += r_atoms[i].atom_count;
        }       
        atom_count = counter;
        atom_loc = (vec3 *)malloc(sizeof(vec3) * atom_count );
        residue_num = (int *)malloc(sizeof(int) * atom_count );
        atom_type = (int *)malloc(sizeof(int) * atom_count );
        temp_factor = (float *)malloc(sizeof(float) * atom_count );     
        counter = 0;
        for (int i= 0; i < count; i++) {
                for (int j = 0; j < r_atoms[i].atom_count; j++) {
                        atom_loc[ counter ][0] = r_atoms[i].atom[j][0];
                        atom_loc[ counter ][1] = r_atoms[i].atom[j][1];
                        atom_loc[ counter ][2] = r_atoms[i].atom[j][2];
                        residue_num[ counter ] = i;
                        temp_factor[ counter ] = r_atoms[i].temp_factor[j];
                        counter++;
                }
        }       
}

*/

AlignFeatures::AlignFeatures(void) {
        for ( int i =0; i < f_count; i++ ) {
                is_set[i] = 0;
        }
}


void AlignFeatures::CalcFeatures( class TemplateStruct *template_data, 
                                                                  class TargetStruct *target_data,  
                                                                  class AlignmentStruct *align ) {      
        
        
        //alignment stuff
        feature[ fn_align ] = 0;
        feature[ fn_ident ] = 0;
        for (int i = 0; i < target_data->len; i++) {
                if (align->GetTargetAlign( i )!= -1 ) {
                        feature[ fn_align ]++;
                        if ( target_data->residue[i].RS == template_data->residue[ align->GetTargetAlign( i ) ].RS ) {
                                feature[fn_ident]++;
                        }
                }
        }
        feature[ fr_align ] = feature[fn_align] / (float)target_data->len;
        feature[ fr_ident ] = feature[fn_ident] / (float)target_data->len;      
        is_set[ fr_align ] = 1;
        is_set[ fn_align ] = 1;
        is_set[ fr_ident ] = 1;
        is_set[ fn_ident ] = 1;
                
        //get number of aligned residue pairs
        feature[ fn_align_pair ] = 0;
        long pair_total = 0;
        for (int i = 0; i < template_data->twobody_count; i++) {
                if ( template_data->twobody_interactions[i].core_number[0] != -1 &&
                         template_data->twobody_interactions[i].core_number[1] != -1 ) {
                        pair_total++;
                        if (align->GetTemplateAlign( template_data->twobody_interactions[i].residue_number[0] ) != -1 &&
                                align->GetTemplateAlign( template_data->twobody_interactions[i].residue_number[1] ) != -1 ) {
                                feature[ fn_align_pair ]++;
                        }
                }
        }
        feature[ fr_align_pair ] = feature[ fn_align_pair ] / (double) pair_total;
        is_set[ fr_align_pair ] = 1;
        is_set[ fn_align_pair ] = 1;
        
        //core alignment stuff
        //BUG: this assumes the whole core is aligned...
        feature[ fn_align_core_res ] = 0;
        feature[ fn_align_core ] = 0;
        for (int j = 0; j < template_data->core_count; j++) {
                if ( align->GetTemplateAlign( template_data->cores[j].head ) != -1 ) {
                        feature[ fn_align_core ]++;
                }               
        }
        for (int j = 0; j < template_data->core_count; j++) {
                feature[ fn_align_core_res ] += template_data->cores[j].len;            
        }
        is_set[ fn_align_core_res ] = 1;
        is_set[ fn_align_core ] = 1;
        
        
        
        
        //contact stuff
        feature[ fn_contact ] = 0;
        feature[ fn_halfcontact ] = 0;
        int total_contact = 0;
        for ( int i = 0; i < template_data->len; i++ ) {
                ResidueAtoms *res_1 = &template_data->structure[i];
                int align_1 = align->GetTemplateAlign(i);
                for ( int j = i+1; j < template_data->len; j++ ) {
                        ResidueAtoms *res_2 = &template_data->structure[j];
                        int align_2 = align->GetTemplateAlign(j);
                        float dist = smDistance( res_1->atom[ res_1->bb_num[ anCa ] ], 
                                                                         res_2->atom[ res_2->bb_num[ anCa ] ] );
                        if ( dist < 7.6 ) {
                                total_contact++;
                                if ( (align_1 == -1 && align_2 != -1) || (align_1 != -1 && align_2 == -1) ) {
                                        feature[ fn_halfcontact ]++;
                                }
                                if ( (align_1 != -1 && align_2 != -1)  ) {
                                        feature[ fn_contact ]++;
                                }                                               
                        }
                }
        }       
        feature[ fr_contact ] = feature[ fn_contact ] / (float) total_contact;
        feature[ fr_halfcontact ] = feature[ fn_halfcontact ] / (float) total_contact;
        is_set[ fr_contact ] = 1;
        is_set[ fn_contact ] = 1;
        is_set[ fr_halfcontact ] = 1;
        is_set[ fn_halfcontact ] = 1;   
}

---