src/common/prospect_output.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"


extern "C" {
#include "foldrec_net.h"
}







//private functions defined in prospect_variables.cc
xmlNodePtr get_first_child(xmlNodePtr parent, char *str);
char *get_first_child_content( xmlNodePtr parent, char *str);

//Functions relateing to Prospect Output files




ProspectOutput::ProspectOutput(char *path) {
        reorder = NULL;
        sort_val = NULL;
        output_xml_handle = NULL;
        if (path != NULL)
                Load(path);
}


ProspectOutput::~ProspectOutput() {
        Free();
}


void ProspectOutput::Free(void){
        if (output_xml_handle != NULL) 
                xmlFreeDoc( (xmlDocPtr) output_xml_handle );
        output_xml_handle = NULL;
        if (reorder != NULL)
                free(reorder);
        reorder = NULL;
        if ( sort_val )
                delete sort_val;
        sort_val = NULL;
}



int ProspectOutput::Load( char *path) {
        output_xml_handle = xmlParseFile( path );
        if (output_xml_handle == NULL) {
                //fprintf(stderr, "Can't open %s\n", argv[1]);
                return 1;
        }       
        //cur_node = xmlDocGetRootElement(doc); 
        return 0;
}


int ProspectOutput::AddThreadingInfo( ProspectThreadingInfo *thread_info ) {    
        if ( output_xml_handle == NULL ) {
                output_xml_handle = (void *)xmlNewDoc( (xmlChar *)"1.0" );
                xmlDocSetRootElement( (xmlDoc *)output_xml_handle, xmlNewNode( NULL, (xmlChar *)"prospectOutput" ));
        }       
        if ( thread_info->xml_node == NULL ) {
                return 1;
        }       
        xmlNodePtr new_node = xmlCopyNode( (xmlNode *)thread_info->xml_node, 1);
        xmlNodePtr root_node = xmlDocGetRootElement( (xmlDoc *)output_xml_handle );
        xmlAddChild(root_node, new_node);       
        return 0;
}



int ProspectOutput::AddThreadingInfo( ProspectParam *param,
                                          WeightArray weight,   
                                          TemplateStruct *template_data,
                                          TargetStruct *target_data, 
                                          ScoreStruct *score,
                                          AlignmentStruct *alignment) {

        if ( output_xml_handle == NULL ) {
                output_xml_handle = (void *)xmlNewDoc( (xmlChar *)"1.0" );
                xmlDocSetRootElement( (xmlDoc *)output_xml_handle, xmlNewNode( NULL, (xmlChar *)"prospectOutput" ));
        }       
        
        char *xml_text = alignment->Write_XML(param,
                                                                                 weight,        
                                                                                 template_data,
                                                                                 target_data, 
                                                                                 score);
        if ( xml_text != NULL ) {
                xmlNodePtr root_node = xmlDocGetRootElement( (xmlDoc *)output_xml_handle );
                xmlDocPtr tmp_doc = xmlParseMemory (xml_text, strlen( xml_text ) ); 
                xmlNodePtr new_node = xmlDocGetRootElement( tmp_doc );
                xmlUnlinkNode(new_node);
                xmlAddChild (root_node,new_node); 
                xmlFreeDoc( tmp_doc );  
                free( xml_text );
        }       
        return 0;
}




int ProspectOutput::Save( char *path ) {        
        return xmlSaveFile( path, (xmlDoc *)output_xml_handle); 
}



int ProspectOutput::GetThreadingCount(void ) {
        xmlDocPtr doc = (xmlDocPtr) output_xml_handle;
        if ( doc == NULL )
                return 0;
        xmlNodePtr cur_node = xmlDocGetRootElement(doc);
        if (cur_node == NULL)
                return 0;
        if ( !strcmp((char *)cur_node->name, "prospectOutput") ) {
                cur_node = cur_node->children;
        }
        if (cur_node == NULL)
                return 0;
        int count = 0;  
        while (cur_node != NULL) {
                if (!strcmp((char *)cur_node->name, "threading")) {
                        count++;
                }
                cur_node = cur_node->next;
        }       
        return count;   
}


char * ProspectOutput::GetTemplateName( int N) {
        xmlDocPtr doc = (xmlDocPtr) output_xml_handle;
        xmlNodePtr cur_node = xmlDocGetRootElement(doc);
        if ( !strcmp((char *)cur_node->name, "prospectOutput") ) {
                cur_node = cur_node->children;
        }
        int count = 0;  
        while (cur_node != NULL) {
                if (!strcmp((char *)cur_node->name, "threading")) {
                        if ( count == N ) {
                                return (char *)xmlGetProp(cur_node, (xmlChar *)"template");                             
                        }
                        count++;
                }
                cur_node = cur_node->next;
        }       
        return NULL;
}


int ProspectOutput::GetTemplateNum( char *name) {
        xmlDocPtr doc = (xmlDocPtr) output_xml_handle;
        xmlNodePtr cur_node = xmlDocGetRootElement(doc);
        if ( !strcmp((char *)cur_node->name, "prospectOutput") ) {
                cur_node = cur_node->children;
        }
        int count = 0, match_num = -1;  
        while (cur_node != NULL) {
                if (!strcmp((char *)cur_node->name, "threading")) {
                        char *cur_name = (char *)xmlGetProp(cur_node, (xmlChar *)"template");
                        if (cur_name) {
                                if ( !strcmp( cur_name, name ) ) {
                                        match_num = count;
                                        cur_node = NULL;
                                }
                                free(cur_name);
                        }
                        count++;
                }
                if ( cur_node != NULL)
                        cur_node = cur_node->next;
        }       
        return match_num;
}


