src/common/prospect_param.cc

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include "open_prospect.h"

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


#ifdef INT_THREADING
//#include "src/int_threading/integer_thread.h"
#include "src/intpro_threading/intpro_threading.h"
#endif

#ifdef DIVCON_THREADING
#include "src/divcon_threading/divcon_threading.h"
#endif


//#include "src/score_functions/score_functions.h"
#include "src/dynamic_threading/dynamic_threading.h"
#ifdef TREEDECOMP_THREADING
#include "src/treedecomp_threading/treedecomp_threading.h"
#endif



//Threading function access




typedef struct {
        threading_function   threading;
        ThreaderDesc     *desc;
        char   name[30];
} threading_procedure_structure;


threading_procedure_structure threaders[] = {
#ifdef INT_THREADING
#ifdef USE_COIN
        {intpro_thread, &intpro_thread_desc, "integer"},
        {intpro_thread_div, &intpro_thread_div_desc, "integer_div"},
        {intpro_thread_sift, &intpro_thread_sift_desc, "integer_sift"},
        {intpro_thread_coredel, &intpro_thread_coredel_desc, "integer_coredel"},
#endif
#ifdef USE_GTMIP
        {intpro_thread_gtmip, &intpro_thread_gtmip_desc, "integer_gtmip"},
#endif
#endif
#ifdef DIVCON_THREADING
        {divcon_thread, &divcon_thread_desc, "divcon"},
#endif
#ifdef TREEDECOMP_THREADING
        {treedecomp_thread, &treedecomp_thread_desc, "treedecomp"},
#endif
        {dynamic_thread, &dynamic_thread_desc, "dynamic"}
};


threading_function ProspectParam::ThreaderFind(char *request_name) {
        for (unsigned int i = 0; i < sizeof(threaders)/sizeof(threaders[0]); i++) {
                if ( !strcmp(threaders[i].name, request_name) ) {
                        return threaders[i].threading;
                }
        }
        return NULL;
}




int ProspectParam::Weight_Find(char *request_name, WeightArray weight) {
        
        char conf_path[2000], buffer[3000];
        sprintf(conf_path, "%s/%s", ProspectPath, WEIGHT_CONF_FILE);
        
        FILE *file = fopen( conf_path, "r" );
        if (file == NULL)
                return 1;
                
        while ( fgets( buffer, sizeof(buffer), file) ) {        
                char *ptr = strtok( buffer, " :\t\n\r" );
                if ( !strcmp( request_name, ptr ) ) {
                        for (int i = 0; i < e_count; i++) 
                                weight[i] = 0.0;
                        while ( NULL != (ptr = strtok( NULL, " :\t\n\r" ) ) ) {
                                if ( !strncmp( ptr, "-w", 2 ) ) {
                                        int cur_e = -1;
                                        for (int i = 0; i < e_count; i++) {
                                                if (!strcmp( &(ptr[2]), e_names[i]) ) {
                                                        cur_e = i;      
                                                }       
                                        }
                                        if (cur_e != -1 && NULL != (ptr = strtok( NULL, " :\t\n\r" ) ) ) {
                                                weight[ cur_e ] = atof( ptr );                                          
                                        } 
                                }
                        }
                        return 0;
                }
        }       
        return 1;
}


PopsParam::PopsParam(ProspectParam *in_parent) : ParamElement( in_parent ) {}


int PopsParam::LoadParam(char *path) {
        FILE *file = fopen(path, "r");
        if ( file == NULL ) {
                return 1;
        }
        char buffer[2000];
        fgets( buffer, sizeof(buffer), file);
        sscanf( buffer, "%d %f %f %f %f", &param_count, &b_ij, &p_j, &p_14, &p_gt14);
        //printf( "%d %f %f %f %f\n", param_count, b_ij, p_j, p_14, p_gt14);
        if (param_count < 0 || param_count > 1000) {
                param_count = 0;
                return 1;
        }
        params = (PopsParamsElement *)malloc(sizeof(PopsParamsElement) * param_count);
        for ( int i = 0; i < param_count; i++) {
                fgets( buffer, sizeof(buffer), file);
                char atom[2000], atom_name[2000], residue[2000];
                float radius, p_i;
                sscanf( buffer, "%s %s %s %f %f ", atom, atom_name, residue, &radius, &p_i );           
                //printf( "%s %s %s %f %f \n", atom, atom_name, residue, radius, p_i );
                params[i].radius = radius;
                params[i].p_i = p_i;
                params[i].atom_num = PDBAtomtoNum( atom_name );
                params[i].rs_num = AA1toAANum( AA3toAA1( residue ) );
        }
        fclose(file);
        return 0;
}


