src/common/score_functions.cc

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

#include <pthread.h>


#include "open_prospect.h"

#include "../dynamic_threading/dynamic_threading.h"



int count_nalign(  TargetStruct *target_data, AlignmentStruct *alignment ) {
        int nalign = 0;
        for (int i = 0; i < target_data->len; i++) {
                if (alignment->GetTargetAlign(i) != -1)
                        nalign++;
        }
        return nalign;  
}



ScoreStruct::ScoreStruct(void ) {
        mask = 0x00000000;      
}

ScoreStruct::~ScoreStruct(void ) {
                
}



char score_names[ScoreStruct::SCORE_COUNT_TOTAL][15] = {
        "raw",
        "z",
        "zfull",
        "nn",
        "user0",
        "user1",
        "user2",
        "user3",
        "zmutationlog",
        "zsingleton",
        "zsec_struct",
        "ztwobody",
        "zdfire"
};


char *ScoreStruct::GetScoreName(int num) {
        return score_names[num];        
}




pValType ScoreStruct::ScoreRaw( ProspectParam *param,
                                                          WeightArray weight,            
                                                          TemplateStruct *template_data,
                                                          TargetStruct *target_data, 
                                                          AlignmentStruct *alignment) {
        EnergyArray energy;
        alignment->GetEnergies(param, template_data, target_data, energy);      
        pValType val = 0;
        for (int i = 0; i< e_count; i++) {
                val += weight[i] * energy[i];
        }       
        this->values[ SCORE_RAW ] = val;
        this->mask |= (0x01 << SCORE_RAW);
        return val;
}

//--------------------------------------------------------------------
void shuffle_seq(TargetStruct *src, TargetStruct *dst) {
        char seq_mask[ src->len ];
        memset( seq_mask, 0x00, src->len );     
        for (int i=0; i < src->len;i++) {
                int j;
                do {
                        j = ((int)(((pValType)src->len * ((pValType)rand()) / ((pValType)RAND_MAX)))) % src->len;
                } while ( seq_mask[j] );        
#if 1
                dst->residue[i] = src->residue[j];
#else
                dst->residue[i].RS = src->residue[j].RS;
                memcpy( dst->residue[i].FREQ, src->residue[j].FREQ, sizeof(pValType) * 20 );    
#endif
                dst->residue[i].orig_num = j;
                seq_mask[j] = 1;
        }
}

typedef struct  {
        ProspectParam *param;
        WeightArray weight; 
        TemplateStruct *template_data;
        TargetStruct *target_data;
        EnergyArray *energy_array;
        AlignmentStruct *original_align;
        int step_start;
        int step_size;  
        bool full;
        long score_mask;
} z_loop_calc_data;


void *z_loop_calc_wrapper(void *arg) {
        z_loop_calc_data *data = (z_loop_calc_data*)arg;
        if ( data->full ) {
                TargetStruct tmp_target;                                                 
                tmp_target.len = data->target_data->len;
                tmp_target.distconst_count = 0;
                tmp_target.twobody = data->target_data->twobody; //this is a hack, original twobody table
                tmp_target.residue = (TargetResidue *)malloc(sizeof(TargetResidue) * data->target_data->len );
                for (int i = 0; i < data->target_data->len; i++) {
                        tmp_target.residue[i] = data->target_data->residue[i];
                }       
                for (int i = 0 + data->step_start; i < data->param->z_cycles; i += data->step_size) {
                        AlignmentStruct  new_alignment;
                        shuffle_seq(data->target_data, &tmp_target );
                        char *status = dynamic_thread(data->param, data->weight, data->template_data, &tmp_target, &new_alignment );
                        new_alignment.GetEnergies(data->param, data->template_data, &tmp_target, data->energy_array[i]);
                }       
                tmp_target.twobody = NULL; //undo the hack, so the twobody table isn't delallocated
        } else {
                TargetStruct tmp_target;
                tmp_target.len = data->target_data->len;
                tmp_target.twobody = data->target_data->twobody;        
                tmp_target.residue = (TargetResidue *)malloc(sizeof(TargetResidue) * data->target_data->len );  
                for (int i = 0; i < data->target_data->len; i++) {
                        tmp_target.residue[i] = data->target_data->residue[i];
                }
                for (int i = 0 + data->step_start; i < data->param->z_cycles; i += data->step_size) {
                        shuffle_seq(data->target_data, &tmp_target );
                        data->original_align->GetEnergies(data->param, data->template_data, &tmp_target, data->energy_array[i],  data->score_mask);
                }       
                tmp_target.twobody = NULL; //undo the hack, so the twobody table isn't delallocated
        }       
        if ( data->step_size > 1 )
                pthread_exit(NULL);
        return NULL;
}



