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





char f_names[f_count][30] = {
        "nalign",
        "nident", 
        "ncontact",
        "nhalfcontact",
        "nalign_pair",
        "nalign_core",
        "nalign_core_res",
        "ralign",
        "rident",
        "rcontact",
        "rhalfcontact", 
        "ralign_pair"
        "rsurface"
};

char *AlignFeatures::GetName(unsigned int i) {
        return f_names[i];
}



ThreadFeatureParse::ThreadFeatureParse() {
        ops = NULL;
        mem = NULL;
}

ThreadFeatureParse::ThreadFeatureParse( feature_parse_op *op_list, int op_count_in, int mem_size ) {
        ops = op_list;
        op_count = op_count_in;
        mem = new double[ mem_size ];
        mem_count = mem_size;
}

ThreadFeatureParse::~ThreadFeatureParse() {
        if (ops)
                delete ops;
        if (mem)
                delete mem;
}


ThreadFeatureParse::ThreadFeatureParse( char *eqn_str ) {
        ThreadFeatureParse *tmp = CompileParse( eqn_str );
        *this = *tmp;
        delete tmp;
}

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

void ThreadFeatureParse::Copy(const ThreadFeatureParse &rhs) {
        op_count = rhs.op_count;
        mem_count = rhs.mem_count;      
        ops = new feature_parse_op[ op_count ];
        mem = new double[ mem_count ];  
        memcpy( ops, rhs.ops, sizeof(feature_parse_op) * op_count );    
}



typedef float (*feature_extract_func)(ProspectThreadingInfo *, int arg);

typedef float * (*feature_array_extract_func)(ProspectOutput *, int arg);

typedef struct feature_struct {
        char tag[6];
        char desc[40];
        feature_extract_func extract_func;
        feature_array_extract_func extract_array_func;
        int extract_arg;
} feature_struct;


#define f_template_len 0
#define f_target_len 1
#define f_core_count 2
#define f_coreres_len 3
#define f_targ_temp_len 4
#define f_temp_targ_len 5