int ProspectParam::InitPops(void) {
        char buffer[3000];
        
        sprintf(buffer, "%s/%s", ProspectPath , "POPS-R-DNA.dat" );
        pops_r_param = new PopsParam( this );
        pops_r_param->LoadParam( buffer );
        pops_param = new PopsParam( this );
        sprintf(buffer, "%s/%s", ProspectPath , "POPS-A-DNA.dat" );
        pops_param->LoadParam( buffer );
        
        return 0;
}


PopsParam* ProspectParam::GetPopsParam(void) {
        if ( pops_param == NULL )
                InitPops();
        return pops_param;
}

PopsParam* ProspectParam::GetPopsRParam(void) {
        if ( pops_param == NULL )
                InitPops();
        return pops_r_param;
}


ProspectParam::~ProspectParam() {
        
        
}



ProspectParam::ProspectParam(char *base_path) {
        
        ProspectPath = NULL;
        residue_info = NULL;
        rotamer_array = NULL;
        atom_data = NULL;
        atom_data_count = 0;
        pops_param = NULL;
        pops_r_param = NULL;

        
        char search_paths[4][4000] = {
                PROSPECT_DATADIR,
                "/usr/local/share/prospect",
                "/usr/share/prospect/",
                "/opt/share/prospect/"
        };
        
        if (base_path !=  NULL) {
                ProspectPath = strdup( base_path );
        } else {        
                char *ptr = (char* ) getenv("PROSPECTPATH");
                if (ptr != NULL) {
                        ProspectPath = strdup(ptr);
                } else {
                        for (int i = 0; i < 4 && ProspectPath == NULL; i++) {
                                char buffer[5000];
                                sprintf(buffer, "%s/%s", search_paths[i], SINGLETON_FILE );
                                FILE *tmpfile = fopen( buffer, "r" );
                                if ( tmpfile ) {
                                        fclose(tmpfile);
                                        ProspectPath = strdup( search_paths[i] );
                                }
                        }
                }
        }
        template_paths = Template_Paths();      
        Init();
}


void ProspectParam::CopyGonnetMatrix( class TemplateStruct *template_data ) {
        template_data->CopyScoreTable( GONNET_MATRIX );
}



char * ProspectParam::TemplateFind(char *template_name) {
        long i = 0;
        FILE *file_test = NULL;
        char file_path[2000];
        
        while ( template_paths[i] && file_test == NULL) {
                //determine the path of the template file
                sprintf(file_path,"%s/%s.xml", template_paths[i], template_name);
                //doc = xmlParseFile(file_path);
                file_test = fopen( file_path, "r" );
                i++;
    }
        
        if (file_test == NULL) {
                return NULL;
        }
        fclose(file_test);
        return strdup(file_path);
}




char ** ProspectParam::Template_Paths(void) {
        FILE *fin;
        char str[1024];
        char outstr[1024], tmp[1024], found;
        char **out_strings, *curptr, *ptr, *ptr2;
        short path_count = 0;
        
        sprintf(str, "%s/template_paths", ProspectPath);
        fin = fopen(str, "r");
        
        path_count = 0;
        
        //check the TEMPLATE_PATH enviromental varible
        
        if (getenv("TEMPLATE_PATH")) {
                ptr = strdup(getenv("TEMPLATE_PATH"));
                curptr = strtok(ptr, ":");
                while(curptr) {
                        path_count++;
                        curptr = strtok(NULL, ":");
                }
                free(ptr);
        }
        
        
        while ( fgets(str, sizeof(str), fin)) {
                found = 0;  
                for (unsigned int i = 0; i < strlen(str); i++) {
                        if (isalnum(str[i]))
                                found = 1;
                }
                if (found)
                        path_count++;
        }
        
        out_strings = (char **)malloc(sizeof(char *) * (path_count+1) ); 
        
        rewind(fin);
        path_count = 0;
        
        if (getenv("TEMPLATE_PATH")) {
                ptr = strdup(getenv("TEMPLATE_PATH"));
                curptr = strtok(ptr, ":");
                while(curptr) {
                        out_strings[path_count++] = strdup(curptr);
                        curptr = strtok(NULL, ":");
                }
                free(ptr);
        }
        
        while ( fgets(str, sizeof(str), fin)) {
                found = 0;  
                for (unsigned int i = 0; i < strlen(str); i++) {
                        if (isalnum(str[i]))
                                found = 1;
                }
                if (found) {
                        for (unsigned int i = 0; i < strlen(str); i++) {
                                if (str[i] == '\n' || str[i] == '\r') {
                                        str[i] = 0;
                                }
                        }
                        //if there is an enviromental varible, translate it
                        memset(outstr, 0, sizeof(outstr));
                        curptr = str;
                        while ( *curptr != 0 && (ptr = strstr( curptr, "$(")) ) {
                                strncat(outstr, curptr, ptr-curptr);
                                ptr += 2;
                                ptr2 = ptr;
                                if ( (ptr2 = strstr( ptr, ")")) ) {
                                        memset(tmp, 0, sizeof(tmp));
                                        strncpy(tmp, ptr, ptr2 - ptr);
                                        if (getenv( tmp )) {
                                                strcat(outstr, getenv(tmp));
                                        }
                                        ptr2++;
                                }
                                curptr = ptr2;
                        }
                        strcat(outstr, curptr);
                        out_strings[path_count++] = strdup( outstr );
                }
        }
        out_strings[path_count++] = NULL;
        fclose(fin);
        
        return out_strings;
}