pValType energy_array_stat( ProspectParam *param, WeightArray weight, 
                                                        TemplateStruct *template_data, TargetStruct *target_data,
                                                        AlignmentStruct *original_align,
                                                        EnergyArray *energy_array, ScoreStruct *scores, 
                                                        bool full = false, long score_mask = 0xFFFFFFFF ) {             
        double raw_thread, total, variance, sd, zscore, mean, raw_original;
        EnergyArray energy_total, energy_mean, energy_sd, energy_z, energy_original;

        
        original_align->GetEnergies(param, template_data, target_data, energy_original, score_mask);
        raw_original = 0;
        for (long i = 0; i < e_count; i++) {
                if ( e_types[i] != et_gap )
                        raw_original += (double)weight[i] * (double)energy_original[i];
        }
        
        double *raw_scores = (double *)malloc(sizeof(double) * param->z_cycles);
        for (long i = 0; i < param->z_cycles; i++) {
                raw_scores[i] = 0;
                for (long j = 0; j < e_count; j++) 
                        raw_scores[i] += (double)weight[j] * (double)energy_array[i][j];
        }
        
        total = 0;
        for (int i = 0; i < param->z_cycles; i++) {
                total += raw_scores[i];
        }
        mean = total / (pValType)param->z_cycles;
        total = 0;
        for (int i = 0; i < param->z_cycles; i++) {
                total += pow( raw_scores[i] - mean, 2);
        }
        variance = total / (pValType)param->z_cycles;
        sd = sqrt(variance);
        if ( sd == 0 )
                sd = 1;
        zscore = (mean - raw_thread)/sd;
        
        
        //total energy zscore   
        total = 0;
        for (long i = 0; i < param->z_cycles; i++) {
                total += (double)raw_scores[i];
        }
        mean = total / (double)param->z_cycles;
        total = 0;
        for (long i = 0; i < param->z_cycles; i++) {
                total += pow( raw_scores[i] - mean, 2);
        }
        variance = total / (double)param->z_cycles;
        sd = sqrt(variance);
        
        zscore = (mean - raw_original)/sd;      
        
        if (param->v_level > 0) {
                printf("raw: %lf  mean: %lf  sd: %lf  zscore: %lf\n", 
                           raw_thread, mean, sd, zscore );
        }       
        //      assert( !isnan(zscore) );
        if ( isnan(zscore) ) {
                zscore = 0;
        }
        free(raw_scores);               
        if ( full ) {
                scores->values[ ScoreStruct::SCORE_ZFULL ] = zscore;
                scores->z_sd[ ScoreStruct::SCORE_ZFULL ] = sd; scores->z_mean[ ScoreStruct::SCORE_ZFULL ] = mean; 
        } else {
                scores->values[ ScoreStruct::SCORE_Z ] = zscore;
                scores->z_sd[ ScoreStruct::SCORE_Z ] = sd; scores->z_mean[ ScoreStruct::SCORE_Z ] = mean; 
        }
        
        
        for (int j = 0; j < e_count; j++) {
                energy_total[j] = 0;
        }
        
        for (long i = 0; i < param->z_cycles; i++) {
                for (int j = 0; j < e_count; j++) {
                        energy_total[j] += energy_array[i][j];
                }
        }
        for (int i = 0; i < e_count; i++) {
                energy_mean[i] = energy_total[i] / (pValType)  param->z_cycles;
        }
        for (int j = 0; j < e_count; j++) {
                energy_total[j] = 0;
        }
        for (long i = 0; i < param->z_cycles; i++) {
                for (int j = 0; j < e_count; j++) {
                        energy_total[j] += pow(energy_array[i][j] - energy_mean[j], 2);
                        //      assert( !isnan(energy_total[j]) );
                }
        }
        for (int i = 0; i < e_count; i++) {
                energy_sd[i] = sqrt( energy_total[i] / (pValType)  param->z_cycles );
        }
        
        for (int i = 0; i < e_count; i++) {
                energy_z[i] = (energy_mean[i] - energy_original[i]) / energy_sd[i];
        }
        
        scores->values[ ScoreStruct::SCORE_ZMUTATIONLOG ] = energy_z[ e_mutationlog ];
        scores->z_sd[ ScoreStruct::SCORE_ZMUTATIONLOG ] = energy_sd[ e_mutationlog ]; scores->z_mean[ ScoreStruct::SCORE_ZMUTATIONLOG ] = energy_mean[ e_mutationlog ]; 
        scores->values[ ScoreStruct::SCORE_ZSINGLETON ] = energy_z[ e_singleton ];
        scores->z_sd[ ScoreStruct::SCORE_ZSINGLETON ] = energy_sd[ e_singleton ]; scores->z_mean[ ScoreStruct::SCORE_ZSINGLETON ] = energy_mean[ e_singleton ]; 
        scores->values[ ScoreStruct::SCORE_ZSEC_STRUCT ] = energy_z[ e_sec_struct ];
        scores->z_sd[ ScoreStruct::SCORE_ZSEC_STRUCT ] = energy_sd[ e_sec_struct ]; scores->z_mean[ ScoreStruct::SCORE_ZSEC_STRUCT ] = energy_mean[ e_sec_struct ]; 
        scores->values[ ScoreStruct::SCORE_ZTWOBODY ] = energy_z[ e_twobody ];
        scores->z_sd[ ScoreStruct::SCORE_ZTWOBODY ] = energy_sd[ e_twobody ]; scores->z_mean[ ScoreStruct::SCORE_ZTWOBODY ] = energy_mean[ e_twobody ]; 
        scores->values[ ScoreStruct::SCORE_ZDFIRE ] = energy_z[ e_dfire ];
        scores->z_sd[ ScoreStruct::SCORE_ZDFIRE ] = energy_sd[ e_dfire ]; scores->z_mean[ ScoreStruct::SCORE_ZDFIRE ] = energy_mean[ e_dfire ]; 
        
        for (int i = ScoreStruct::SCORE_ZENERGY_START; i < ScoreStruct::SCORE_ZENERGY_END; i++) {
                if ( isnan( scores->values[i] ) ) {
                        scores->values[i] = 0.0;
                }
        }               
        return zscore;
}