feature_struct feature_array[] = {
        //scores
        //      {"sa", "score average", ThreadFeatureParse::calc_averaw, NULL, 0},
        {"sr", "score raw", ThreadFeatureParse::get_score, ThreadFeatureParse::get_score_array, ScoreStruct::SCORE_RAW},
        {"srz", "score raw", ThreadFeatureParse::get_raw_zscore, ThreadFeatureParse::get_raw_zscore_array, 0},
        {"sz", "score z", ThreadFeatureParse::get_score, ThreadFeatureParse::get_score_array, ScoreStruct::SCORE_Z},
        {"sf", "score zfull", ThreadFeatureParse::get_score, ThreadFeatureParse::get_score_array, ScoreStruct::SCORE_ZFULL},
        //      {"sn", "score neural network", ThreadFeatureParse::get_score, NULL, SCORE_NN},
        {"sel", "score z mutation log", ThreadFeatureParse::get_score, ThreadFeatureParse::get_score_array, ScoreStruct::SCORE_ZMUTATIONLOG},
        {"se1", "score z singleton", ThreadFeatureParse::get_score, ThreadFeatureParse::get_score_array, ScoreStruct::SCORE_ZSINGLETON},
        {"ses", "score z secondary structure", ThreadFeatureParse::get_score, ThreadFeatureParse::get_score_array, ScoreStruct::SCORE_ZSEC_STRUCT},
        {"se2", "score z twobody", ThreadFeatureParse::get_score, ThreadFeatureParse::get_score_array, ScoreStruct::SCORE_ZTWOBODY},
        {"sed", "score z dfire", ThreadFeatureParse::get_score, ThreadFeatureParse::get_score_array, ScoreStruct::SCORE_ZDFIRE},
        /*
        {"su0",  "score user defined", ThreadFeatureParse::get_user_score, ThreadFeatureParse::get_user_score_array, 0},
        {"su1",  "score user defined", ThreadFeatureParse::get_user_score, ThreadFeatureParse::get_user_score_array, 1},
        {"su2",  "score user defined", ThreadFeatureParse::get_user_score, ThreadFeatureParse::get_user_score_array, 2},
        {"su3",  "score user defined", ThreadFeatureParse::get_user_score, ThreadFeatureParse::get_user_score_array, 3},
        */
        {"su0",  "score user defined", ThreadFeatureParse::get_score, ThreadFeatureParse::get_score_array, ScoreStruct::SCORE_USER0},
        {"su1",  "score user defined", ThreadFeatureParse::get_score, ThreadFeatureParse::get_score_array, ScoreStruct::SCORE_USER1},
        {"su2",  "score user defined", ThreadFeatureParse::get_score, ThreadFeatureParse::get_score_array, ScoreStruct::SCORE_USER2},
        {"su3",  "score user defined", ThreadFeatureParse::get_score, ThreadFeatureParse::get_score_array, ScoreStruct::SCORE_USER3},
        
        {"szm", "score z mean", ThreadFeatureParse::get_score_mean, ThreadFeatureParse::get_score_array_mean, ScoreStruct::SCORE_Z},
        {"szs", "score z sd", ThreadFeatureParse::get_score_sd, ThreadFeatureParse::get_score_array_sd, ScoreStruct::SCORE_Z},
        {"sfm", "score zfull mean", ThreadFeatureParse::get_score_mean, ThreadFeatureParse::get_score_array_mean, ScoreStruct::SCORE_ZFULL},
        {"sfs", "score zfull sd", ThreadFeatureParse::get_score_sd, ThreadFeatureParse::get_score_array_sd, ScoreStruct::SCORE_ZFULL},
        {"selm", "score z mutation log mean", ThreadFeatureParse::get_score_mean, ThreadFeatureParse::get_score_array_mean, ScoreStruct::SCORE_ZMUTATIONLOG},
        {"sels", "score z mutation log sd", ThreadFeatureParse::get_score_sd, ThreadFeatureParse::get_score_array_sd, ScoreStruct::SCORE_ZMUTATIONLOG},
        {"se1m", "score z singleton mean", ThreadFeatureParse::get_score_mean, ThreadFeatureParse::get_score_array_mean, ScoreStruct::SCORE_ZSINGLETON},
        {"se1s", "score z singleton sd", ThreadFeatureParse::get_score_sd, ThreadFeatureParse::get_score_array_sd, ScoreStruct::SCORE_ZSINGLETON},
        {"sesm", "score z secondary structure mean", ThreadFeatureParse::get_score_mean, ThreadFeatureParse::get_score_array_mean, ScoreStruct::SCORE_ZSEC_STRUCT},
        {"sess", "score z secondary structure sd", ThreadFeatureParse::get_score_sd, ThreadFeatureParse::get_score_array_sd, ScoreStruct::SCORE_ZSEC_STRUCT},
        {"se2m", "score z twobody mean", ThreadFeatureParse::get_score_mean, ThreadFeatureParse::get_score_array_mean, ScoreStruct::SCORE_ZTWOBODY},
        {"se2s", "score z twobody sd", ThreadFeatureParse::get_score_sd, ThreadFeatureParse::get_score_array_sd, ScoreStruct::SCORE_ZTWOBODY},
        {"sedm", "score z dfire mean", ThreadFeatureParse::get_score_mean, ThreadFeatureParse::get_score_array_mean, ScoreStruct::SCORE_ZDFIRE},
        {"seds", "score z dfire sd", ThreadFeatureParse::get_score_sd, ThreadFeatureParse::get_score_array_sd, ScoreStruct::SCORE_ZDFIRE},
        
        //features
        {"ft", "template length",  ThreadFeatureParse::get_struct_feature, ThreadFeatureParse::get_struct_len_array, f_template_len },
        {"fq", "target length", ThreadFeatureParse::get_struct_feature, ThreadFeatureParse::get_struct_len_array, f_target_len },
        {"fc", "core count", ThreadFeatureParse::get_struct_feature, NULL, f_core_count },
        //      {"fm", "core residue count", ThreadFeatureParse::get_struct_feature, f_coreres_len },
        {"fr", "length ratio (fq/ft)", ThreadFeatureParse::get_struct_feature, NULL, f_targ_temp_len},
        {"fR", "length ratio (ft/fq)", ThreadFeatureParse::get_struct_feature, NULL, f_temp_targ_len},
        
        //alignment info
        {"an", "num align res", ThreadFeatureParse::get_align_feature, ThreadFeatureParse::get_align_feature_array, AlignFeatures::fn_align},
        {"aN", "% align res", ThreadFeatureParse::get_align_feature, ThreadFeatureParse::get_align_feature_array, AlignFeatures::fr_align },
        {"ac", "num align core", ThreadFeatureParse::get_align_feature,ThreadFeatureParse::get_align_feature_array,  AlignFeatures::fn_align_core },
        //{"aC", "% align core", ThreadFeatureParse::get_align_feature, ThreadFeatureParse::get_align_feature_array, AlignFeatures::fr_align_align_core },
        {"ai", "num align ident res", ThreadFeatureParse::get_align_feature,ThreadFeatureParse::get_align_feature_array,  AlignFeatures::fn_ident},
        {"aI", "% align ident res", ThreadFeatureParse::get_align_feature, ThreadFeatureParse::get_align_feature_array, AlignFeatures::fr_ident},
        {"am", "num align core res", ThreadFeatureParse::get_align_feature,ThreadFeatureParse::get_align_feature_array,  AlignFeatures::fn_align_core_res },
        //{"aM", "% align core res", ThreadFeatureParse::get_align_feature, ThreadFeatureParse::get_align_feature_array, AlignFeatures::fr_align_core_res},
        
        {"a2c", "num contact pairs",  ThreadFeatureParse::get_align_feature, ThreadFeatureParse::get_align_feature_array, AlignFeatures::fn_contact },
        {"a2d", "num half contact pairs",  ThreadFeatureParse::get_align_feature, ThreadFeatureParse::get_align_feature_array, AlignFeatures::fn_halfcontact },
        
        {"af", "num target alignment fragments", ThreadFeatureParse::get_target_frag_count, NULL, 0},
        
        //packing info
        //{"pc", "packing atom count", "d"},
        //{"pl", "packing LRET energy", "f"},
        
        //weights
        {"wm", "mutation weight", ThreadFeatureParse::get_weight, ThreadFeatureParse::get_weight_array,  e_mutation },
        {"wl", "mutation log weight", ThreadFeatureParse::get_weight, ThreadFeatureParse::get_weight_array, e_mutationlog },
        {"w1", "singleton weight", ThreadFeatureParse::get_weight, ThreadFeatureParse::get_weight_array,  e_singleton },
        {"ws", "secondary struct weight", ThreadFeatureParse::get_weight, ThreadFeatureParse::get_weight_array, e_sec_struct },
        {"w2", "twobody weight", ThreadFeatureParse::get_weight, ThreadFeatureParse::get_weight_array,  e_twobody },
        {"wd", "dfire weight", ThreadFeatureParse::get_weight, ThreadFeatureParse::get_weight_array,  e_dfire },
        {"wo", "gap open weight", ThreadFeatureParse::get_weight, ThreadFeatureParse::get_weight_array,  e_gapopen },
        {"we", "gap const weight", ThreadFeatureParse::get_weight, ThreadFeatureParse::get_weight_array,  e_gapconst },
        {"wt", "coredel weight", ThreadFeatureParse::get_weight, ThreadFeatureParse::get_weight_array,  e_coredel },
        {"wa", "template singledel weight", ThreadFeatureParse::get_weight, ThreadFeatureParse::get_weight_array,  e_tsingledel },
        {"wi", "template singleins weight",  ThreadFeatureParse::get_weight, ThreadFeatureParse::get_weight_array,  e_tsingleins },
        {"wqa", "target singledel weight", ThreadFeatureParse::get_weight, ThreadFeatureParse::get_weight_array,  e_qsingledel },
        {"wqi", "target singleins weight",  ThreadFeatureParse::get_weight, ThreadFeatureParse::get_weight_array,  e_qsingleins },
        {"wc", "contactdel weight",  ThreadFeatureParse::get_weight, ThreadFeatureParse::get_weight_array,  e_contactdel },
        
        //energy
        {"em", "mutation energy", ThreadFeatureParse::get_energy, ThreadFeatureParse::get_energy_array,  e_mutation },
        {"el", "mutation log energy", ThreadFeatureParse::get_energy, ThreadFeatureParse::get_energy_array,  e_mutationlog },
        {"e1", "singleton energy", ThreadFeatureParse::get_energy, ThreadFeatureParse::get_energy_array,  e_singleton},
        {"es", "secondary struct energy", ThreadFeatureParse::get_energy, ThreadFeatureParse::get_energy_array,  e_sec_struct},
        {"e2", "twobody energy", ThreadFeatureParse::get_energy, ThreadFeatureParse::get_energy_array,  e_twobody },
        {"ed", "dfire energy", ThreadFeatureParse::get_energy, ThreadFeatureParse::get_energy_array,  e_dfire},
        {"et", "coredel energy", ThreadFeatureParse::get_energy, ThreadFeatureParse::get_energy_array,  e_coredel},
        {"eo", "gap open energy", ThreadFeatureParse::get_energy, ThreadFeatureParse::get_energy_array,  e_gapopen},
        {"ee", "gap const energy", ThreadFeatureParse::get_energy, ThreadFeatureParse::get_energy_array,  e_gapconst},
        
        {"ea", "template singledel energy", ThreadFeatureParse::get_energy, ThreadFeatureParse::get_energy_array,  e_tsingledel},
        {"ei", "template singleins energy", ThreadFeatureParse::get_energy, ThreadFeatureParse::get_energy_array,  e_tsingleins},
        {"eqa", "target singledel energy", ThreadFeatureParse::get_energy, ThreadFeatureParse::get_energy_array,  e_qsingledel},
        {"eqi", "target singleins energy", ThreadFeatureParse::get_energy, ThreadFeatureParse::get_energy_array,  e_qsingleins},
        {"ec", "contactdel energy", ThreadFeatureParse::get_energy, ThreadFeatureParse::get_energy_array,  e_contactdel},
        
        //energy * weight
        {"emw", "mutation energy * weight", ThreadFeatureParse::get_energy_weight, ThreadFeatureParse::get_energy_weight_array,  e_mutation },
        {"elw", "mutation log energy * weight", ThreadFeatureParse::get_energy_weight, ThreadFeatureParse::get_energy_weight_array,  e_mutationlog },
        {"e1w", "singleton energy * weight", ThreadFeatureParse::get_energy_weight, ThreadFeatureParse::get_energy_weight_array,  e_singleton},
        {"esw", "secondary struct energy * weight", ThreadFeatureParse::get_energy_weight,  ThreadFeatureParse::get_energy_weight_array, e_sec_struct},
        {"e2w", "twobody energy * weight", ThreadFeatureParse::get_energy_weight, ThreadFeatureParse::get_energy_weight_array,  e_twobody },
        {"edw", "dfire energy * weight", ThreadFeatureParse::get_energy_weight, ThreadFeatureParse::get_energy_weight_array,  e_dfire},
        {"etw", "coredel energy * weight", ThreadFeatureParse::get_energy_weight, ThreadFeatureParse::get_energy_weight_array,  e_coredel},
        {"eow", "gap open energy * weight", ThreadFeatureParse::get_energy_weight, ThreadFeatureParse::get_energy_weight_array,  e_gapopen},
        {"ecw", "gap const energy * weight", ThreadFeatureParse::get_energy_weight, ThreadFeatureParse::get_energy_weight_array,  e_gapconst},
        
        {"eaw", "template singledel energy * weight", ThreadFeatureParse::get_energy_weight, ThreadFeatureParse::get_energy_weight_array,  e_tsingledel},
        {"eiw", "template singleins energy * weight", ThreadFeatureParse::get_energy_weight, ThreadFeatureParse::get_energy_weight_array,  e_tsingleins},
        {"eqaw", "query singledel energy * weight", ThreadFeatureParse::get_energy_weight, ThreadFeatureParse::get_energy_weight_array,  e_qsingledel},
        {"eqiw", "query singleins energy * weight", ThreadFeatureParse::get_energy_weight, ThreadFeatureParse::get_energy_weight_array,  e_qsingleins},
        {"ecw", "contactdel energy * weight", ThreadFeatureParse::get_energy_weight, ThreadFeatureParse::get_energy_weight_array,  e_contactdel}
        
};