//#define DEBUG 1 


#define kwAminoAcid     "AMINOACID"
#define kwTorsionAngle  "TORSIONANGLE"
#define kwTorsion       "TORSION"
#define kwATOM          "ATOM"


ResidueInfoElement * Prospect_ResidueInfo_LoadFile(ProspectParam *param, char *path) {
        FILE *input;
        ResidueInfoElement *residues;
        char buffer[4048], buffer2[4048], tmp[100], res_name[25], atom_name[25], atom_type[25],
                res_char;
        int atom_count, res_num;
        float atom_loc[3];
        char *cur_ptr;
        
        input = fopen(path, "r");
        if ( input == NULL)
                return NULL;
        
        residues = (ResidueInfoElement *)malloc(sizeof(ResidueInfoElement) * 20);
        
        for (int i = 0; i < 20; i++) {
                residues[ i ].structure.atom_count = 0;
                residues[ i ].torsion_count = 0;
                residues[ i ].torsion_a = 0;
                residues[ i ].torsion_b = 0;
        }
        
        while (fgets(buffer, sizeof(buffer), input)) {
                if ( !strncmp(buffer, kwAminoAcid, strlen(kwAminoAcid)) ) {
                        sscanf(buffer, "%s %s %c %d", tmp,
                                   res_name, &res_char, &atom_count);
                        res_num = AA1toAANum( res_char );
                        if (res_num >= 0 && res_num < 20) {
                                residues[ res_num ].RS = res_char;
                                residues[ res_num ].structure.atom_count = 0;
                        } else {
                                //                  fprintf(stderr, "Formatting error in residue.lib: %s\n", buffer); 
                        }
                } else if ( !strncmp(buffer, kwATOM, strlen(kwATOM)) ) {
                        if (res_num >= 0 && res_num < 20) {
                                sscanf(buffer, "%s %s %s %f %f %f", 
                                           tmp, atom_name, atom_type, 
                                           &atom_loc[0], &atom_loc[1], &atom_loc[2] ); 
                                int atom_namenum = PDBAtomtoNum( atom_name );
                                if ( atom_namenum != -1) {
                                        residues[res_num].structure.atom[ residues[res_num].structure.atom_count ][0] = atom_loc[0];
                                        residues[res_num].structure.atom[ residues[res_num].structure.atom_count ][1] = atom_loc[1];
                                        residues[res_num].structure.atom[ residues[res_num].structure.atom_count ][2] = atom_loc[2];
                                        residues[res_num].structure.atom_name[ residues[res_num].structure.atom_count ] = atom_namenum;
                                        /*
                                         for (int i = 0; i < param->atom_data_count; i++) {
                                                 if ( atom_namenum == PDBAtomtoNum( param->atom_data[i].name) ) {
                                                         residues[res_num].structure.radius[residues[res_num].structure.atom_count] = param->atom_data[i].radius;
                                                 }
                                         }
                                         */
                                        residues[res_num].structure.radius[residues[res_num].structure.atom_count] = param->atom_data->FindRadius( atom_namenum );
                                        if ( atom_namenum == anN || atom_namenum == anCa || 
                                                 atom_namenum == anC || atom_namenum == anO  ) {
                                                residues[res_num].structure.bb_num[ atom_namenum ] = residues[res_num].structure.atom_count;
                                        }                                                
                                        residues[res_num].structure.atom_count++;
                                }
                                assert( residues[ res_num ].structure.atom_count < kMaxResAtoms );
                        }
                } else if ( !strncmp(buffer, kwTorsionAngle, strlen(kwTorsionAngle)) ) {
                        if (res_num >= 0 && res_num < 20) {
                                sscanf(buffer, "%s %d:", tmp, &residues[ res_num ].torsion_count);
                                if ( residues[ res_num ].torsion_count > 0 ) {
                                        residues[ res_num ].torsion_val = (pValType *)malloc(sizeof(pValType) * residues[ res_num ].torsion_count);
                                        residues[ res_num ].torsion_a = (char *)malloc(sizeof(char) * residues[ res_num ].torsion_count);
                                        residues[ res_num ].torsion_b = (char *)malloc(sizeof(char) * residues[ res_num ].torsion_count);
                                        
                                        strcpy(buffer2, buffer);
                                        cur_ptr = strtok( buffer2, " \t\n");
                                        cur_ptr = strtok( NULL, " \t\n:");
                                        cur_ptr = strtok( NULL, " \t\n");
                                        for (int i = 0; i < residues[ res_num ].torsion_count && cur_ptr; i++) {
                                                residues[ res_num ].torsion_val[i] = atof(cur_ptr);
                                                cur_ptr = strtok( NULL, " \t\n");
                                        }
                                        residues[ res_num ].torsion_mask = (char **)malloc(sizeof(char *) * residues[ res_num ].torsion_count);
                                        residues[ res_num ].torsion_mask[0] = (char *)malloc(sizeof(char) * residues[ res_num ].torsion_count * residues[ res_num ].structure.atom_count);
                                        for (int i = 0; i < residues[ res_num ].torsion_count; i++) {
                                                residues[ res_num ].torsion_mask[i] = residues[ res_num ].torsion_mask[0] + i * residues[ res_num ].structure.atom_count;
                                        }
                                        memset( residues[ res_num ].torsion_mask[0], 0, residues[ res_num ].torsion_count * residues[ res_num ].structure.atom_count );
                                } else {
                                        residues[ res_num ].torsion_val = NULL;
                                        residues[ res_num ].torsion_a = NULL;
                                        residues[ res_num ].torsion_b = NULL;
                                        residues[ res_num ].torsion_mask = NULL;
                                }
                        }
                } else if ( !strncmp(buffer, kwTorsion, strlen(kwTorsion)) ) {
                        if (res_num >= 0 && res_num < 20) {
                                int atom_count, torsion_num;
                                sscanf(buffer, "%s %d: %d", tmp, &torsion_num, &atom_count );
                                if ( torsion_num < residues[ res_num ].torsion_count ) {
                                        strcpy(buffer2, buffer);
                                        cur_ptr = strtok( buffer2, " \t\n\r");
                                        cur_ptr = strtok( NULL, " \t\n\r:");
                                        cur_ptr = strtok( NULL, " \t\n\r");
                                        cur_ptr = strtok( NULL, " \t\n\r");
                                        for (int i = 0; i < atom_count && cur_ptr; i++) {
                                                for (int j = 0; j < residues[ res_num ].structure.atom_count; j++) {
                                                        if ( residues[ res_num ].structure.atom_name[j] == PDBAtomtoNum( cur_ptr ) ) {
                                                                if (i == 0) {
                                                                        residues[res_num].torsion_a[ torsion_num ] = j;
                                                                } else if (i == 1) {
                                                                        residues[res_num].torsion_b[ torsion_num ] = j;
                                                                } else {
                                                                        residues[ res_num ].torsion_mask[ torsion_num ][j] = 1;
                                                                }
                                                        }
                                                }
                                                cur_ptr = strtok( NULL, " \t\n\r");
                                        }                       
                                } else {
                                        fprintf(stderr, "format error in rotamer.lib: %s\n", buffer);
                                }
                        }
                }
        }
        fclose(input);
        return residues;
}