void ProspectOutput::GetRawNormalize( float &mean, float &sd ) {
        xmlDocPtr doc = (xmlDocPtr) output_xml_handle;  
        xmlXPathContextPtr xpathCtx = xmlXPathNewContext(doc); /* Create xpath evaluation context */
        char searchpath[200] = "/prospectOutput/threading/score/raw";
        xmlXPathObjectPtr xpathObj = xmlXPathEvalExpression( (const xmlChar*)searchpath, 
                                                                                                                 xpathCtx);             
        if(xpathObj == NULL) {
                xmlXPathFreeContext(xpathCtx); 
                return;
        }       
        int size = xpathObj->nodesetval->nodeNr; 
        if ( size > 0 ) {       
          float *raw_array = new float[ size ]; 
          for(int i = 0; i < size; i++) {
            char * text = (char *)xmlNodeGetContent(xpathObj->nodesetval->nodeTab[i]);
            raw_array[i] = atof( text );
            free(text);
          }     
          /* Cleanup of XPath data */
          xmlXPathFreeObject(xpathObj);
          xmlXPathFreeContext(xpathCtx);        
          double total = 0.0;
          for (int i = 0; i < size; i++) {
            total += raw_array[i];
        }
          mean = total / (double)size;  
          total = 0.0;
          for (int i = 0; i < size; i++) {
            total += pow( raw_array[i] - mean, 2 );
          }
          sd = sqrt( (total / (double)size) );
          delete raw_array;
        } else {
          mean = 0;
          sd = 0;
        }
}


xmlNodePtr Prospect_OutGetThreadingNode( xmlDocPtr output_xml_handle, int N) {
        xmlDocPtr doc = (xmlDocPtr) output_xml_handle;
        xmlNodePtr cur_node = xmlDocGetRootElement(doc);
        if ( !strcmp((char *)cur_node->name, "prospectOutput") ) {
                cur_node = cur_node->children;
        }
        int count = 0;  
        while (cur_node != NULL) {
                if (!strcmp((char *)cur_node->name, "threading")) {
                        if ( count == N ) {
                                return cur_node;
                        }
                        count++;
                }
                cur_node = cur_node->next;
        }       
        return NULL;
}


typedef struct {
        int rank;
        float value;    
} sort_feature; 

int cmp_sort_feature(const void *arg1, const void *arg2) {
        sort_feature *stat1 = (sort_feature*)arg1, 
    *stat2 = (sort_feature*)arg2;
        if ( stat1->value < stat2->value ) 
                return 1;
        else 
                return -1;
        return 0;       
}


void ProspectOutput::Sort(char *eqn, bool minimize) {
        ThreadFeatureParse *feature_eqn = ThreadFeatureParse::CompileParse( eqn );
        int N = GetThreadingCount();
        sort_feature *features = (sort_feature *)malloc(sizeof(sort_feature) * N );     
/*
        for (int i = 0; i < N; i++) { 
                ProspectThreadingInfo my_thread_info = GetThreadingInfo( i );
                features[i].value = feature_eqn->RunParse(&my_thread_info) ;
                if (minimize)
                        features[i].value = -(features[i].value);
                features[i].rank = i;
        }
 */
        if ( sort_val == NULL ) {
                sort_val = new float[ N ];
        }
        float *tmp_array = feature_eqn->RunParse( this );
        for (int i = 0; i < N; i++) { 
                features[i].value = tmp_array[i];
                //SetSortVal( i, tmp_array[i]);
                sort_val[i] = tmp_array[i];
                if (minimize)
                        features[i].value = -(features[i].value);
                features[i].rank = i;
        }
        free( tmp_array );
        
        qsort( features, N, sizeof(sort_feature), cmp_sort_feature );   
        if ( reorder )
                free( reorder );
        reorder = (int *)malloc(sizeof(int) * N );
        for (int i = 0; i < N; i++) {
                reorder[i] = features[i].rank;
        }
        free(features); 
        delete feature_eqn;
}


void ProspectOutput::Sort(ThreadFeatureParse *sort_func, bool minimize) {
        int N = GetThreadingCount();
        sort_feature *features = (sort_feature *)malloc(sizeof(sort_feature) * N );     
/*
        ProspectThreadingInfo my_thread_info;
        for (int i = 0; i < N; i++) { 
                GetThreadingInfo( i, &my_thread_info );
                features[i].value = sort_func->RunParse(&my_thread_info) ;
                if (minimize)
                        features[i].value = -(features[i].value);
                features[i].rank = i;
        }
 */
        float *tmp_array = sort_func->RunParse( this );
        for (int i = 0; i < N; i++) { 
                features[i].value = tmp_array[i];
                if (minimize)
                        features[i].value = -(features[i].value);
                features[i].rank = i;
        }
        free( tmp_array );
        
        qsort( features, N, sizeof(sort_feature), cmp_sort_feature );   
        if ( reorder )
                free( reorder );
        reorder = (int *)malloc(sizeof(int) * N );
        for (int i = 0; i < N; i++) {
                reorder[i] = features[i].rank;
        }
        free(features); 
}



