src/common/dynamic_align.cc

#include "open_prospect.h"


pValType **_2d_array_malloc( int x, int y ) {
        pValType **array;
        
        array = (pValType **)malloc(sizeof(pValType *) * x);
        array[0] = (pValType *)malloc(sizeof(pValType) * x * y );
        
        for ( int i = 0; i < x; i++ ) {
                array[i] = array[0] + i * y;
                for (int j = 0; j < y; j++) 
                        array[i][j] = 0.0;
        }
        return array;
}

void _2d_array_free( pValType **array ) {
        free( array[0] );
        free( array );  
}


void DynamicAlign::SetDim( int x, int y ) {
        if ( energy_matrix ) {
                free( energy_matrix[0] );
                free( energy_matrix );
        }
        if ( score_matrix ) {
                free( score_matrix[0] );
                free( score_matrix );
        }
        if ( source_matrix ) {
                free( source_matrix[0] );
                free( source_matrix );
        }       
        size_x = x;
        size_y = y;     
        energy_matrix = (pValType **)malloc(sizeof(pValType*) * (x + 1) );
        energy_matrix[0] = (pValType *)malloc(sizeof(pValType) * (x + 1) * (y + 1) );
        
        score_matrix = (pValType **)malloc(sizeof(pValType*) * (x + 1) );
        score_matrix[0] = (pValType *)malloc(sizeof(pValType) * (x + 1) * (y + 1) );
        
        source_matrix = (char **)malloc(sizeof(char*) * (x + 1) );
        source_matrix[0] = (char *)malloc(sizeof(char) * (x + 1) * (y + 1) );
        
        for (int i = 0; i < size_x+1; i++) {
                score_matrix[i] = score_matrix[0] + (size_y+1) * i;
                source_matrix[i] = source_matrix[0] + (size_y+1) * i;
                energy_matrix[i] = energy_matrix[0] + (size_y+1) * i;
                for (int j = 0; j < size_y + 1; j++) {
                        score_matrix[i][j] = NAN;
                        energy_matrix[i][j] = NAN;
                        source_matrix[i][j] = 0;                
                }
        }               
}


void ProspectDynamicAlign::CalcEnergyTable( ProspectParam *param,
                                                                                        WeightArray weights,   
                                                                                        TemplateStruct *template_data,
                                                                                        TargetStruct *target_data )  {  
        SetDim( target_data->len, template_data->len ); 
        gap_open = weights[ e_gapopen ];
        gap_extend = weights[ e_gapconst ];
        for (int i = 0; i < target_data->len; i++) {
                for (int j = 0; j < template_data->len; j++) {
                        int template_residue_no = j;
                        int target_residue_no = i;
                        energy_matrix[i][j] = calc_single_energy(param, weights, template_data, target_data, template_residue_no, target_residue_no);
                        assert(!isnan(energy_matrix[i][j]));
                }
        }
}



void DynamicAlign::ContinueAlign( int start_x, int start_y ) {
        pValType val;
        char state;
        if ( start_x < size_x && start_y < size_y ) {
                if ( start_x > 0 && start_y > 0 ) {
                        val = GetCellVal( start_x-1, start_y-1 );
                        state = GetCellState( start_x-1, start_y-1 );
                        //assert( !isnan( val ) );
                } else {
                        if ( start_x == 0 && start_y == 0 ) {
                                state = vMatch;
                                val = 0;
                        } else {
                                if ( start_x > 0 ) {
                                        state = vLeft;
                                        val = gap_open + gap_extend * ( start_x );
                                }
                                if ( start_y > 0 ) {
                                        state = vUp;
                                        val = gap_open + gap_extend * ( start_y );
                                }
                        }
                }
                Align( start_x, start_y, size_x, size_y, state, val );
        }
}



void DynamicAlign::ForceXDelete( int start_x, int len ) {
        for ( int x = start_x; x < start_x + len; x++) {
                for ( int y = 0; y < size_y; y++ ) {
                        pValType score = score_matrix[ x ][ y+1 ];
                        if ( source_matrix[x][y+1] != vLeft )
                                score += gap_open;
                        score += gap_extend;
                        source_matrix[x+1][y+1] = vLeft;
                }
        }
}