RotamerArray * Prospect_RotamerArray_LoadBbind(char *path) {
        FILE *input;
        char buffer[4000], res_type[100], res_char, res_num, *cur_ptr;
        RotamerArray *rotamer_array;
        int cur_rotamer[20], r[4];
        
        input = fopen(path, "r");
        if (input == NULL)
                return NULL;
        
        rotamer_array = (RotamerArray *)malloc(sizeof(RotamerArray) * 20);
        for (int i = 0; i < 20; i++) {
                rotamer_array[i].rotamer_count = 0;
                rotamer_array[i].psi_bin = 720;
                rotamer_array[i].phi_bin = 720;
        }
        while (fgets( buffer, sizeof(buffer), input)) {
                if ( 5 == sscanf(buffer, "%s %d %d %d %d", res_type, &r[0], &r[1], &r[2], &r[3]) ) {
                        res_char = AA3toAA1( res_type );
                        res_num = AA1toAANum( res_char );                       
                        /*
                         if (rotamer_array[res_num].rotamer_count == 0) {
                                 for (int j = 0; j < 4; j++) {
                                         if ( r[j] != 0) 
                                                 rotamer_array[res_num].chi_count = j+1;
                                 }
                         }
                         */
                        rotamer_array[ res_num ].rotamer_count++;
                        rotamer_array[ res_num ].RS = res_char;
                }               
        }
        fseek(input, 0, SEEK_SET);
        
        for (int i = 0; i < 20; i++) {
                if (  rotamer_array[i].rotamer_count > 0 ) 
                        rotamer_array[i].rotamer = (RotamerElement *)malloc(sizeof(RotamerElement) * rotamer_array[i].rotamer_count );
                else 
                        rotamer_array[i].rotamer = NULL;
                cur_rotamer[i] = 0;
        }       
        
        while (fgets( buffer, sizeof(buffer), input)) {
                if ( 5 == sscanf(buffer, "%s %d %d %d %d", res_type, &r[0], &r[1], &r[2], &r[3]) ) {
                        res_char = AA3toAA1( res_type );
                        res_num = AA1toAANum( res_char );
                        cur_ptr = strtok(buffer, " \t\n");
                        for (int i = 0; i < 9; i++) {
                                cur_ptr = strtok(NULL, " \t\n");
                        }
                        rotamer_array[(int)res_num].rotamer[ cur_rotamer[ (int)res_num ] ].prob = atof( cur_ptr );
                        cur_ptr = strtok(NULL, " \t\n");
                        cur_ptr = strtok(NULL, " \t\n");                        
                        //                      for (int i = 0;  i < rotamer_array[res_num].chi_count && cur_ptr; i++) {
                        for (int i = 0; cur_ptr != NULL && i < 4; i++) {
                                rotamer_array[(int)res_num].rotamer[ cur_rotamer[(int)res_num ] ].chi[i] = atof(cur_ptr);
                                cur_ptr = strtok(NULL, " \t\n");
                                cur_ptr = strtok(NULL, " \t\n");
                        }
                        cur_rotamer[ (int)res_num ]++;
                        }
                }
        fclose(input);  
        return rotamer_array;
        }