int ThreadFeatureParse::GetFeatureID( char *feature_key) {      
        //check numerical features
        for (int i= 0; i < sizeof(feature_array) / sizeof(feature_array[0]); i++) {
                if ( !strcmp( feature_array[i].tag, feature_key ) ) {
                        return i;
                }
        }       
        return -1;
}



float ThreadFeatureParse::GetFeatureFloat(  ProspectThreadingInfo *thread, char *feature_key) { 
        //check numerical features
        for (int i= 0; i < sizeof(feature_array) / sizeof(feature_array[0]); i++) {
                if ( !strcmp( feature_array[i].tag, feature_key ) ) {
                        return feature_array[i].extract_func( thread, feature_array[i].extract_arg );
                }
        }       
        return 0;
}


float ThreadFeatureParse::GetFeatureFloat(  ProspectThreadingInfo *thread, int feature_id) {    
        if ( feature_id < 0 || feature_id > sizeof(feature_array) / sizeof(feature_array[0]) )
                return 0;
        return feature_array[feature_id].extract_func( thread, feature_array[feature_id].extract_arg );
}

char *ThreadFeatureParse::GetFeatureStr( ProspectThreadingInfo *thread, char *feature_key ) {
        
        //check numerical features
        for (int i= 0; i < sizeof(feature_array) / sizeof(feature_array[0]); i++) {
                if ( !strcmp( feature_array[i].tag, feature_key ) ) {
                        char buffer[100];
                        sprintf(buffer, "%g", feature_array[i].extract_func( thread, feature_array[i].extract_arg ) );
                        return strdup(buffer);
                }
        }
        if ( !strcmp( "ss", feature_key ) ) {
                if ( thread->parent == NULL )
                        return NULL;
                char buffer[100];
                int num = thread->parent->GetTemplateNum( thread->template_name ); 
                float value = thread->parent->GetSortVal(num);
                sprintf(buffer, "%g", value );
                return strdup( buffer );        
        }       
        //check str features
        if ( !strcmp( "n", feature_key ) ) {
                return strdup( thread->template_name );
        }       
        return NULL;
}