void DynamicAlign::ForceYDelete( int start_y, int len ) {
        for ( int y = start_y; y < start_y + len; y++) {
                for ( int x = 0; x < size_x; x++ ) {
                        pValType score = score_matrix[ x+1 ][ y ];
                        if ( source_matrix[x+1][y] != vUp )
                                score += gap_open;
                        score += gap_extend;
                        source_matrix[x+1][y+1] = vUp;
                }
        }       
}

void DynamicAlign::ForceAlign( int start_x, int start_y, int len ) {
        for (int i = 0; i < len; i++) {
                if ( start_x+i < size_x && start_y+i < size_y &&
                         start_x+i >= 0 && start_y+i >= 0) {
                        score_matrix[start_x+i+1][start_y+i+1] = score_matrix[start_x+i][start_y+i] + energy_matrix[start_x+i][start_y+i];
                        source_matrix[start_x+i+1][start_y+i+1] = vMatch;
                        //assert( !isnan( score_matrix[start_x+i+1][start_y+i+1] ) );
                }
        }
}




void DynamicAlign::Align(int start_x, int start_y, int end_x, int end_y, char start_state, pValType start_val) {        
        int x,y;
        pValType score_left, score_up, score_match;
                
        assert(start_x >= 0 && start_x < size_x);
        assert(start_y >= 0 && start_y < size_y);
        assert(end_x >= 0 && end_x <= size_x);
        assert(end_y >= 0 && end_y <= size_y);
        assert(start_x <= end_x);
        assert(start_y <= end_y);       
        
        score_matrix[start_x][start_y] = start_val;
        source_matrix[start_x][start_y] = start_state;
        
        //fill in gap edges
        for (y = start_y+1; y < end_y+1; y++) {
                if ( start_state != vUp )
                        score_matrix[start_x][y] = start_val + (gap_open) + (gap_extend)*(y-(start_y));
                else 
                        score_matrix[start_x][y] = start_val + (gap_extend)*(y-(start_y));
                source_matrix[start_x][y] = vUp;
        }
        for (x = start_x+1; x < end_x+1; x++) {
                if ( start_state != vLeft )
                        score_matrix[x][start_y] = start_val + (gap_open) + (gap_extend)*(x-(start_x));
                else 
                        score_matrix[x][start_y] = start_val + (gap_extend)*(x-(start_x));
                source_matrix[x][start_y] = vLeft;
        }
        
        for (x = start_x; x < end_x; x++) {
                for (y = start_y; y < end_y; y++) {
                        if ( x >= 0  && x < size_x  ) {
                                score_matrix[x+1][y+1] = energy_matrix[x][y];
                                assert(!isnan( score_matrix[x+1][y+1] ) ); 
                        } else {
                                score_matrix[x+1][y+1] = 0;
                        }
                }
        }       
        for (x = start_x; x < end_x; x++) {
                for (y = start_y; y < end_y; y++) {
                        score_match = score_matrix[x][y] + score_matrix[x+1][y+1];
                        
                        score_up = score_matrix[x+1][y];
                        if (source_matrix[x+1][y] == vMatch || source_matrix[x+1][y] == vLeft ) 
                                score_up += gap_open;
                        score_up += gap_extend;
                        
                        score_left = score_matrix[x][y+1];
                        if (source_matrix[x][y+1] == vMatch || source_matrix[x][y+1] == vUp ) 
                                score_left += gap_open;
                        score_left += gap_extend;
                        
                        if ( score_match <= score_left && score_match <= score_up) {
                                score_matrix[x+1][y+1] = score_match;
                                source_matrix[x+1][y+1] = vMatch;
                        } else {
                                if ( score_left <= score_up ) {
                                        score_matrix[x+1][y+1] = score_left;
                                        source_matrix[x+1][y+1] = vLeft;
                                } else {
                                        score_matrix[x+1][y+1] = score_up;
                                        source_matrix[x+1][y+1] = vUp;
                                }
                        }
                        assert(!isnan(score_matrix[x+1][y+1]));
                }
        }       
}               