int ProspectParam::InitRotamer(void) {
        char buffer[4000];
        
        atom_data = new AtomData( this );
        
        sprintf(buffer, "%s/%s", ProspectPath, RESIDUELIB_FILE);        
        residue_info = Prospect_ResidueInfo_LoadFile(this, buffer);
        if (residue_info == NULL)
                return 1;
        sprintf(buffer, "%s/%s", ProspectPath, BBIND_FILE);
        rotamer_array = Prospect_RotamerArray_LoadBbind(buffer);        
        if (rotamer_array == NULL) 
                return 1;
        
        return 0;
}

AtomData::AtomData(ProspectParam *in_parent) : ParamElement(in_parent) {
        char buffer[4000];
        
        sprintf(buffer, "%s/%s", parent->ProspectPath, ATOMDATA_FILE);  
        
        atom_data_count = 0;
        FILE *atom_file = fopen(buffer, "r");
        if (!atom_file)
                return;
        char tmp[100];
        float tmp_float;
        while (fgets(buffer, sizeof(buffer), atom_file)) {
                if ( 2 == sscanf(buffer, "%s %f", tmp, &tmp_float)) {
                        atom_data_count++;
                }               
        }
        fseek(atom_file, 0,  SEEK_SET);
        atoms = (AtomDataElement *)malloc(sizeof(AtomDataElement) * atom_data_count);
        int i = 0;
        while (fgets(buffer, sizeof(buffer), atom_file)) {
                if ( 2 == sscanf(buffer, "%s %f", tmp, &tmp_float)) {
                        strncpy(atoms[i].name, tmp, 4);
                        atoms[i].radius = tmp_float;
                        i++;
                }               
        }
        fclose(atom_file);              
}