float ThreadFeatureParse::calc_averaw( ProspectThreadingInfo *thread, int foo) {
        pValType total = 0;
        for (int i = 0; i < e_count; i++) {
                pValType tmp = thread->energy[i] * thread->weight[i];
                if ( e_types[i] == et_single )
                        tmp /= (pValType) thread->features.Get( AlignFeatures::fn_align);
                if ( e_types[i] == et_pair )
                        tmp /= (pValType) thread->features.Get( AlignFeatures::fn_align_pair);
                if (!isnan(tmp))
                        total += tmp;
        }
        return total;   
}       


float ThreadFeatureParse::get_raw_zscore( ProspectThreadingInfo *thread, int feature ) {
        if (thread->parent == NULL)
                return thread->scores.values[ ScoreStruct::SCORE_RAW ];
        float mean, sd;
        thread->parent->GetRawNormalize( mean, sd );    
        return (thread->scores.values[ ScoreStruct::SCORE_RAW ] - mean) / sd;   
}



float ThreadFeatureParse::get_struct_feature( ProspectThreadingInfo *thread, int feature ) {
        switch (feature) {
                case f_template_len:
                        return thread->template_len;
                        break;
                case  f_target_len:
                        return thread->target_len;
                        break;
                case  f_core_count:
                        return thread->core_count;
                        break; 
                        //case  f_coreres_len:
                        //                      break; 
                case  f_targ_temp_len:
                        return ((float)thread->target_len)/((float)thread->template_len);
                        break; 
                case  f_temp_targ_len:
                        return ((float)thread->template_len)/((float)thread->target_len);
                        break;
        }
        return 0;       
}

int get_alpha_count( char *str) {
        int i = 0;
        int count = 0;
        while ( str[i] ) {
                if ( isalpha( str[i] ) )
                        count++;
                i++;
        }
        return count;
}

/*
float *ThreadFeatureParse::get_user_score_array( ProspectOutput *threads, int feature ) {
        float *tmp = (float *)malloc(sizeof(float) * threads->GetThreadingCount() ) ;
        for (int i = 0; i < threads->GetThreadingCount(); i++ ) {
                tmp[i] = threads->GetScoreExtra(i, feature);            
        }
        return tmp;
}
*/


float *ThreadFeatureParse::get_struct_len_array( ProspectOutput *threads, int feature ) {
        char searchpath[200] = "/prospectOutput/threading/alignment/templ_seq";
        if (  f_target_len == feature )
                sprintf(searchpath, "/prospectOutput/threading/alignment/query_seq" );
        
        int N = threads->GetThreadingCount();
        float *scores = (float *)malloc(sizeof(float) * N );    
        
        xmlDocPtr doc = (xmlDocPtr) threads->output_xml_handle; 
        xmlXPathContextPtr xpathCtx = xmlXPathNewContext(doc); /* Create xpath evaluation context */
        if(xpathCtx == NULL) {
                return NULL;
        }                       
        /* Evaluate xpath expression */
        xmlXPathObjectPtr xpathObj = xmlXPathEvalExpression( (const xmlChar*)searchpath, xpathCtx);             
        if(xpathObj == NULL) {
                xmlXPathFreeContext(xpathCtx); 
                return NULL;
        }       
        int size = xpathObj->nodesetval->nodeNr;        
        for(int i = 0; i < size; i++) {
                char * text;
                //if ( threads->reorder )
                //      text = (char *)xmlNodeGetContent(xpathObj->nodesetval->nodeTab[ threads->reorder[i] ]);
                //else 
                text = (char *)xmlNodeGetContent(xpathObj->nodesetval->nodeTab[ i ]);
                scores[i] = get_alpha_count( text );
                free(text);
        }
        /* Cleanup of XPath data */
        xmlXPathFreeObject(xpathObj);
        xmlXPathFreeContext(xpathCtx); 
        return scores;
}


float *ThreadFeatureParse::get_raw_zscore_array( ProspectOutput *threads, int feature ) {
        float mean, sd;
        threads->GetRawNormalize( mean, sd );
        float *array = get_score_array( threads, ScoreStruct::SCORE_RAW );
        int N = threads->GetThreadingCount();
        for ( int i = 0; i< N; i++) {
                array[i] = (array[i] - mean) / sd;
        }       
        return array;
}