#if 0
//#ifdef ALTIVEC_test
vector float cur_score_vector;
vector unsigned char source_vec, tmp_vec1, tmp_vec2, tmp_vec3, vUp_vec, vMatch_vec;
float *score_start, *score_row_start, *score_ptr, *top_ptr, *left_ptr, *top_left_ptr;
char *source_start, *source_ptr, *source_row_start;
vector bool select_vec1, select_vec2, select_vec3, select_vec4;
vector float open_vec, extend_vec, 
score_vec1, score_vec2, score_vec3, score_vec4, 
top_vec1, top_vec2, top_vec3, top_vec4,
top_left_vec1, top_left_vec2, top_left_vec3, top_left_vec4,
left_vec1, left_vec2, left_vec3, left_vec4,
score_match_vec1, score_match_vec2, score_match_vec3, score_match_vec4,
score_up_vec1, score_up_vec2, score_up_vec3, score_up_vec4;
vector unsigned char vAlign,
vSourceAlign1 = (vector unsigned char)(0x00,0x04,0x08,0x0C,0x10,0x14,0x18,0x1C,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00),
vSourceAlign2 = (vector unsigned char)(0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x04,0x08,0x0C,0x10,0x14,0x18,0x1C),
vSourceAlign3 = (vector unsigned char)(0x00,0x01,0x02,0x03,0x04,0x05,0x06,0x07,0x18,0x19,0x20,0x21,0x22,0x23,0x24,0x25);
long row_size = (((template_data->len/16)+2)*16); 
unsigned char blksize=4; /* number of 16-byte chunks, 0-31 */
unsigned char blkcount=1; /* number of blks, 0=256 */
short blkstride = 64;         /* stride in bytes between blocks 0=+32768 */
int cwrd = blkstride | (blkcount<<16) | (blksize<<24);

vUp_vec = vec_splat_u8( vUp );
vMatch_vec = vec_splat_u8( vMatch ); 

// permute vector required to move it to element 0
open_vec = vec_lde( 0, &open_score );
vAlign = vec_lvsl( 0, &open_score );
open_vec = vec_perm( open_vec, open_vec, vAlign );
open_vec = vec_splat( open_vec, 0 );

extend_vec = vec_lde( 0, &extend_score );
vAlign = vec_lvsl( 0, &extend_score );
extend_vec = vec_perm( extend_vec, extend_vec, vAlign );
extend_vec = vec_splat( extend_vec, 0 );

score_matrix = (float **)malloc(sizeof(float *) * (target_data->len+1) );
score_start = (float *)valloc(sizeof(float) * (row_size) * (target_data->len+1) );
score_matrix[0] = &score_start[3];

source_matrix = (char **)malloc(sizeof(char *) * (target_data->len+1));
source_start = (char *)valloc(sizeof(char) * (row_size) * (target_data->len+1) );
source_matrix[0] = &source_start[15];
//  source_matrix[0] = (char *)memalign(16, sizeof(char) * (row_size) * (target_data->len+1) );

for (i = 1; i < target_data->len+1; i++) {
        score_matrix[i] = score_matrix[0] + (row_size) * i;
        source_matrix[i] = source_matrix[0] + (row_size) * i;
}

score_ptr = score_start + row_size + 4;
source_ptr = source_start + row_size + 16;
top_ptr = score_start + 4;
top_left_ptr = score_start + 3;