float AtomData::FindRadius( char *name ) {
        return FindRadius( PDBAtomtoNum( name ) );
}

float AtomData::FindRadius( int name_num ) {
        for (int i = 0; i < atom_data_count; i++) {
                if ( name_num == PDBAtomtoNum( atoms[i].name) ) {
                        return atoms[i].radius;
                }
        }
        return -1;
}





float * get_float_file(char *path, long *count) {
        FILE *file;
        long i;
        float temp, *returnval;
        char templine[300];
        
        file = fopen(path, "r");
        if (file == NULL) {
                fprintf(stderr, "Unable to open: %s\n", path);
                return NULL;
        }
        *count = 0;
        while ( fgets( templine, sizeof(templine), file) ){
                if (strlen(templine) != 0) {
                        if ( sscanf(templine, "%f", &temp) ) {
                                (*count)++;
                        }
                }
        }       
        returnval = (float *)malloc(sizeof(float) * (*count) );
        fseek(file, 0, SEEK_SET);
        i = 0;
        while ( fgets( templine, sizeof(templine), file) ){
                if (strlen(templine) != 0) {
                        if ( sscanf(templine, "%f", &temp) ) {
                                returnval[i] = temp;
                                i++;
                        }
                }
        }
        fclose(file);   
        return returnval;
}