float* ThreadFeatureParse::get_score_array( ProspectOutput *threads, int feature ) {
        char searchpath[200] = "/prospectOutput/threading/score/z";
        
        sprintf(searchpath, "/prospectOutput/threading/score/%s", ScoreStruct::GetScoreName( feature ) );       
        /*
        switch ( feature ) {
                case ScoreStruct::SCORE_RAW:
                        strcpy( searchpath, "/prospectOutput/threading/score/raw" );
                        break;
                case ScoreStruct::SCORE_Z:
                        strcpy( searchpath, "/prospectOutput/threading/score/z" );
                        break;
                case ScoreStruct::SCORE_ZFULL:
                        strcpy( searchpath, "/prospectOutput/threading/score/zfull" );
                        break;
                case ScoreStruct::SCORE_NN:
                        strcpy( searchpath, "/prospectOutput/threading/score/nn" );
                        break;
                case ScoreStruct::SCORE_ZMUTATIONLOG:
                        strcpy( searchpath, "/prospectOutput/threading/score/zmutationlog" );
                        break;
                case ScoreStruct::SCORE_ZSINGLETON:
                        strcpy( searchpath, "/prospectOutput/threading/score/zsingleton" );
                        break;
                case ScoreStruct::SCORE_ZSEC_STRUCT:
                        strcpy( searchpath, "/prospectOutput/threading/score/zsec_struct" );
                        break;
                case ScoreStruct::SCORE_ZTWOBODY:
                        strcpy( searchpath, "/prospectOutput/threading/score/ztwobody" );
                        break;
                case ScoreStruct::SCORE_ZDFIRE:
                        strcpy( searchpath, "/prospectOutput/threading/score/zdfire" );
                        break;
                default:
                        return NULL;
        }       
         */
        int N = threads->GetThreadingCount();
        float *scores = (float *)malloc(sizeof(float) * N );    
        xmlDocPtr doc = (xmlDocPtr) threads->output_xml_handle; 
        xmlXPathContextPtr xpathCtx = xmlXPathNewContext(doc); /* Create xpath evaluation context */
        if(xpathCtx == NULL) {
                return NULL;
        }                       
        /* Evaluate xpath expression */
        xmlXPathObjectPtr xpathObj = xmlXPathEvalExpression( (const xmlChar*)searchpath, xpathCtx);             
        if(xpathObj == NULL) {
                xmlXPathFreeContext(xpathCtx); 
                return NULL;
        }       
        int size = xpathObj->nodesetval->nodeNr;        
        for(int i = 0; i < size; i++) {
                char * text;
                //if ( threads->reorder )
                //      text = (char *)xmlNodeGetContent(xpathObj->nodesetval->nodeTab[ threads->reorder[i] ]);
                //else 
                text = (char *)xmlNodeGetContent(xpathObj->nodesetval->nodeTab[ i ]);
                scores[i] = atof( text );
                free(text);
        }       
        /* Cleanup of XPath data */
        xmlXPathFreeObject(xpathObj);
        xmlXPathFreeContext(xpathCtx); 
        return scores;
}


float* ThreadFeatureParse::get_score_array_mean( ProspectOutput *threads, int feature ) {
        char searchpath[200] = "/prospectOutput/threading/score/z@mean";
        switch ( feature ) {
                case ScoreStruct::SCORE_RAW:
                        strcpy( searchpath, "/prospectOutput/threading/score/raw@mean" );
                        break;
                case ScoreStruct::SCORE_Z:
                        strcpy( searchpath, "/prospectOutput/threading/score/z@mean" );
                        break;
                case ScoreStruct::SCORE_ZFULL:
                        strcpy( searchpath, "/prospectOutput/threading/score/zfull@mean" );
                        break;
                case ScoreStruct::SCORE_NN:
                        strcpy( searchpath, "/prospectOutput/threading/score/nn@mean" );
                        break;
                case ScoreStruct::SCORE_ZMUTATIONLOG:
                        strcpy( searchpath, "/prospectOutput/threading/score/zmutationlog@mean" );
                        break;
                case ScoreStruct::SCORE_ZSINGLETON:
                        strcpy( searchpath, "/prospectOutput/threading/score/zsingleton@mean" );
                        break;
                case ScoreStruct::SCORE_ZSEC_STRUCT:
                        strcpy( searchpath, "/prospectOutput/threading/score/zsec_struct@mean" );
                        break;
                case ScoreStruct::SCORE_ZTWOBODY:
                        strcpy( searchpath, "/prospectOutput/threading/score/ztwobody@mean" );
                        break;
                case ScoreStruct::SCORE_ZDFIRE:
                        strcpy( searchpath, "/prospectOutput/threading/score/zdfire@mean" );
                        break;
                default:
                        return NULL;
        }       
        int N = threads->GetThreadingCount();
        float *scores = (float *)malloc(sizeof(float) * N );    
        xmlDocPtr doc = (xmlDocPtr) threads->output_xml_handle; 
        xmlXPathContextPtr xpathCtx = xmlXPathNewContext(doc); /* Create xpath evaluation context */
        if(xpathCtx == NULL) {
                return NULL;
        }                       
        /* Evaluate xpath expression */
        xmlXPathObjectPtr xpathObj = xmlXPathEvalExpression( (const xmlChar*)searchpath, xpathCtx);             
        if(xpathObj == NULL) {
                xmlXPathFreeContext(xpathCtx); 
                return NULL;
        }       
        int size = xpathObj->nodesetval->nodeNr;        
        for(int i = 0; i < size; i++) {
                char * text;
                //if ( threads->reorder )
                //      text = (char *)xmlNodeGetContent(xpathObj->nodesetval->nodeTab[ threads->reorder[i] ]);
                //else 
                text = (char *)xmlNodeGetContent(xpathObj->nodesetval->nodeTab[ i ]);
                scores[i] = atof( text );
                free(text);
        }       
        /* Cleanup of XPath data */
        xmlXPathFreeObject(xpathObj);
        xmlXPathFreeContext(xpathCtx); 
        return scores;
}