void ProspectOutput::Sort(int method) { 
        int N = GetThreadingCount();
        sort_feature *features = (sort_feature *)malloc(sizeof(sort_feature) * N );
        
        xmlDocPtr doc = (xmlDocPtr) output_xml_handle;  
        xmlXPathContextPtr xpathCtx = xmlXPathNewContext(doc); /* Create xpath evaluation context */
        if(xpathCtx == NULL) {
                return;
        }
        
        char searchpath[200] = "/prospectOutput/threading/score/z";
        
        if (method == ScoreStruct::SCORE_RAW)
                strcpy( searchpath, "/prospectOutput/threading/score/raw" );
        if (method == ScoreStruct::SCORE_Z)
                strcpy( searchpath, "/prospectOutput/threading/score/z" );
        if (method == ScoreStruct::SCORE_ZFULL)
                strcpy( searchpath, "/prospectOutput/threading/score/zfull" );
        if (method == ScoreStruct::SCORE_NN)
                strcpy( searchpath, "/prospectOutput/threading/score/nn" );
        if (method == ScoreStruct::SCORE_ZMUTATIONLOG)
                strcpy( searchpath, "/prospectOutput/threading/score/zmutationlog" );
        if (method == ScoreStruct::SCORE_ZSINGLETON)
                strcpy( searchpath, "/prospectOutput/threading/score/zsingleton" );
        if (method == ScoreStruct::SCORE_ZSEC_STRUCT)
                strcpy( searchpath, "/prospectOutput/threading/score/zsec_struct" );
        if (method == ScoreStruct::SCORE_ZTWOBODY)
                strcpy( searchpath, "/prospectOutput/threading/score/ztwobody" );
        if (method == ScoreStruct::SCORE_ZDFIRE)
                strcpy( searchpath, "/prospectOutput/threading/score/zdfire" );
        
        /* Evaluate xpath expression */
        xmlXPathObjectPtr xpathObj = xmlXPathEvalExpression( (const xmlChar*)searchpath, 
                                                                                                                 xpathCtx);             
        if(xpathObj == NULL) {
                xmlXPathFreeContext(xpathCtx); 
                return;
        }
        
        int size = xpathObj->nodesetval->nodeNr;
        
        for(int i = 0; i < size; i++) {
                char * text = (char *)xmlNodeGetContent(xpathObj->nodesetval->nodeTab[i]);
                features[i].value = atof( text );
                if ( method == ScoreStruct::SCORE_RAW)
                        features[i].value = -(features[i].value);
                features[i].rank = i;
                free(text);
        }       
        /* Cleanup of XPath data */
        xmlXPathFreeObject(xpathObj);
        xmlXPathFreeContext(xpathCtx); 
        qsort( features, size, sizeof(sort_feature), cmp_sort_feature );
        if ( reorder )
                free( reorder );
        reorder = (int *)malloc(sizeof(int) * size );
        for (int i = 0; i < size; i++) {
                reorder[i] = features[i].rank;
        }
        free(features);
}

/*
void ProspectOutput::SortScoreUser( bool minimize ) {
        int N = GetThreadingCount();
        sort_feature *features = (sort_feature *)malloc(sizeof(sort_feature) * N );     
        for (int i = 0; i < N; i++) { 
                features[i].value = sort_val[i];
                if (minimize)
                        features[i].value = -(features[i].value);
                features[i].rank = i;
        }
        qsort( features, N, sizeof(sort_feature), cmp_sort_feature );   
        if ( reorder )
                free( reorder );
        reorder = (int *)malloc(sizeof(int) * N );
        for (int i = 0; i < N; i++) {
                reorder[i] = features[i].rank;
        }
        free(features);         
}
*/



void ProspectOutput::SetSortVal( int i, float val  ) {
        if ( sort_val == NULL ) {
                sort_val = new float[ GetThreadingCount() ];
        }
        int num = i;
        sort_val[ num ] = val;
}



float ProspectOutput::GetSortVal( int i  ) {
        if (sort_val == NULL)
                return 0;
        int num = i;
        return sort_val[ num ]; 
}



double get_original_energy( xmlNodePtr node, char *name ) {
        xmlNodePtr child;
        char *text;
        double value;
        
        child = get_first_child( node, name );
        if (child) {
                text = (char *)xmlGetProp(child, (xmlChar *)"e");
                value = atof(text);
                free(text);
                return value;   
        }
        return 0;
}

double get_weights( xmlNodePtr node, char *name ) {
        xmlNodePtr child;
        char *text;
        double value;
        
        child = get_first_child( node, name );
        if (child) {
                text = (char *)xmlNodeGetContent(child);
                value = atof(text);
                free(text);
                return value;
        }
        return 0;
}