int ProspectParam::Init(void) {
        pValType val1, val2;
        int i,j, k;
        char command[2048];
        char dummy[1024], tmp1[1024], tmp2[1024];
        char dummyc, *ptr;
        
        second_guess = 0;
        filter = 0;
        no_loop = 0;
        v_level = 0;
        max_bytes = -1;
        output_timeinfo = 0;
        output_weights = 0;
        output_pdb = 0;
        output_backbone = 0;
        output_template_loc = 1;
        no_loop = 0;
        z_cycles = 1000;
        simple_twobody = 0;
        rotamer_array = NULL;
        residue_info = NULL;
        atom_data = NULL;
        atom_data_count = 0;
        ncpus = 1;
        full_twobody = false;
        //  READING IN FITNESS POTENTIAL:  SINGLETON  //
        
        FILE *file_in;
        sprintf(command,"%s/%s",ProspectPath, SINGLETON_FILE);
        file_in = fopen(command, "r"); 
        if (file_in == NULL) {
                //    sprintf(dummy, "Unable to open file %s", command);
                //    Prospect_Error_Exit(conf, dummy);
                return 1;
        }
        char residue[200];
        
        for (i = 0; i < 21; i++) {
                for (j = 0; j < 9; j++) {
                        SINGLETON[i][j] = 0;
                }
        }
        
        for (i=0; i < 20;i++) {  
                fscanf(file_in,"%s",residue);
                for (j=0; j < 9;j++) {
                        //     fscanf(file_in,"%lf",&conf->param.SINGLETON[i][j]);
                        fscanf(file_in,"%f",&SINGLETON[i][j]);
                        
                }
        }
        fclose(file_in);
        
        //  READING IN PAIR POTENTIAL:  PAIRWISE  //
        
        sprintf(command,"%s/%s",ProspectPath,PAIRWISE_FILE);
        file_in = fopen(command, "r");
        if (file_in == NULL) {
                //      sprintf(dummy, "Unable to open file %s", command);
                //       Prospect_Error_Exit(conf, dummy);
                return 1;
        }
        
        PAIR_POT = (pValType **)malloc(sizeof(pValType *) * 21);
#ifdef ALTIVEC
        PAIR_POT[0] = (pValType *)valloc(sizeof(pValType) * 21*21);
#else
        PAIR_POT[0] = (pValType *)malloc(sizeof(pValType) * 21*21);
#endif
        
        for (i = 0; i < 21; i++) {
                PAIR_POT[i] = PAIR_POT[0] + i*21;
                for (j = 0; j < 21; j++)
                        PAIR_POT[i][j] = 0;
        }
        
        for (i=0; i < 20; i++) {  
                fscanf(file_in,"%s",residue);
                for (j=0; j <= i; j++) {
                        float tmp_float;
                        fscanf(file_in, "%f", &tmp_float);
                        PAIR_POT[j][i]=PAIR_POT[i][j]=tmp_float;
                }
        }
        fclose(file_in);
        
        
        //  READING IN Contact count matrix  //
        for (int i=0; i < 21; i++) {
                for (int j = 0; j < MaxContactCount; j++) {
                        CONTACT_MATRIX[i][j] = 3200;
                }
        }
        sprintf(command,"%s/%s",ProspectPath,CONTACT_FILE);
        file_in = fopen(command, "r");
        for (int i=0; i < 20; i++) {  
                char buffer[2000];
                fgets(buffer, sizeof(buffer), file_in);
                char *str = strtok( buffer, "\t \n" );
                int j = 0;
                while ( (str = strtok( NULL, "\t \n" )) ) {
                        CONTACT_MATRIX[i][j] = atof( str );
                        j++;
                }
        }
        fclose(file_in);
        
        //  READING IN alignment matrix  //
        
        for (i=0; i < 21; i++) {
                for (j = 0; j < 21; j++) {
                        ALIGN_MATRIX[i][j] = 0;
                }
        }  
        
        sprintf(command,"%s/%s",ProspectPath,ALIGNMENT_FILE);
        file_in = fopen(command, "r");
        if (file_in == NULL) {
                //      sprintf(dummy, "Unable to open file %s", command);
                //       Prospect_Error_Exit(conf, dummy);
                return 1;
        }
        fgets(dummy,2000,file_in);
        for (i=0;i<20;i++) {
                dummyc=fgetc(file_in);
                for (j=0;j<20;j++) {
                        ALIGN_MATRIX[i][j]=fgetc(file_in);
                }
                fgets(dummy,2000,file_in);
        }
        fclose(file_in);
        
        
        //Reading in the GONNET mutation matrix
        sprintf(command, "%s/%s", ProspectPath, GONNETMATRIX_FILE);
        file_in = fopen(command, "r");
        
        for (i = 0; i < 21; i++) {
                for (j = 0; j < 21; j++) {
                        GONNET_MATRIX[i][j] = 0;
                }
        }
        
        while (fgets(dummy, sizeof(dummy), file_in)) {
                ptr = strtok(dummy, "\t \n");
                if (strcmp(ptr, "REM") ) {
                        int cur_aa = AA1toAANum( AA3toAA1( ptr ) );     
                        int i = 0;
                        while ( NULL != (ptr = strtok(NULL, "\t \n") ) ) {
                                GONNET_MATRIX[i][cur_aa] = -(atof( ptr ));
                                GONNET_MATRIX[cur_aa][i] = -(atof( ptr ));
                                i++;
                        }
                }
        } 
        fclose(file_in);  
        
        
        
        //Reading in the defire pairwise potential information
        
        //read the bin info
        
        
        j = -1;
        sprintf(command, "%s/%s", ProspectPath, DEFIRE_BIN_DATA);
        file_in = fopen(command, "r");
        while ( fgets(dummy, sizeof(dummy), file_in) ) {
                sscanf(dummy, "%d", &i);
                if (i > j) {
                        j = i;
                }
        }
        fseek(file_in, SEEK_SET, 0);
        dfire.bin_count = j;
        dfire.bin_min = (pValType*)malloc(sizeof(pValType) * j );
        dfire.bin_max = (pValType*)malloc(sizeof(pValType) * j );
        dfire.dist_max = 0;
        
        i = 0;
        while ( fgets(dummy, sizeof(dummy), file_in) ) {
                sscanf(dummy, "%d %f - %f", &i, &val1, &val2);
                dfire.bin_min[i-1] = val1;
                dfire.bin_max[i-1] = val2;
                if ( val2 > dfire.dist_max)
                        dfire.dist_max = val2;
        }
        fclose(file_in);
        
        dfire.dist_hash_size = (long)(dfire.dist_max * 10.0);
        dfire.dist_hash = (short *)malloc(sizeof(short) * dfire.dist_hash_size);
        for (i = 0; i < dfire.dist_hash_size; i++) 
                dfire.dist_hash[i] = -1;
        for (i = 0; i < dfire.bin_count; i++) {
                for (val1 = dfire.bin_min[i]; val1 < dfire.bin_max[i]; val1 += .1) {
                        dfire.dist_hash[ (long)(val1 * 10.0) ] = i;
                }
        }
        
        dfire.bin = (pValType ***)malloc(sizeof(pValType**) * dfire.bin_count);
        dfire.bin[0] = (pValType **)malloc(sizeof(pValType*) * dfire.bin_count * 21);
        dfire.bin[0][0] = (pValType *)malloc(sizeof(pValType) * 21 * 21 * dfire.bin_count);
        for (i = 0; i < dfire.bin_count; i++) {
                dfire.bin[i] = dfire.bin[0] + (i * 21);
                for (j = 0; j < 21; j++) {
                        dfire.bin[i][j] = dfire.bin[0][0] + (i * 21 * 21) + j * 21; 
                        for (k = 0; k < 21; k++) {
                                dfire.bin[i][j][k] = 0;
                        }
                }
        }
        
        sprintf(command, "%s/%s", ProspectPath, DEFIRE_CB_DATA);
        file_in = fopen(command, "r");
        
        while (fgets( dummy, sizeof(dummy), file_in) ) {
                sscanf(dummy, "%s %s %f %d", tmp1, tmp2, &val1, &i);
                j = AA1toAANum( AA3toAA1( tmp1 ) );
                k = AA1toAANum( AA3toAA1( tmp2 ) );
                dfire.bin[i-1][j][k] = val1;
                dfire.bin[i-1][k][j] = val1;    
        }
        fclose(file_in);
        
        
        
        //Read in information for secondary structure prediction
    /*
         FILE *file;
         long i;
         long temp;
         char s[2048];
         
         len_weight = the_len_weight;
         trans_weight = the_trans_weight;
         cap_weight = the_cap_weight;
         
         sprintf(s, "%s/data/parameters/pred/c_len.dist", ProspectPath );
         c_hist = get_float_file(s, &c_hist_count);
         sprintf(s, "%s/data/parameters/pred/e_len.dist", ProspectPath );
         e_hist = get_float_file(s, &e_hist_count);
         sprintf(s, "%s/data/parameters/pred/h_len.dist", ProspectPath );
         h_hist = get_float_file(s, &h_hist_count);
         sprintf(s,"%s/data/parameters/pred/hc.dat" , ProspectPath );
         
         hc_aa_freq = get_float_file(s, &temp);
         //assert( temp == (20 * 4) );
         sprintf(s,"%s/data/parameters/pred/ch.dat" , ProspectPath );
         ch_aa_freq = get_float_file(s, &temp);
         //assert( temp == (20 * 4) );
         sprintf(s,"%s/data/parameters/pred/ce.dat" , ProspectPath );
         ce_aa_freq = get_float_file(s, &temp);
         //assert( temp == (20 * 4) );
         sprintf(s, "%s/data/parameters/pred/ec.dat", ProspectPath );
         ec_aa_freq = get_float_file(s, &temp);
         //assert( temp == (20 * 4) );  
         */
        
        
        return 0;
}