for (x = 0; x < target_data->len; x++) {
    score_row_start = score_ptr;
    source_row_start = source_ptr;
    for (y = 0; y < template_data->len; y+=16) {
                vec_dstst( score_ptr, cwrd, 0);
                vec_dstst( top_ptr, cwrd, 0);
                vec_dstst( (unsigned char *)source_ptr, cwrd, 0);
                
                score_vec1 = vec_ld(0, score_ptr);
                score_vec2 = vec_ld(4, score_ptr);
                score_vec3 = vec_ld(8, score_ptr);
                score_vec4 = vec_ld(12, score_ptr);
                
                source_vec = vec_ld(0, (unsigned char*)source_ptr);
                
                top_vec1 = vec_ld(0, top_ptr);
                top_vec2 = vec_ld(4, top_ptr);
                top_vec3 = vec_ld(8, top_ptr);
                top_vec4 = vec_ld(12, top_ptr);
                
                top_left_vec1 = vec_ld(0, top_left_ptr);
                vAlign = vec_lvsl ( 0 , top_left_ptr);   
                top_left_vec1 = vec_perm ( top_left_vec1, top_vec1, vAlign);
                top_left_vec2 = vec_perm ( top_vec1, top_vec2, vAlign );
                top_left_vec3 = vec_perm ( top_vec2, top_vec3, vAlign );
                top_left_vec4 = vec_perm ( top_vec3, top_vec4, vAlign );
                
                score_match_vec1 = vec_add(score_vec1, top_left_vec1);
                score_match_vec2 = vec_add(score_vec2, top_left_vec2);
                score_match_vec3 = vec_add(score_vec3, top_left_vec3);
                score_match_vec4 = vec_add(score_vec4, top_left_vec4);
                
                score_up_vec1 = vec_add(score_vec1, top_vec1);
                score_up_vec2 = vec_add(score_vec1, top_vec2);
                score_up_vec3 = vec_add(score_vec1, top_vec3);
                score_up_vec4 = vec_add(score_vec1, top_vec4);
                
                select_vec1 = vec_cmple( score_match_vec1, score_up_vec1 );
                select_vec2 = vec_cmple( score_match_vec2, score_up_vec2 );
                select_vec3 = vec_cmple( score_match_vec3, score_up_vec3 );
                select_vec4 = vec_cmple( score_match_vec4, score_up_vec4 );
                
                tmp_vec1 = vec_perm( (vector unsigned char)select_vec1, (vector unsigned char)select_vec2, vSourceAlign1);
                tmp_vec2 = vec_perm( (vector unsigned char)select_vec3, (vector unsigned char)select_vec4, vSourceAlign2);
                
                tmp_vec3 = vec_perm( tmp_vec1, tmp_vec2, vSourceAlign3);
                
                source_vec = vec_sel( vMatch_vec, vUp_vec, tmp_vec3 );
                
                score_vec1 = vec_sel( score_match_vec1, score_up_vec1, select_vec1 );
                score_vec2 = vec_sel( score_match_vec2, score_up_vec2, select_vec2 );
                score_vec3 = vec_sel( score_match_vec3, score_up_vec3, select_vec3 );
                score_vec4 = vec_sel( score_match_vec4, score_up_vec4, select_vec4 );
                
                vec_st(source_vec, 0, (unsigned char *)source_ptr );
                vec_st(score_vec1, 0, score_ptr );
                vec_st(score_vec2, 4, score_ptr );
                vec_st(score_vec3, 8, score_ptr );
                vec_st(score_vec4, 12, score_ptr );
                
                score_ptr += 16;
                source_ptr += 16;
                top_ptr += 16;
                top_left_ptr += 16;
    }
        
    score_ptr += row_size-y;
    source_ptr += row_size-y;
    top_ptr += row_size-y;
    top_left_ptr += row_size-y;
        
    for (i = 0; i < template_data->len; i++) {
                if (score_row_start[i] > score_row_start[i-1] + open_score) {
                        source_row_start[i] = vLeft;
                        score_row_start[i] = score_row_start[i-1] + open_score;
                }
                //     assert(source_row_start[i] != 0);
    }
} 
#endif





void DynamicAlign::TraceAlign(  int start_x, int start_y, int end_x, int end_y,
                                                                int *align_x, int *align_y) {   
        int x = end_x - 1;
        int y = end_y - 1;                      
        while( x >= start_x || y >= start_y) {
                if (source_matrix[x+1][y+1] == vUp) {
                        align_y[y] = -1;
                        y--;
                } else if (source_matrix[x+1][y+1] == vLeft) {
                        align_x[x] = -1;
                        x--;
                } else  if (source_matrix[x+1][y+1] == vMatch ) {
                        if (x >= 0 && x < size_x && y >= 0 && y <size_y ) {
                                align_x[x] = y;
                                align_y[y] = x;
                        }
                        x--;
                        y--;
                } else {
                        assert(0);
                        return;
                }
        }       
}


void DynamicAlign::TraceAlign( AlignmentStruct *out_align ) {
        int *target_align = (int *)malloc(sizeof(int) * size_x);
        int *template_align = (int *)malloc(sizeof(int) * size_y);
        TraceAlign( target_align, template_align );
        out_align->SetAlign( size_x, target_align, size_y, template_align);
}