int ExtractNodeThreadingInfo( xmlNodePtr threading_node, ProspectThreadingInfo *thread_info) {
        thread_info->template_name =  (char *)xmlGetProp(threading_node, (xmlChar *)"template");
        
        thread_info->xml_node = (void *)threading_node;
        //get score information
        xmlNodePtr score_node = get_first_child( threading_node, "score");
        
        thread_info->scores.mask = 0x00;
        for ( int i = 0; i < ScoreStruct::SCORE_COUNT_TOTAL; i++) {
                xmlNodePtr score_val_node = get_first_child( score_node, ScoreStruct::GetScoreName(i) );
                if (score_val_node != NULL) {
                        char *ptr =  (char *)xmlNodeGetContent(score_val_node);
                        thread_info->scores.values[ i ] = atof(  ptr ) ;
                        free(ptr);
                        
                        ptr =  (char *)xmlGetProp(score_val_node, (xmlChar *)"mean");
                        if ( ptr ) {
                                thread_info->scores.z_mean[ i ] = atof(  ptr ) ;
                                free(ptr);
                        }
                        ptr =  (char *)xmlGetProp(score_val_node, (xmlChar *)"sd");
                        if ( ptr ) {
                                thread_info->scores.z_sd[ i ] = atof(  ptr ) ;
                                free(ptr);
                        }
                        thread_info->scores.mask |= (0x01<<i);
                } 
        }               
        
        //get alignment information
        xmlNodePtr alignment_node = get_first_child( threading_node, "alignment" );
        xmlNodePtr target_ss_node = get_first_child( alignment_node, "query_ss_" );
        xmlNodePtr target_seq_node = get_first_child( alignment_node, "query_seq" );
        xmlNodePtr match_node = get_first_child( alignment_node, "match_str" );
        xmlNodePtr template_ss_node = get_first_child( alignment_node, "templ_ss_" );
        xmlNodePtr template_seq_node = get_first_child( alignment_node, "templ_seq" );
        xmlNodePtr template_energy_node = get_first_child( threading_node, "energy");
        xmlNodePtr template_weight_node = get_first_child(threading_node, "weights");
        xmlNodePtr template_core_node = get_first_child(alignment_node, "core_str_");
        
        thread_info->align_target_seq = (char *)xmlNodeGetContent( target_seq_node );
        thread_info->align_target_ss = (char *)xmlNodeGetContent( target_ss_node );
        thread_info->align_match = (char *)xmlNodeGetContent( match_node );
        thread_info->align_template_seq = (char *)xmlNodeGetContent( template_seq_node );
        thread_info->align_template_ss = (char *)xmlNodeGetContent( template_ss_node );
        
        char *core_align = (char *)xmlNodeGetContent( template_core_node );
        
        for (int j = 0; j < e_count; j++) {
                thread_info->energy[j] = get_original_energy( template_energy_node, e_names[j]);
        }
        for (int j = 0; j < e_count; j++) {
                thread_info->weight[j] = get_weights( template_weight_node, e_names[j]);
        }
        
        int target_pos = 0;
        thread_info->target_len = 0;
        while( thread_info->align_target_seq[ target_pos ] != 0 ) {
                if ( thread_info->align_target_seq[ target_pos ] != '-' ) {
                        thread_info->target_len++;
                }       
                target_pos++;
        }
        
        int template_pos = 0;
        thread_info->template_len = 0;
        while( thread_info->align_template_seq[ template_pos ] != 0 ) {
                if ( thread_info->align_template_seq[ template_pos ] != '-' ) {
                        thread_info->template_len++;
                }       
                template_pos++;
        }
        
        
        thread_info->target_rs = (char *)malloc( thread_info->target_len );
        thread_info->target_ss = (char *)malloc( thread_info->target_len );
        thread_info->template_rs = (char *)malloc( thread_info->template_len );
        thread_info->template_ss = (char *)malloc( thread_info->template_len );
        
        char *core_rs = (char *)malloc( thread_info->template_len );
        
        //thread_info->align.target_align = (int *)malloc(sizeof(int) * thread_info->target_len );
        //thread_info->align.template_align = (int *)malloc(sizeof(int) * thread_info->template_len );  
        //thread_info->align.target_len = thread_info->target_len;
        //thread_info->align.template_len = thread_info->template_len;
        int * target_align = (int *)malloc(sizeof(int) * thread_info->target_len );
        int * template_align = (int *)malloc(sizeof(int) * thread_info->template_len ); 
        
        template_pos = 0;
        target_pos = 0;
        int read_len = strlen(thread_info->align_target_seq);
        for (int i = 0; i < read_len ; i++) {
                if ( target_pos < thread_info->target_len )
                        target_align[ target_pos ] = -1;
                if (  template_pos < thread_info->template_len )
                        template_align[ template_pos ] = -1;
                
                if ( thread_info->align_target_seq[i] != '-' && thread_info->align_template_seq[i] != '-' ) {
                        assert( target_pos < thread_info->target_len );
                        assert(  template_pos < thread_info->template_len );
                        target_align[ target_pos ] = template_pos;
                        template_align[ template_pos ] = target_pos;            
                }
                if ( thread_info->align_target_seq[i] != '-' ) {
                        assert( target_pos < thread_info->target_len );
                        thread_info->target_rs[ target_pos ] = thread_info->align_target_seq[i];
                        thread_info->target_ss[ target_pos ] = thread_info->align_target_ss[i];
                        target_pos++;
                }
                if ( thread_info->align_template_seq[i] != '-' ) {
                        assert(  template_pos < thread_info->template_len );
                        thread_info->template_rs[ template_pos ] = thread_info->align_template_seq[i];          
                        thread_info->template_ss[ template_pos ] = thread_info->align_template_ss[i];           
                        core_rs[ template_pos ] = core_align[i];
                        template_pos++;
                }       
        }
        thread_info->align.SetAlign( thread_info->target_len, target_align, thread_info->template_len, template_align );
        
        //count the cores
        thread_info->core_count = 0;
        for (int i = 0; i < thread_info->template_len; i++) {
                if ( core_rs[i] == '-' ) {
                        if (i == 0 || core_rs[i-1] == ' ') {
                                thread_info->core_count++;
                        }
                }
        }
        
        thread_info->core = (TemplateCoreDefinition *)malloc(sizeof(TemplateCoreDefinition) * thread_info->core_count );
        
        template_pos = 0;
        for (int i = 0; i < thread_info->template_len; i++) {
                if ( core_rs[i] == '-' ) {
                        if (i == 0 || core_rs[i-1] == ' ') {
                                thread_info->core[ template_pos ].head = i;
                                template_pos++;
                        }
                }
                if (core_rs[i] == ' ') {
                        if (i > 0 && core_rs[i-1] == '-') {
                                thread_info->core[ template_pos-1 ].len = i - thread_info->core[ template_pos-1 ].head;
                        }
                }
        }
        if ( core_rs[ thread_info->template_len-1 ] == '-' ) {
                thread_info->core[ template_pos-1 ].len = 
                thread_info->template_len - thread_info->core[ template_pos-1 ].head;
        }
        
        free(core_rs);
        free(core_align);
        
        thread_info->core_res_count = 0;
        for (int i = 0; i < thread_info->core_count; i++) {
                thread_info->core_res_count += thread_info->core[ i ].len;
        }
        
        /*      
                thread_info->nalign_core_res = 0;
         thread_info->nalign_core = 0;
         thread_info->nalign = 0;
         thread_info->nident = 0;
         for (int j = 0; j < thread_info->core_count; j++) {
                 if ( thread_info->align.GetTemplateAlign( thread_info->core[j].head ) != -1 ) {
                         thread_info->nalign_core++;
                 }              
         }
         for (int j = 0; j < thread_info->core_count; j++) {
                 thread_info->nalign_core_res += thread_info->core[j].len;              
         }
         for (int j = 0; j < thread_info->target_len; j++) {
                 if ( thread_info->align.GetTargetAlign(j) != -1 )  {
                         thread_info->nalign ++;
                         if ( thread_info->target_rs[j] == thread_info->template_rs[ thread_info->align.GetTargetAlign(j) ] ) {
                                 thread_info->nident++;
                         }
                 }
         }      
         
         
         {
                 char *nalign_ptr = (char *)xmlGetProp(alignment_node, (xmlChar *)"nalign");
                 char *npair_align_ptr = (char *)xmlGetProp(alignment_node, (xmlChar *)"npair_align");
                 
                 if ( nalign_ptr ) {
                         thread_info->nalign = atoi( nalign_ptr );
                         free( nalign_ptr );
                 }
                 if ( npair_align_ptr ) {
                         thread_info->npair_align = atoi( npair_align_ptr );
                         free( npair_align_ptr );
                 }
                 
         }      
         */
        {
                xmlNodePtr feature_node = get_first_child( threading_node, "features" );                
                if ( feature_node ) {
                        for ( int i = 0; i < thread_info->features.Count(); i++ ) {
                                char *text = get_first_child_content( feature_node, thread_info->features.GetName(i) );
                                if ( text != NULL ) {
                                        thread_info->features.Set(i, atof( text ) );
                                        free( text );
                                }
                        }               
                }
        }
        
        
        xmlNodePtr sidechain_node = get_first_child( threading_node, "sidechain_pack" );
        if ( sidechain_node != NULL ) {
                xmlNodePtr lret_node = get_first_child( sidechain_node, "lret" );
                if ( lret_node ) {
                        char *tmp_ptr = (char *)xmlNodeGetContent( lret_node );
                        thread_info->lret = atof(tmp_ptr);
                        free(tmp_ptr);
                } else {
                        thread_info->lret = 100000;
                }               
                xmlNodePtr atom_count_node = get_first_child( sidechain_node, "atom_count" );
                if ( lret_node ) {
                        char *tmp_ptr = (char *)xmlNodeGetContent( atom_count_node );
                        thread_info->atom_count = atoi(tmp_ptr);
                        free(tmp_ptr);
                } else {
                        thread_info->atom_count = 0;
                }
        } else {
                thread_info->lret = 100000;
                thread_info->atom_count = 0;
        }
        
        xmlNodePtr pdb_node = get_first_child( threading_node, "pdb" );
        if (pdb_node != NULL) {
                thread_info->residue_atoms = (ResidueAtoms *)malloc(sizeof(ResidueAtoms) * thread_info->target_len);
                for (int z = 0; z < thread_info->target_len; z++)
                        thread_info->residue_atoms[z].atom_count = 0;
                char *read_buffer = (char *)xmlNodeGetContent(pdb_node);
                char *cur_line = read_buffer;
                int res_num = -1;
                char cur_residue_name[30];
                strcpy(cur_residue_name, "");   
                char tmp_x_str[30], tmp_y_str[30], tmp_z_str[30];
                strcpy(tmp_x_str, "");
                strcpy(tmp_y_str, "");
                strcpy(tmp_z_str, "");
                while(cur_line[0] != 0) {
                        if ( !strncmp(cur_line, "ATOM", 4) ) {
                                char tmp_str[50];
                                strncpy( tmp_str, &(cur_line[22]), 5 );
                                tmp_str[5] = 0;
                                if (strcmp( cur_residue_name, tmp_str ) ) {
                                        //                                      res_num++;
                                        //residue id//////////////////////////////////////
                                        char residue_id[10];
                                        strncpy( residue_id, &(cur_line[22]), 5 );
                                        residue_id[5] = 0;
                                        //                                      clean_space( residue_id );
                                        res_num = atoi( residue_id );
                                        thread_info->residue_atoms[res_num].atom_count = 0;
                                        strcpy( cur_residue_name, tmp_str );
                                }
                                if (res_num >= 0) {                             
                                        strncpy(tmp_x_str, &(cur_line[30]), 8);
                                        tmp_x_str[8] = 0;
                                        strncpy(tmp_y_str, &(cur_line[38]), 8);
                                        tmp_y_str[8] = 0;
                                        strncpy(tmp_z_str, &(cur_line[46]), 8);
                                        tmp_z_str[8] = 0;
                                        //location information/////////////////////////////
                                        thread_info->residue_atoms[res_num].atom[thread_info->residue_atoms[res_num].atom_count][0] = atof(tmp_x_str);
                                        thread_info->residue_atoms[res_num].atom[thread_info->residue_atoms[res_num].atom_count][1] = atof(tmp_y_str);
                                        thread_info->residue_atoms[res_num].atom[thread_info->residue_atoms[res_num].atom_count][2] = atof(tmp_z_str);
                                        
                                        //Temperature factor/////////////////////////////
                                        strncpy( tmp_str, &(cur_line[60]), 6);
                                        tmp_str[6] = 0;
                                        thread_info->residue_atoms[res_num].temp_factor[thread_info->residue_atoms[res_num].atom_count] = atof( tmp_str );                                      
                                        //strncpy( tmp_str, &(cur_line[17]), 3 );
                                        //tmp_str[3] = 0;
                                        //target_data->residue[res_num].RS = AA3toAA1( tmp_str );                                       
                                        //atom name//////////////////////////////////////
                                        //strncpy( thread_info->residue_atoms[res_num].name[thread_info->residue_atoms[res_num].atom_count], &(cur_line[13]), 3);
                                        //thread_info->residue_atoms[res_num].name[thread_info->residue_atoms[res_num].atom_count][3] = 0;
                                        int atom_namenum = PDBAtomtoNum( &(cur_line[13]) );
                                        thread_info->residue_atoms[res_num].atom_name[ thread_info->residue_atoms[res_num].atom_count ] = 
                                                atom_namenum;                                   
                                        if ( atom_namenum == anN || atom_namenum == anCa || 
                                                 atom_namenum == anC || atom_namenum == anO  ) {
                                                thread_info->residue_atoms[res_num].bb_num[ atom_namenum ] = thread_info->residue_atoms[res_num].atom_count;
                                        }                               
                                        
                                        //FIXME:  need to fill in bb_num stuff
                                        //      clean_space( the_data->atom[i].atom_name );
                                        
                                        /*
                                         the_data->atom[i].chain = str[21];
                                         the_data->atom[i].color[0] = ((float)random()) / (float)LONG_MAX;
                                         the_data->atom[i].color[1] = ((float)random()) / (float)LONG_MAX;
                                         the_data->atom[i].color[2] = ((float)random()) / (float)LONG_MAX;
                                         */
                                        thread_info->residue_atoms[res_num].atom_count++;
                                }
                        }
                        if ( !strncmp(cur_line, "ENDMDL", 6) ) {
                                break;
                        }
                        
                        while (cur_line[0] != '\n' && cur_line[0] != 0) {
                                cur_line++;
                        }
                        if (cur_line[0] == '\n')
                                cur_line++;             
                }
                free(read_buffer);              
        } else {
                thread_info->residue_atoms = NULL;
        }       
        return 0;       
}



