Merge pull request #3 from criniguez/dev

Adding vectorized banded bpm implementation.

Author: maxdoblas

.github/workflows/build_and_test.yaml

@@ -9,6 +9,9 @@ on:
       - '**/LICENSE'
   pull_request:
     branches: [ main ]
+    paths-ignore: # Skip the workflow if only documentation files are changed
+      - '**/*.md'
+      - '**/LICENSE'
   workflow_dispatch:
 
 jobs:

CMakeLists.txt

@@ -24,6 +24,7 @@ if (CMAKE_BUILD_TYPE STREQUAL "Debug")
 
     # Ignore known non-concerning warnings
     add_compile_options("$<$<COMPILE_LANGUAGE:C>:-Wno-override-init>")
+    add_compile_options(-Wno-error=cpp) # #warning are not errors
 
     # Enable code coverage
     add_compile_options(--coverage)

README.md

@@ -230,7 +230,7 @@ The `quicked_params_t` configuration struct has the following parameters:
 * **unsigned int** `hew_threshold[2]`: The error percentage threshold inside a window to be considered a high error window (HEW). This parameter is only used inside Quicked. Position [0] refers to the WindowEd(S) step and position [1] to the WindowEd(L) step.
 * **unsigned int** `hew_percentage[2]`: percentage of HEW in a particular WindowEd alignment to consider that the estimation is not fitted. This parameter is only used inside Quicked. Position [0] refers to the WindowEd(S) step and position [1] to the WindowEd(L) step.
 * **bool** `only_score`: If set to true, turn off the CIGAR generation for the WindowEd and BandEd methods.
-* **bool** `force_scalar`: If set to true, it forces the WindowEd implementation to use the scalar code.
+* **bool** `force_scalar`: If set to true, it forces WindowEd and BandEd implementation to use the scalar code.
 
 > [!WARNING]
 > **Experimental Parameters**

examples/bindings/basic.cpp

@@ -42,8 +42,8 @@ int main(void) {
 
         aligner.align(&pattern, &text);  // Align the sequences!
 
-        aligner.getScore();              // Get the score
-        aligner.getCigar();              // Get the CIGAR string
+        score = aligner.getScore();              // Get the score
+        cigar = aligner.getCigar();              // Get the CIGAR string
     } catch (quicked::QuickedException &e) {
         cerr << e.what() << endl;
         return 1;
@@ -53,4 +53,4 @@ int main(void) {
     cout << "Cigar: " << cigar << endl;  // Print the CIGAR string
 
     return 0;
-}
+}

examples/bindings/basic.py

@@ -40,6 +40,6 @@
 except QuickedEception as e:
     print(e)
 
-print(f"Score: {aligner.getScore()}") # Print the score
-print(f"Cigar: {aligner.getCigar()}")  # Print the CIGAR string
+print(f"Score: {score}") # Print the score
+print(f"Cigar: {cigar}")  # Print the CIGAR string
 

examples/bindings/params.cpp

@@ -47,8 +47,8 @@ int main(void) {
 
         aligner.align(&pattern, &text);         // Align the sequences!
 
-        aligner.getScore();                     // Get the score
-        aligner.getCigar();                     // Get the CIGAR string
+        score = aligner.getScore();                     // Get the score
+        cigar = aligner.getCigar();                     // Get the CIGAR string
     } catch (quicked::QuickedException &e) {
         cerr << e.what() << endl;
         return 1;
@@ -58,4 +58,4 @@ int main(void) {
     cout << "Cigar <Expecting NULL>: " << cigar << endl;   // We didn't compute the CIGAR, so it's NULL
 
     return 0;
-}
+}

quicked/include/bpm_banded.h

@@ -90,6 +90,7 @@ void banded_compute(
     const char* text,
     const int64_t text_length,
     const int64_t text_finish_pos,
-    const bool only_score);
+    const bool only_score,
+    const bool force_scalar);
 
 #endif /* BPM_BANDED_H_ */

quicked/include/bpm_commons.h

@@ -100,4 +100,53 @@
     Pv = Mh | ~(Xv | Ph);                                               \
     Mv = Ph & Xv
 