void ProspectDynamicCoreAlign:: CalcEnergyTable(ProspectParam *param,
                                                                                                WeightArray weights,   
                                                                                                TemplateStruct *template_data,
                                                                                                TargetStruct *target_data ) {
        ProspectDynamicAlign::CalcEnergyTable( param, weights, template_data, target_data );
        extra_del_score = _2d_array_malloc( target_data->len, template_data->len ); //(pValType *)malloc(sizeof(pValType) * template_data->len )
        extra_ins_score = _2d_array_malloc( target_data->len, template_data->len ); //(pValType *)malloc(sizeof(pValType) * (template_data->len+1) )                    
        for (int i = 0; i < template_data->len; i++) {
                for (int j = 0; j < target_data->len; j++) {
                        extra_del_score[j][i] = calc_del_val( param, weights, template_data, target_data, i, j );
                }
        }
        for (int i = 0; i < template_data->len; i++) {
                for (int j = 0; j < target_data->len; j++) {
                        extra_ins_score[j][i] = calc_ins_val( param, weights, template_data, target_data, i, j );
                }       
        }
}


void ProspectDynamicCoreAlign::Align(int start_x, int start_y, int end_x, int end_y, char start_state, pValType start_val) {    
        int x,y;
        pValType score_left, score_up, score_match;
        
        assert(start_x >= 0 && start_x < size_x);
        assert(start_y >= 0 && start_y < size_y);
        assert(end_x >= 0 && end_x <= size_x);
        assert(end_y >= 0 && end_y <= size_y);
        assert(start_x <= end_x);
        assert(start_y <= end_y);       
        
        score_matrix[start_x][start_y] = start_val;
        source_matrix[start_x][start_y] = start_state;
        
        //fill in gap edges
        for (y = start_y+1; y < end_y+1; y++) {
                if ( start_state != vUp )
                        score_matrix[start_x][y] = start_val + (gap_open) + (gap_extend)*(y-(start_y));
                else 
                        score_matrix[start_x][y] = start_val + (gap_extend)*(y-(start_y));
                for ( int y1 = start_y+1; y1 < y; y1++) 
                        score_matrix[start_x][y] += extra_del_score[ start_x ][ y1 ];  //add in the template point deletion score
                source_matrix[start_x][y] = vUp;
        }
        for (x = start_x+1; x < end_x+1; x++) {
                if ( start_state != vLeft )
                        score_matrix[x][start_y] = start_val + (gap_open) + (gap_extend)*(x-(start_x));
                else 
                        score_matrix[x][start_y] = start_val + (gap_extend)*(x-(start_x));

                for ( int x1 = start_x; x1 < x; x1++) 
                        score_matrix[x][start_y] += extra_ins_score[ x1 ][ start_y ];  //add in the template point insertion score

                source_matrix[x][start_y] = vLeft;
        }
        
        for (x = start_x; x < end_x; x++) {
                for (y = start_y; y < end_y; y++) {
                        if ( x >= 0  && x < size_x  ) {
                                score_matrix[x+1][y+1] = energy_matrix[x][y];
                                assert(!isnan( score_matrix[x+1][y+1] ) ); 
                        } else {
                                score_matrix[x+1][y+1] = 0;
                        }
                }
        }       
        for (x = start_x; x < end_x; x++) {
                for (y = start_y; y < end_y; y++) {
                        score_match = score_matrix[x][y] + score_matrix[x+1][y+1];                      
                        score_up = score_matrix[x+1][y] + extra_del_score[ x ][ y ];
                        if (source_matrix[x+1][y] == vMatch || source_matrix[x+1][y] == vLeft ) 
                                score_up += gap_open;
                        score_up += gap_extend;
                        
                        score_left = score_matrix[x][y+1] + extra_ins_score[ x ][ y ];
                        if (source_matrix[x][y+1] == vMatch || source_matrix[x][y+1] == vUp ) 
                                score_left += gap_open;
                        score_left += gap_extend;
                        
                        if ( score_match <= score_left && score_match <= score_up) {
                                score_matrix[x+1][y+1] = score_match;
                                source_matrix[x+1][y+1] = vMatch;
                        } else {
                                if ( score_left <= score_up ) {
                                        score_matrix[x+1][y+1] = score_left;
                                        source_matrix[x+1][y+1] = vLeft;
                                } else {
                                        score_matrix[x+1][y+1] = score_up;
                                        source_matrix[x+1][y+1] = vUp;
                                }
                        }
                        assert(!isnan(score_matrix[x+1][y+1]));
                }
        }       
}               

---