float* ThreadFeatureParse::get_score_array_sd( ProspectOutput *threads, int feature ) {
        char searchpath[200] = "/prospectOutput/threading/score/z@sd";
        switch ( feature ) {
                case ScoreStruct::SCORE_RAW:
                        strcpy( searchpath, "/prospectOutput/threading/score/raw@sd" );
                        break;
                case ScoreStruct::SCORE_Z:
                        strcpy( searchpath, "/prospectOutput/threading/score/z@sd" );
                        break;
                case ScoreStruct::SCORE_ZFULL:
                        strcpy( searchpath, "/prospectOutput/threading/score/zfull@sd" );
                        break;
                case ScoreStruct::SCORE_NN:
                        strcpy( searchpath, "/prospectOutput/threading/score/nn@sd" );
                        break;
                case ScoreStruct::SCORE_ZMUTATIONLOG:
                        strcpy( searchpath, "/prospectOutput/threading/score/zmutationlog@sd" );
                        break;
                case ScoreStruct::SCORE_ZSINGLETON:
                        strcpy( searchpath, "/prospectOutput/threading/score/zsingleton@sd" );
                        break;
                case ScoreStruct::SCORE_ZSEC_STRUCT:
                        strcpy( searchpath, "/prospectOutput/threading/score/zsec_struct@sd" );
                        break;
                case ScoreStruct::SCORE_ZTWOBODY:
                        strcpy( searchpath, "/prospectOutput/threading/score/ztwobody@sd" );
                        break;
                case ScoreStruct::SCORE_ZDFIRE:
                        strcpy( searchpath, "/prospectOutput/threading/score/zdfire@sd" );
                        break;
                default:
                        return NULL;
        }       
        int N = threads->GetThreadingCount();
        float *scores = (float *)malloc(sizeof(float) * N );    
        xmlDocPtr doc = (xmlDocPtr) threads->output_xml_handle; 
        xmlXPathContextPtr xpathCtx = xmlXPathNewContext(doc); /* Create xpath evaluation context */
        if(xpathCtx == NULL) {
                return NULL;
        }                       
        /* Evaluate xpath expression */
        xmlXPathObjectPtr xpathObj = xmlXPathEvalExpression( (const xmlChar*)searchpath, xpathCtx);             
        if(xpathObj == NULL) {
                xmlXPathFreeContext(xpathCtx); 
                return NULL;
        }       
        int size = xpathObj->nodesetval->nodeNr;        
        for(int i = 0; i < size; i++) {
                char * text;
                //if ( threads->reorder )
                //      text = (char *)xmlNodeGetContent(xpathObj->nodesetval->nodeTab[ threads->reorder[i] ]);
                //else 
                text = (char *)xmlNodeGetContent(xpathObj->nodesetval->nodeTab[ i ]);
                scores[i] = atof( text );
                free(text);
        }       
        /* Cleanup of XPath data */
        xmlXPathFreeObject(xpathObj);
        xmlXPathFreeContext(xpathCtx); 
        return scores;
}


float* ThreadFeatureParse::get_weight_array( ProspectOutput *threads, int feature ) {
        char searchpath[200];
        sprintf( searchpath,  "/prospectOutput/threading/weights/%s", e_names[ feature ] );
        int N = threads->GetThreadingCount();
        float *scores = (float *)malloc(sizeof(float) * N );    
        xmlDocPtr doc = (xmlDocPtr) threads->output_xml_handle; 
        xmlXPathContextPtr xpathCtx = xmlXPathNewContext(doc); /* Create xpath evaluation context */
        if(xpathCtx == NULL) {
                return NULL;
        }                       
        /* Evaluate xpath expression */
        xmlXPathObjectPtr xpathObj = xmlXPathEvalExpression( (const xmlChar*)searchpath, xpathCtx);             
        if(xpathObj == NULL) {
                xmlXPathFreeContext(xpathCtx); 
                return NULL;
        }       
        int size = xpathObj->nodesetval->nodeNr;        
        for(int i = 0; i < size; i++) {
                char * text;
                //if ( threads->reorder )
                //      text = (char *)xmlNodeGetContent(xpathObj->nodesetval->nodeTab[ threads->reorder[i] ]);
                //else 
                text = (char *)xmlNodeGetContent(xpathObj->nodesetval->nodeTab[ i ]);
                scores[i] = atof( text );
                free(text);
        }       
        /* Cleanup of XPath data */
        xmlXPathFreeObject(xpathObj);
        xmlXPathFreeContext(xpathCtx); 
        return scores;
}

float* ThreadFeatureParse::get_energy_array( ProspectOutput *threads, int feature ) {
        char searchpath[200];
        sprintf( searchpath,  "/prospectOutput/threading/energy/%s", e_names[ feature ] );
        int N = threads->GetThreadingCount();
        float *scores = (float *)malloc(sizeof(float) * N );    
        xmlDocPtr doc = (xmlDocPtr) threads->output_xml_handle; 
        xmlXPathContextPtr xpathCtx = xmlXPathNewContext(doc); /* Create xpath evaluation context */
        if(xpathCtx == NULL) {
                return NULL;
        }                       
        /* Evaluate xpath expression */
        xmlXPathObjectPtr xpathObj = xmlXPathEvalExpression( (const xmlChar*)searchpath, xpathCtx);             
        if(xpathObj == NULL) {
                xmlXPathFreeContext(xpathCtx); 
                return NULL;
        }       
        int size = xpathObj->nodesetval->nodeNr;        
        for(int i = 0; i < size; i++) {
                char * text;
                //if ( threads->reorder )
                //      text = (char *)xmlNodeGetContent(xpathObj->nodesetval->nodeTab[ threads->reorder[i] ]);
                //else 
                text = (char *)xmlNodeGetContent(xpathObj->nodesetval->nodeTab[ i ]);
                scores[i] = atof( text );
                free(text);
        }       
        /* Cleanup of XPath data */
        xmlXPathFreeObject(xpathObj);
        xmlXPathFreeContext(xpathCtx); 
        return scores;
}