ProspectThreadingInfo ProspectOutput::GetThreadingInfo(int N) { 
        ProspectThreadingInfo my_thread;
        GetThreadingInfo(  N, &my_thread );     
        return my_thread;
}


int ProspectOutput::GetThreadingInfo(  const char *name, ProspectThreadingInfo *thread_info) {  
        xmlDocPtr doc = (xmlDocPtr) output_xml_handle;  
        xmlXPathContextPtr xpathCtx = xmlXPathNewContext(doc); /* Create xpath evaluation context */
        if(xpathCtx == NULL) {
                return 1;
        }       
        char searchpath[2000];
        sprintf(searchpath, "/prospectOutput/threading[@template='%s']", name); 
        xmlXPathObjectPtr xpathObj = xmlXPathEvalExpression( (const xmlChar*)searchpath, xpathCtx);             
        if(xpathObj == NULL) {
                xmlXPathFreeContext(xpathCtx); 
                return 1;
        }       
        int size = xpathObj->nodesetval->nodeNr;        
        if ( size == 0 )
                return 1;       
        ExtractNodeThreadingInfo( xpathObj->nodesetval->nodeTab[0], thread_info);               
        thread_info->parent = this;
        xmlXPathFreeObject(xpathObj);
        xmlXPathFreeContext(xpathCtx);          
        return 0;
}