int ProspectParam::DoThreading(
                                                           char *threader_name,
                                                           WeightArray weight,           
                                                           TemplateStruct *template_data,
                                                           TargetStruct *target_data,
                                                           AlignmentStruct *alignment,
                                                           ScoreStruct *scores ) {
        threading_function threader = this->ThreaderFind( threader_name );
        if ( target_data->twobody == NULL )
                target_data->Setup( this );
        DoThreading(   threader,
                                   weight,               
                                   template_data,
                                   target_data,
                                   alignment,
                                   scores );
}


int ProspectParam::DoThreading(
                                                           threading_function threader,
                                                           WeightArray weight,           
                                                           TemplateStruct *template_data,
                                                           TargetStruct *target_data,
                                                           AlignmentStruct *alignment,
                                                           ScoreStruct *scores ) {
        
        
        //call threading function
        char * status = (*threader)(this, weight, template_data, target_data, alignment );       
        
        if (status) {
                return -1;
        }       
        if ( scores != NULL ) {
                if ( scores->mask & ScoreStruct::DO_SCORE_RAW )
                        scores->ScoreRaw( this,
                                                          weight,                
                                                          template_data,
                                                          target_data, 
                                                          alignment);           
                if ( scores->mask & ScoreStruct::DO_SCORE_Z || ( !(scores->mask & ScoreStruct::DO_SCORE_ZFULL) && (scores->mask & ScoreStruct::DO_SCORE_ZENERGY) ) )
                        scores->ScoreZ( this,
                                                        weight,                  
                                                        template_data,
                                                        target_data, 
                                                        alignment);             
                if ( scores->mask & ScoreStruct::DO_SCORE_ZFULL )
                        scores->ScoreZFull( this,
                                                                weight,                  
                                                                template_data,
                                                                target_data, 
                                                                alignment);
        }
        return 0;
}

---