float* ThreadFeatureParse::get_energy_weight_array( ProspectOutput *threads, int feature ) {
        float *tmp1 = get_energy_array( threads, feature );
        float *tmp2 = get_weight_array( threads, feature );
        int N = threads->GetThreadingCount();
        for (int i = 0; i < N; i++) {
                tmp1[i] *= tmp2[i];
        }
        free(tmp2);
        return tmp1;
}



float* ThreadFeatureParse::get_align_feature_array( ProspectOutput *threads, int feature ) {
        char searchpath[200];
        sprintf( searchpath,  "/prospectOutput/threading/features/%s", AlignFeatures::GetName( feature ) );
        int N = threads->GetThreadingCount();
        float *scores = (float *)malloc(sizeof(float) * N );    
        xmlDocPtr doc = (xmlDocPtr) threads->output_xml_handle; 
        xmlXPathContextPtr xpathCtx = xmlXPathNewContext(doc); /* Create xpath evaluation context */
        if(xpathCtx == NULL) {
                return NULL;
        }                       
        /* Evaluate xpath expression */
        xmlXPathObjectPtr xpathObj = xmlXPathEvalExpression( (const xmlChar*)searchpath, xpathCtx);             
        if(xpathObj == NULL) {
                xmlXPathFreeContext(xpathCtx); 
                return NULL;
        }       
        int size = xpathObj->nodesetval->nodeNr;        
        for(int i = 0; i < size; i++) {
                char * text;
                //if ( threads->reorder )
                //      text = (char *)xmlNodeGetContent(xpathObj->nodesetval->nodeTab[ threads->reorder[i] ]);
                //else 
                text = (char *)xmlNodeGetContent(xpathObj->nodesetval->nodeTab[ i ]);
                scores[i] = atof( text );
                free(text);
        }       
        /* Cleanup of XPath data */
        xmlXPathFreeObject(xpathObj);
        xmlXPathFreeContext(xpathCtx); 
        return scores;
}


//declaration to a private function found in seq_io.cc
void append_to_str(char **str1, char *str2);

char *ThreadFeatureParse::ParseFeatureStr( ProspectThreadingInfo *thread_info, char *format_str ) {
        
        char *out_str = NULL;
        
        for (int j = 0; j < strlen(format_str); j++) {
                if (format_str[j] == '%') {
                        j++;
                        int k = j;
                        while ( format_str[k] != '\n' && 
                                        format_str[k] != ' '  && 
                                        format_str[k] != '\t' && 
                                        format_str[k] != '\\' && 
                                        format_str[k] != '%'  && 
                                        format_str[k] != 0)
                                k++;
                        char format_word[201];
                        for (int l = 0; l < k-j && l < 200; l++) {
                                format_word[l] = format_str[l+j];
                                format_word[l+1] = 0;                                                   
                        }
                        //print_feature( &threads[i], format_word );
                        char *tmp = GetFeatureStr( thread_info, format_word);
                        if ( tmp ) {
                                append_to_str( &out_str, tmp );
                                free( tmp );
                        }
                        j = k-1;
                } else if (format_str[j] == '\\') {
                        j++;
                        switch (format_str[j]) {
                                case 't':
                                        //printf("\t");
                                        append_to_str( &out_str, "\t" );
                                        break;
                        }
                } else {
                        char buffer[10];
                        sprintf(buffer, "%c", format_str[j]);
                        append_to_str( &out_str, buffer );
                }
        }
        return out_str;
}


float ThreadFeatureParse::ParseFeatureEqn(ProspectThreadingInfo* thread, char* eqn_str) {
        ThreadFeatureParse*  tmp_parse = ThreadFeatureParse::CompileParse( eqn_str );
        float tmp_float = tmp_parse->RunParse( thread );
        delete tmp_parse;
        return tmp_float; 
}


//BUG:  ParseFeatureEqn links into BISON code with alot of global variables, and therefore is not re-entrant
// be aware when doing parallel processes

#include <list>

extern int mem_stack;
extern std::list<feature_parse_op> op_stack;

extern int yy_scan_string(const char *);
int yyparse (void);



ThreadFeatureParse* ThreadFeatureParse::CompileParse( char *eqn_str ) {
        char *buffer = strdup( eqn_str );
        mem_stack = 0;
        op_stack.empty();  
        yy_scan_string (buffer);
        yyparse();  
        int c = 0;
        feature_parse_op *ops_tmp = new feature_parse_op[ op_stack.size() ];
        for ( std::list<feature_parse_op>::iterator i = op_stack.begin(); i != op_stack.end(); ++i) {
                ops_tmp[c] = (*i);
                c++;
        }
        free( buffer );
        ThreadFeatureParse *out_parse = new ThreadFeatureParse( ops_tmp, c, mem_stack);  
        mem_stack = 0;
        op_stack.empty(); 
        return out_parse;
}