int ProspectOutput::GetThreadingInfo(  int N, ProspectThreadingInfo *thread_info) {
        xmlNodePtr threading_node = Prospect_OutGetThreadingNode( (xmlDocPtr)output_xml_handle, N);
        if ( NULL == threading_node ) {
                return 1;
        }
        ExtractNodeThreadingInfo( threading_node, thread_info);         
        thread_info->parent = this;
        return 0;
}


int ProspectOutput::Get_PDBText(int N, ProspectThreadingInfo *thread_info) {
        xmlNodePtr threading_node = Prospect_OutGetThreadingNode((xmlDocPtr)output_xml_handle, N);
        if ( NULL == threading_node ) {
                return 1;
        }
        xmlNodePtr sidechain_node = get_first_child( threading_node, "sidechain_pack" );
        if ( sidechain_node != NULL ) {
                xmlNodePtr pdb_node = get_first_child( sidechain_node, "pdb" );
                if ( pdb_node != NULL ) {
                        thread_info->pdb_text = (char *)xmlNodeGetContent( pdb_node );
                }               
        }       
        return 0;
}


ProspectThreadingInfo::ProspectThreadingInfo(){
        parent = NULL;
        align_target_seq = NULL;
        align_target_ss = NULL;
        align_match = NULL;
        align_template_seq = NULL;
        align_template_ss = NULL;
        template_name = NULL;
        template_rs = NULL;
        pdb_text = NULL;
        target_rs = NULL;
        template_ss = NULL;
        target_ss = NULL;
        core = NULL;
        xml_node = NULL;
}