pValType ScoreStruct::ScoreZ(ProspectParam *param,
                                                 WeightArray weight,             
                                                 TemplateStruct *template_data, 
                                                 TargetStruct *target_data,
                                                 AlignmentStruct *original_align) {
        long score_mask = (e_mask( e_mutationlog ) | e_mask( e_singleton ) | e_mask( e_sec_struct ) | e_mask( e_twobody ) | e_mask( e_dfire ) );
        EnergyArray  *energy_array = (EnergyArray *)malloc(sizeof(EnergyArray) * param->z_cycles);
        z_loop_calc_data *zloop_data = (z_loop_calc_data *)malloc(sizeof(z_loop_calc_data) * param->ncpus );
        
        for (int i = 0; i < param->ncpus; i++) {        
                for (int j= 0; j < e_count; j++)
                        zloop_data[i].weight[j] = weight[j];
                zloop_data[i].param = param;
                zloop_data[i].template_data = template_data; 
                zloop_data[i].target_data = target_data;
                zloop_data[i].energy_array = energy_array;
                zloop_data[i].step_size = param->ncpus;                         
                zloop_data[i].original_align = original_align;
                zloop_data[i].step_start = i;
                zloop_data[i].full = false;
                zloop_data[i].score_mask = score_mask;
        }       
        if ( param->ncpus > 1 ) {               
                pthread_t *threads;
                threads = (pthread_t *)malloc(sizeof(pthread_t) * param->ncpus);
                for (int i = 0; i < param->ncpus; i++) {        
                        pthread_create(&threads[i], NULL, z_loop_calc_wrapper, &zloop_data[i]);
                }
                for (int i = 0; i < param->ncpus; i++) {        
                        pthread_join(threads[i], NULL);
                }
                free(threads);
                free(zloop_data);               
        } else {
                z_loop_calc_wrapper( zloop_data );
        }       
        pValType zscore = energy_array_stat( param, weight, 
                                                                                 template_data, target_data,
                                                                                 original_align,
                                                                                 energy_array, this, 
                                                                                 false, score_mask );
        free( energy_array );   
        this->values[ SCORE_Z ] = zscore;
        this->mask |= (0x01 << SCORE_Z) ;

        return zscore;
}


pValType ScoreStruct::ScoreZFull(ProspectParam *param,
                                                                 WeightArray weight,             
                                                                 TemplateStruct *template_data, 
                                                                 TargetStruct *target_data,
                                                                 AlignmentStruct *original_align) {
        EnergyArray  *energy_array = (EnergyArray *)malloc(sizeof(EnergyArray) * param->z_cycles);
        z_loop_calc_data *zloop_data = (z_loop_calc_data *)malloc(sizeof(z_loop_calc_data) * param->ncpus );
        for (int i = 0; i < param->ncpus; i++) {        
                for (int j= 0; j < e_count; j++)
                        zloop_data[i].weight[j] = weight[j];
                zloop_data[i].param = param;
                zloop_data[i].template_data = template_data; 
                zloop_data[i].target_data = target_data;
                zloop_data[i].energy_array = energy_array;
                zloop_data[i].step_size = param->ncpus;                 
                zloop_data[i].original_align = original_align;
                zloop_data[i].step_start = i;
                zloop_data[i].full = true;              
                zloop_data[i].score_mask = 0xFFFFFFFF;
        }
        
        if ( param->ncpus > 1 ) {               
                pthread_t *threads;
                threads = (pthread_t *)malloc(sizeof(pthread_t) * param->ncpus);
                for (int i = 0; i < param->ncpus; i++) {        
                        pthread_create(&threads[i], NULL, z_loop_calc_wrapper, &zloop_data[i]);
                }
                for (int i = 0; i < param->ncpus; i++) {        
                        pthread_join(threads[i], NULL);
                }
                free(threads);
                free(zloop_data);               
        } else {
                z_loop_calc_wrapper( zloop_data );
        }
        pValType zscore = energy_array_stat( param, weight, 
                                                                                 template_data, target_data,
                                                                                 original_align,
                                                                                 energy_array, this, 
                                                                                 true, 0xFFFFFFFF );
        free( energy_array );   
        this->values[ SCORE_ZFULL ] = zscore;
        this->mask |= (0x01 << SCORE_ZFULL) ;

        return zscore;
}

---