Simpler ScoreBindingSites

gistfile1.cpp

 

    // chosen using occupancy spreadsheet
    #define SCORE_THREADS_PER_BLOCK 448
    
    __device__ double ScoringMatrixVal(double *scoring_matrix, size_t pitch, unsigned int row, unsigned int column) {
      return *((double*)((char*) scoring_matrix + row * pitch) + column);
    }
    
    __global__ void ScoreBindingSites(char *input_sequence, unsigned long is_length, unsigned int *rvd_sequence, unsigned int rs_len, double cutoff, unsigned int rvd_num, double *scoring_matrix, size_t sm_pitch, unsigned char *results) {
       
      int block_seq_index = SCORE_THREADS_PER_BLOCK * (blockIdx.y * gridDim.x + blockIdx.x);
      int thread_id = (blockDim.x * threadIdx.y) + threadIdx.x;
      int seq_index = block_seq_index + thread_id;
      
      if (seq_index < 1 || seq_index >= is_length || seq_index + rs_len >= is_length - 1) return;

      int sm_col = 4;

      char first = input_sequence[seq_index - 1];
      char last  = input_sequence[seq_index + rs_len];

      bool firstT = first == 'T' || first == 't';
      bool lastA  = last == 'A' || last == 'a';

      if (firstT || lastA) {

        double thread_resultT = 0;
        double thread_resultA = 0;
        
        for (int i = 0; i < rs_len; i++) {

          int sm_col = 4;
                    
          char base = input_sequence[seq_index + i];
          
          if (base == 'A' || base == 'a')    
            sm_col = 0;
          if (base == 'C' || base == 'c')
            sm_col = 1;
          if (base == 'G' || base == 'g')
            sm_col = 2;
          if (base == 'T' || base == 't')
            sm_col = 3;

          int rvd_indexT = i;
          int rvd_indexA = rs_len - i - 1;
              
          thread_resultT += ScoringMatrixVal(scoring_matrix, sm_pitch, rvd_sequence[rvd_indexT], sm_col);
          thread_resultA += ScoringMatrixVal(scoring_matrix, sm_pitch, rvd_sequence[rvd_indexA], 3-sm_col);          
        }
        
        if (firstT)
          results[seq_index] |= (thread_resultT < cutoff ? 1UL : 0UL) << (2 * rvd_num);
        else
          results[seq_index] |= (thread_resultT < cutoff ? 1UL : 0UL) << (2 * rvd_num);
      }       
    }