ProspectThreadingInfo::~ProspectThreadingInfo() {
        Free(); 
}


void ProspectThreadingInfo::Free(void) {
        if (align_target_seq != NULL)
                free( align_target_seq );
        align_target_seq = NULL;
        if (align_target_ss != NULL)
                free( align_target_ss );
        align_target_ss = NULL;
        if (align_match != NULL)
                free( align_match );
        align_match = NULL;
        if (align_template_seq != NULL)
                free( align_template_seq );
        align_template_seq = NULL;
        if (align_template_ss != NULL)
                free( align_template_ss );      
        align_template_ss = NULL;
        if (template_name != NULL)
                free( template_name );
        template_name = NULL;
        if (template_rs != NULL)
                free(template_rs);
        template_rs = NULL;
        if (target_rs != NULL)
                free(target_rs);
        target_rs = NULL;
        if (template_ss != NULL)
                free(template_ss);
        template_ss = NULL;
        if (target_ss != NULL)
                free(target_ss);
        target_ss = NULL;
        if (core != NULL)
                free(core);
        core = NULL;
        if (pdb_text != NULL)
                free(pdb_text);
        pdb_text = NULL;
        if (residue_atoms != NULL)
                free(residue_atoms);
        residue_atoms = NULL;
        
        parent = NULL;
}



void ProspectThreadingInfo::SetScoreExtra(int extra_num, pValType val ) {

        int score_num = -1;
        switch (extra_num) {
                case 0: {
                        score_num = ScoreStruct::SCORE_USER0;
                        break;
                }
                case 1: {
                        score_num = ScoreStruct::SCORE_USER1;
                        break;
                }
                case 2: {
                        score_num = ScoreStruct::SCORE_USER2;
                        break;
                }
                case 3: {
                        score_num = ScoreStruct::SCORE_USER3;
                        break;
                }
        };
        if ( score_num == -1 )
                return;
        
        scores.SetScore( score_num, val );

        if ( xml_node == NULL )
                return;
        
        //write into the XML, so the change will be saved
        xmlNodePtr threading_node = (xmlNodePtr) xml_node;
        xmlNodePtr score_node = get_first_child( threading_node, "score");
        xmlNodePtr score_val_node = get_first_child( score_node, ScoreStruct::GetScoreName( score_num ) );
        
        if ( score_val_node != NULL ) {
                xmlUnlinkNode(score_val_node);
                xmlFreeNode(score_val_node);    
        }
        char buffer[1000];
        sprintf(buffer, "%g", val);
        score_val_node = xmlNewChild(score_node, NULL, (xmlChar *)ScoreStruct::GetScoreName( score_num ), (xmlChar *)buffer);
}



pValType ProspectThreadingInfo::GetScoreExtra(int extra_num ) {
        if (extra_num == 0)     return scores.GetScore( ScoreStruct::SCORE_USER0 );
        if (extra_num == 1)     return scores.GetScore( ScoreStruct::SCORE_USER1 );
        if (extra_num == 2)     return scores.GetScore( ScoreStruct::SCORE_USER2 );
        if (extra_num == 3)     return scores.GetScore( ScoreStruct::SCORE_USER3 );
        return 0;
}


ProspectThreadingInfo::ProspectThreadingInfo(const ProspectThreadingInfo &rhs) {
        Free();
        Copy(rhs);
}

ProspectThreadingInfo & ProspectThreadingInfo::operator =(const ProspectThreadingInfo &rhs) {
        Free();
        Copy(rhs);
        return *this;
}