+#define BPM_ADVANCE_BLOCK_SI256(Eq, mask, Pv, Mv, PHin, MHin, PHout, MHout) 	       \
+    __m256i Xv    = _mm256_or_si256(Eq, Mv);      /*Eq | Mv*/                          \
+    __m256i _Eq   = _mm256_or_si256(Eq, MHin);    /*Eq | MHin*/      	       	       \
+    __m256i Xh    = _mm256_and_si256(_Eq, Pv);    /*(((_Eq & Pv) + Pv) ^ Pv) | _Eq*/   \
+    	    Xh    = _mm256_add_epi64(Xh, Pv);					                       \
+    	    Xh    = _mm256_xor_si256(Xh, Pv); 					                       \
+    	    Xh    = _mm256_or_si256(Xh, _Eq);					                       \
+    __m256i Ph    = _mm256_or_si256(Xh, Pv);      /*Mv | ~(Xh | Pv)*/                  \
+    	    Ph    = _mm256_or_si256(Mv, ~Ph);					                       \
+    __m256i Mh 	  = _mm256_and_si256(Pv, Xh);     /*Pv & Xh*/			               \
+            PHout = _mm256_and_si256(Ph, mask);   /*(Ph & mask) != 0*/	        	   \
+    	    PHout = _mm256_cmpeq_epi64(PHout, _mm256_setzero_si256());		           \
+            PHout = _mm256_andnot_si256(PHout, _mm256_set1_epi64x(1));				   \
+            MHout = _mm256_and_si256(Mh, mask);   /*(Mh & mask) != 0 */	               \
+    	    MHout = _mm256_cmpeq_epi64(MHout, _mm256_setzero_si256());				   \
+            MHout = _mm256_andnot_si256(MHout, _mm256_set1_epi64x(1));				   \
+      	    Ph 	  = _mm256_slli_epi64(Ph, 1);     /*Ph <<= 1*/		            	   \
+      	    Mh 	  = _mm256_slli_epi64(Mh, 1);     /*Mh <<= 1*/		            	   \
+      	    Ph    = _mm256_or_si256(Ph, PHin);    /*Ph |= PHin*/	        		   \
+      	    Mh    = _mm256_or_si256(Mh, MHin);    /*Mh |= MHin*/		        	   \
+      	    Pv    = _mm256_or_si256(Xv, Ph);      /*Mh | ~(Xv | Ph)*/      		       \
+      	    Pv    = _mm256_or_si256(Mh, ~Pv);                                          \
+            Mv    = _mm256_and_si256(Ph, Xv);                                          \
+
+
+#define BPM_ADVANCE_BLOCK_SI256_2(Eq2, mask2, Pv2, Mv2, PHin2, MHin2, PHout2, MHout2) 	   \
+    __m256i Xv2    = _mm256_or_si256(Eq2, Mv2);      /*Eq | Mv*/                           \
+    __m256i _Eq2   = _mm256_or_si256(Eq2, MHin2);    /*Eq | MHin*/      	       	       \
+    __m256i Xh2    = _mm256_and_si256(_Eq2, Pv2);    /*(((_Eq & Pv) + Pv) ^ Pv) | _Eq*/    \
+    	    Xh2    = _mm256_add_epi64(Xh2, Pv2);					                       \
+    	    Xh2    = _mm256_xor_si256(Xh2, Pv2); 					                       \
+    	    Xh2    = _mm256_or_si256(Xh2, _Eq2);					                       \
+    __m256i Ph2    = _mm256_or_si256(Xh2, Pv2);      /*Mv | ~(Xh | Pv)*/                   \
+    	    Ph2    = _mm256_or_si256(Mv2, ~Ph2);					                       \
+    __m256i Mh2 	  = _mm256_and_si256(Pv2, Xh2);     /*Pv & Xh*/			               \
+            PHout2 = _mm256_and_si256(Ph2, mask2);   /*(Ph & mask) != 0*/	        	   \
+    	    PHout2 = _mm256_cmpeq_epi64(PHout2, _mm256_setzero_si256());		           \
+            PHout2 = _mm256_andnot_si256(PHout2, _mm256_set1_epi64x(1));				   \
+            MHout2 = _mm256_and_si256(Mh2, mask2);   /*(Mh & mask) != 0 */	               \
+    	    MHout2 = _mm256_cmpeq_epi64(MHout2, _mm256_setzero_si256());				   \
+            MHout2 = _mm256_andnot_si256(MHout2, _mm256_set1_epi64x(1));				   \
+      	    Ph2 	  = _mm256_slli_epi64(Ph2, 1);     /*Ph <<= 1*/		            	   \
+      	    Mh2 	  = _mm256_slli_epi64(Mh2, 1);     /*Mh <<= 1*/		            	   \
+      	    Ph2    = _mm256_or_si256(Ph2, PHin2);    /*Ph |= PHin*/	        		       \
+      	    Mh2    = _mm256_or_si256(Mh2, MHin2);    /*Mh |= MHin*/		        	       \
+      	    Pv2    = _mm256_or_si256(Xv2, Ph2);      /*Mh | ~(Xv | Ph)*/      		       \
+      	    Pv2    = _mm256_or_si256(Mh2, ~Pv2);                                           \
+            Mv2    = _mm256_and_si256(Ph2, Xv2);                                           \
+
 #endif /* BPM_COMMON_H_ */

quicked/include/bpm_hirschberg.h

@@ -37,6 +37,7 @@ quicked_status_t bpm_compute_matrix_hirschberg(
     const int64_t pattern_length,
     const int64_t cutoff_score,
     cigar_t* cigar_out,
+    const bool force_scalar,
     mm_allocator_t* const mm_allocator);
 
 #endif /* BPM_HIRSCHBERG_H_ */