float ThreadFeatureParse::RunFeatureProgram( ProspectThreadingInfo *thread, float *feature_vals ) {
        if ( thread == NULL && feature_vals == NULL  )
                return 0;       
        for (int i = 0; i < op_count; i++) {
                switch (ops[i].op_code) {
                        case FPOP_SET:
#ifdef PRINT_FPOPS
                                printf("SET %d\t%f\n", ops[i].mem_store, ops[i].op_arg.val );
#endif
                                mem[ ops[i].mem_store ] = ops[i].op_arg.val;
                                break;
                        case FPOP_GET:
                                if ( thread ) {
                                        mem[ ops[i].mem_store ] = ThreadFeatureParse::GetFeatureFloat( thread, ops[i].op_arg.reg[0] );
                                } else {
                                        mem[ ops[i].mem_store ] = feature_vals[ ops[i].op_arg.reg[0] ];
                                }
#ifdef PRINT_FPOPS
                                printf("GET %d\t%d(%f)\n", ops[i].mem_store, ops[i].op_arg.reg[0], mem[ ops[i].mem_store ]  );
#endif
                                break;
                        case FPOP_ADD:
                                mem[ ops[i].mem_store ] = mem[ ops[i].op_arg.reg[0] ] + mem[ ops[i].op_arg.reg[1] ];
#ifdef PRINT_FPOPS
                                printf("ADD %d\t%d\t%d\t(%f)\n", ops[i].mem_store, ops[i].op_arg.reg[0], ops[i].op_arg.reg[1], mem[ ops[i].mem_store ]  );                              
#endif
                                break;
                        case FPOP_SUB:
                                mem[ ops[i].mem_store ] = mem[ ops[i].op_arg.reg[0] ] - mem[ ops[i].op_arg.reg[1] ];
#ifdef PRINT_FPOPS
                                printf("SUB %d\t%d\t%d\n", ops[i].mem_store, ops[i].op_arg.reg[0], ops[i].op_arg.reg[1] );                              
#endif
                                break;
                        case FPOP_MUL:
                                mem[ ops[i].mem_store ] = mem[ ops[i].op_arg.reg[0] ] * mem[ ops[i].op_arg.reg[1] ];
#ifdef PRINT_FPOPS
                                printf("MUL %d\t%d\t%d\n", ops[i].mem_store, ops[i].op_arg.reg[0], ops[i].op_arg.reg[1] );                              
#endif
                                break;
                        case FPOP_DIV:
                                mem[ ops[i].mem_store ] = mem[ ops[i].op_arg.reg[0] ] / mem[ ops[i].op_arg.reg[1] ];
#ifdef PRINT_FPOPS
                                printf("DIV %d\t%d\t%d\n", ops[i].mem_store, ops[i].op_arg.reg[0], ops[i].op_arg.reg[1] );                              
#endif
                                break;
                        case FPOP_POW:
                                mem[ ops[i].mem_store ] = pow( mem[ ops[i].op_arg.reg[0] ],  mem[ ops[i].op_arg.reg[1] ] );
#ifdef PRINT_FPOPS
                                printf("POW %d\t%d\t%d\n", ops[i].mem_store, ops[i].op_arg.reg[0], ops[i].op_arg.reg[1] );      
#endif
                                break;
                        case FPOP_SRT:
                                mem[ ops[i].mem_store ] = sqrt( mem[ ops[i].op_arg.reg[0] ] );
                                break;
                        case FPOP_LOG:
                                mem[ ops[i].mem_store ] = log( mem[ ops[i].op_arg.reg[0] ] );
                                break;
                        case FPOP_SIN:
                                mem[ ops[i].mem_store ] = sin( mem[ ops[i].op_arg.reg[0] ] );
                                break;
                        case FPOP_COS:
                                mem[ ops[i].mem_store ] = cos( mem[ ops[i].op_arg.reg[0] ] );
                                break;
                        case FPOP_NEG:
                                mem[ ops[i].mem_store ] = -(mem[ ops[i].op_arg.reg[0] ]);
#ifdef PRINT_FPOPS
                                printf("NEG %d\t%d\n", ops[i].mem_store, ops[i].op_arg.reg[0] );        
#endif
                                break;
                        case FPOP_ABS:
                                mem[ ops[i].mem_store ] = fabs( mem[ ops[i].op_arg.reg[0] ] );
                                break;
                        case FPOP_RES:                                  
#ifdef PRINT_FPOPS
                                printf("RES\t%d\t%f\n", ops[i].op_arg.reg[0], mem[ ops[i].op_arg.reg[0] ]  );
#endif
                                return mem[ ops[i].op_arg.reg[0] ]  ;
                                break;
                }
        }
        return 0;
        
        
}


float ThreadFeatureParse::RunParse( ProspectThreadingInfo *thread ) {
        return RunFeatureProgram( thread );
}


float *ThreadFeatureParse::RunParse( ProspectOutput *threads ) {
        int N = threads->GetThreadingCount();
        int feature_id_count = (sizeof(feature_array) / sizeof(feature_array[0]) );     
        float **feature_val = (float **)malloc( sizeof(float *) * feature_id_count );
        memset( feature_val, 0, sizeof(float *) * feature_id_count);
    //figure out which of the features need to be extracted
        for (int i = 0; i < op_count; i++) {
                if ( ops[i].op_code == FPOP_GET ) {
                        int feature_id = ops[i].op_arg.reg[0];
                        if ( feature_array[feature_id].extract_array_func != NULL ) {
                                feature_val[ feature_id ] = feature_array[feature_id].extract_array_func( threads, feature_array[feature_id].extract_arg );
                        }
                }
        }
        float *array = new float[ N ];  
        float tmp_array[ feature_id_count ];
        for (int i = 0; i < N; i++) {
                for (int j = 0; j < feature_id_count; j++) {
                        if ( feature_val[j] != NULL )
                                tmp_array[j] = feature_val[j][i];
                        else 
                                tmp_array[j] = 0;
                }
                array[i] = RunFeatureProgram( NULL, tmp_array );
        }       
        for (int j = 0; j < feature_id_count; j++) {
                if ( feature_val[j] != NULL )
                        free( feature_val[j] );
        }
        free(feature_val);
        return array;
}

---