void ProspectThreadingInfo::Copy( const ProspectThreadingInfo &rhs) {
        memcpy( this, &rhs, sizeof( ProspectThreadingInfo ) );
        if ( rhs.align_target_seq )
                align_target_seq = strdup( rhs.align_target_seq );
        if ( rhs.align_target_ss )
                align_target_ss = strdup( rhs.align_target_ss );
        if( rhs.align_match )
                align_match = strdup( rhs.align_match );
        if( rhs.align_template_seq )
                align_template_seq = strdup( rhs.align_template_seq );
        if( rhs.align_template_ss )
                align_template_ss = strdup( rhs.align_template_ss );
        if( rhs.template_name )
                template_name = strdup( rhs.template_name );
        if( rhs.template_rs )
                template_rs = strdup( rhs.template_rs );
        if( rhs.target_rs )
                target_rs = strdup( rhs.target_rs );
        if( rhs.template_ss )
                template_ss = strdup( rhs.template_ss );
        if( rhs.target_ss )
                target_ss = strdup( rhs.target_ss );
        if ( rhs.pdb_text )
                pdb_text = strdup( rhs.pdb_text );
        
        if ( core_count > 0 ) {
                core = (TemplateCoreDefinition *)malloc(sizeof(TemplateCoreDefinition) * core_count);
                if ( core ) 
                        memcpy( core, rhs.core, sizeof(TemplateCoreDefinition) * core_count );
        } else {
                core = NULL;
        }
        if ( atom_count > 0 ) {
                residue_atoms = (ResidueAtoms *)malloc(sizeof(ResidueAtoms) * atom_count);
                if ( residue_atoms ) 
                        memcpy( residue_atoms, rhs.residue_atoms, sizeof(ResidueAtoms) * atom_count );
        } else {
                residue_atoms = NULL;
        }
        align = (AlignmentStruct) rhs.align;
        parent = rhs.parent;
        for ( int i=0; i < e_count; i++ ) {
                weight[i] = rhs.weight[i];
                energy[i] = rhs.energy[i];
        }
}


pValType ProspectThreadingInfo::Score_Calc_NN(void ) {
        //format string "%sr %el %e1 %eg %es %ed %sr %fq %ft %fc %am %ac %an %ai"       
        float nn_input[14];
        nn_input[0] = this->scores.values[ ScoreStruct::SCORE_RAW ];
        nn_input[1] = this->energy[ e_mutationlog ];
        nn_input[2] = this->energy[ e_singleton ];
        nn_input[3] = this->energy[ e_gapconst ] + this->energy[ e_gapopen ];
        nn_input[4] = this->energy[ e_sec_struct ];
        nn_input[5] = this->energy[ e_dfire ];
        nn_input[6] = this->scores.values[ ScoreStruct::SCORE_Z ];
        nn_input[7] = this->target_len;
        nn_input[8] = this->template_len;
        nn_input[9] = this->core_count;         
        nn_input[10] = this->features.Get( AlignFeatures::fn_align_core_res );
        nn_input[11] = this->features.Get( AlignFeatures::fn_align_core );
        nn_input[12] = this->features.Get( AlignFeatures::fn_align );
        nn_input[13] = this->features.Get( AlignFeatures::fn_ident );   
        float nn_output;
        foldrec_net(  nn_input, &nn_output, 1 );
        return nn_output;
}


pValType ProspectThreadingInfo::GetAveRaw(void) {
        pValType total = 0;
        pValType nalign = features.Get( AlignFeatures::fn_align );      
        pValType nalign_pair = features.Get( AlignFeatures::fn_align_pair );    
        if ( ((int)nalign) == 0)
                nalign = 1;
        if ( ((int)nalign_pair) == 0)
                nalign_pair = 1;        
        for (int i = 0; i < e_count; i++) {
                if ( e_types[i] == et_single ) {
                        total += (( energy[i] * weight[i] ) / nalign);
                }
                if ( e_types[i] == et_gap ) {
                        total += ( energy[i] * weight[i] ) / target_len;
                }
                if ( e_types[i] == et_pair ) {
                        total += (( energy[i] * weight[i] ) / nalign_pair);                     
                }
        }
        //assert( abs( total ) < 10000 );
        return total;
}


void ProspectThreadingInfo::ExtractTarget(TargetStruct *target_in, int start_in, int stop_in) {
        int start = start_in;
        int stop = stop_in;
        if (target_in == NULL)
                return;
        if (start < 0) 
                start = 0;
        if (stop < 0 || stop > target_len)
                stop = target_len;
        if ( stop <= start )
                return;
        
        int len = stop - start;
        target_in->residue = (TargetResidue *)malloc(sizeof(TargetResidue) * len);
        target_in->len = len;
        target_in->residue_offset = start;
        for (int i = start; i < stop; i++) {
                int cur = i - start;            
                int template_pos = align.target_align[ i ];
                target_in->residue[ cur ].RS = target_rs[i];
                for (int j = 0; j < 21; j++)
                        target_in->residue[ cur ].FREQ[j] = 0;
                target_in->residue[ cur ].FREQ[ AA1toAANum( target_rs[i] ) ] = 1;
                
                if ( template_pos != -1)
                        target_in->residue[ cur ].SS = template_ss[template_pos];
                else 
                        target_in->residue[ cur ].SS = 'L'; //set residues without matching template structure to loop
                for (int j = 0; j < 3; j++) {
                        target_in->residue[ cur ].SS_PROB[j] = 0.0;
                }
                target_in->residue[ cur ].SS_PROB[ SS1toSSNum(target_in->residue[ cur ].SS ) ] = 1.0;
                target_in->residue[ cur ].orig_num = cur;
        }
        if ( residue_atoms ) {
                target_in->structure = new ResidueAtoms[ len ];
                for (int i = start; i < stop; i++) {
                        int cur = i - start;            
                        target_in->structure[ cur ] = residue_atoms[i] ;
                }
        } else {
                target_in->structure = NULL;
        }       
}

---