Full refactorisation and proper use of M*LIB.

Should be as fast as before.

bench/nucleotide/nucleotide.c

@@ -3,50 +3,50 @@
 //
 // Original version Contributed by Jeremy Zerfas
 // M*LIB port contributed by Patrick Pélissier
-// gcc -O2 -march=native nucleotide.c -I../.. -std=c11 -lpthread -DNDEBUG
-
 #define NDEBUG
-
 #include <stdint.h>
 #include <stdlib.h>
 #include <stdio.h>
 #include <string.h>
-#include <time.h>
-#include <sys/time.h>
 #include <assert.h>
 
 #include "m-worker.h"
 #include "m-dict.h"
+#include "m-bptree.h"
 
 // This controls the maximum length for each set of oligonucleotide frequencies
 // and each oligonucleotide count output by this program.
 #define MAXIMUM_OUTPUT_LENGTH 4096
 
-// Define Open Addressing dictionnary. Need to define how to handle out of bound data.
+// Number of outputs to display
-static inline bool oor_equal_p(uint64_t k, unsigned char n)
+#define MAXMUM_NUMBER_OUTPUT 7
+
+/* Define the needed containers */
+/* Define Open Addressing dictionary of uint64_t --> uint32_t. Need to define how to handle out of bound data. */
+static inline bool my_oor_equal_p(uint64_t k, unsigned char n)
 {
   return k == ~(uint64_t) n;
 }
-static inline void oor_set(uint64_t *k, unsigned char n)
+static inline void my_oor_set(uint64_t *k, unsigned char n)
 {
   *k = ~(uint64_t)n;
 }
-// Custom hash function.
 static inline size_t my_hash(uint64_t key)
 {
   return ((key) ^ (key)>>7);
 }
-// Define a dictionnary of uint64_t --> uint32_t
 DICT_OA_DEF2(dict_oligonucleotide,
-             uint64_t, M_OPEXTEND(M_DEFAULT_OPLIST, HASH(my_hash), OOR_EQUAL(oor_equal_p), OOR_SET(oor_set M_IPTR)),
+             uint64_t, M_OPEXTEND(M_DEFAULT_OPLIST, HASH(my_hash), OOR_EQUAL(my_oor_equal_p), OOR_SET(my_oor_set M_IPTR)),
              uint32_t, M_DEFAULT_OPLIST)
-// Define oplist of dictionary
+#define M_OPL_dict_oligonucleotide_t() DICT_OPLIST(dict_oligonucleotide, M_OPEXTEND(M_DEFAULT_OPLIST, HASH(my_hash), OOR_EQUAL(my_oor_equal_p), OOR_SET(my_oor_set M_IPTR)), M_DEFAULT_OPLIST)
-#define OLIGONUCLEOTIDE_OPLIST DICT_OPLIST(dict_oligonucleotide, M_OPEXTEND(M_DEFAULT_OPLIST, HASH(oor_hash), OOR_EQUAL(oor_equal_p), OOR_SET(oor_set M_IPTR)), M_DEFAULT_OPLIST)
 
-typedef struct {
+// Define a B*TREE-21 of uint32_t --> uint64_t and register its oplist
-  uint64_t	key;
+#define TREE_CMP(a,b)      ((a) > (b) ? -1 : (a) < (b))
-  uint32_t	value;
+BPTREE_DEF2(tree_dict_oligonucleotide, 21, uint32_t, M_OPEXTEND(M_DEFAULT_OPLIST, CMP(TREE_CMP)), uint64_t, M_DEFAULT_OPLIST)
-} element;
+#define M_OPL_tree_dict_oligonucleotide_t() BPTREE_OPLIST2(tree_dict_oligonucleotide, M_DEFAULT_OPLIST, M_DEFAULT_OPLIST)
+
+/* Define the data to analyse (a plain array) */
+ARRAY_DEF(polynucleotide, unsigned char)
 
 // Macro to convert a nucleotide character to a code. Note that upper and lower
 // case ASCII letters only differ in the fifth bit from the right and we only
@@ -55,143 +55,105 @@ typedef struct {
 // characters other than 'A', 'C', 'G', and 'T' aren't used.
 #define CODE_FOR_NUCLEOTIDE(nucleotide) (" \0 \1\3  \2"[nucleotide & 0x7])
 
-
 // And one more macro to convert the codes back to nucleotide characters.
 #define NUCLEOTIDE_FOR_CODE(code) ("ACGT"[code & 0x3])
 
-
+static void polynucleotide_init_stream(polynucleotide_t p, FILE *stream)
-// Function to use when sorting elements with qsort() later. Elements with
-// larger values will come first and in cases of identical values then elements
-// with smaller keys will come first.
-static int element_Compare(const void *a, const void *b)
 {
-  const element * const left_Element = a;
+  char buffer[MAXIMUM_OUTPUT_LENGTH], c;
-  const element * const right_Element = b;
+
-  // Sort based on element values.
+  polynucleotide_init(p);
-  if(left_Element->value < right_Element->value) return 1;
+  polynucleotide_reserve(p, MAXIMUM_OUTPUT_LENGTH);
-  if(left_Element->value > right_Element->value) return -1;
+  // Find the start of the third polynucleotide.
-  
+  while(fgets(buffer, sizeof(buffer), stream) && memcmp(">THREE", buffer, 6));
-  // If we got here then both items have the same value so then sort based on
+  // Start reading and encoding the third polynucleotide.
-  // key.
+  while(fgets(buffer, sizeof(buffer), stream) && buffer[0]!='>'){
-  return left_Element->key > right_Element->key ? 1 : -1;
+    for(size_t i = 0; (c = buffer[i]) !='\0'; i++) {
+      if (c != '\n') {
+          polynucleotide_push_back(p, CODE_FOR_NUCLEOTIDE(c));
+      }
+    }
+  }
 }
 
+static void compute_hash(dict_oligonucleotide_t hash_Table, const polynucleotide_t p, const size_t desiredLength)
+{
+  uint64_t key = 0;
+  const uint64_t mask = (1ULL<<(2*desiredLength))-1;
+  assert (2*desiredLength < 64 && desiredLength < polynucleotide_size(p));
 
-// Generate frequencies for all oligonucleotides in polynucleotide that are of
-// desired_Length_For_Oligonucleotides and then save it to output.
-static void
-generate_Frequencies_For_Desired_Length_Oligonucleotides(const char polynucleotide[],
-                                                         const int polynucleotide_Length,
-                                                         const int desired_Length_For_Oligonucleotides,
-                                                         char * const output)
-{  
-  assert (2*desired_Length_For_Oligonucleotides <= 64);
-
-  dict_oligonucleotide_t hash_Table;
-  dict_oligonucleotide_init(hash_Table);
-  
-  uint64_t key=0;
-  const uint64_t mask=(1ULL<<(2*desired_Length_For_Oligonucleotides))-1;
-  
   // For the first several nucleotides we only need to append them to key in
   // preparation for the insertion of complete oligonucleotides to hash_Table.
-  for(int i=0; i<desired_Length_For_Oligonucleotides-1; i++)
+  for(size_t i = 0; i < desiredLength - 1; i++)
-    key=(key<<2 & mask) | polynucleotide[i];
+    key = (key<<2) | *polynucleotide_get(p, i);
 
   // Add all the complete oligonucleotides of
-  // desired_Length_For_Oligonucleotides to hash_Table and update the count
+  // desiredLength to hash_Table, updating their count
-  // for each oligonucleotide.
+  for(size_t i=desiredLength-1 ; i < polynucleotide_size(p); i++) {
-  for(int i=desired_Length_For_Oligonucleotides-1
+    key = (key<<2 & mask) | *polynucleotide_get(p, i);
-        ; i<polynucleotide_Length; i++){
-    
-    key=(key<<2 & mask) | polynucleotide[i];
-    
     uint32_t *p = dict_oligonucleotide_get(hash_Table, key);
-    if (p != NULL)
+    if (p != NULL) {
       (*p) ++;
-    else {
+    } else {
       dict_oligonucleotide_set_at(hash_Table, key, 1);
     }
   }
-  
-  // Create an array of elements from hash_Table.
-  int elements_Array_Size=dict_oligonucleotide_size(hash_Table), i=0;
-  element * elements_Array=malloc(elements_Array_Size*sizeof(element));
-  if (!elements_Array) abort();
-  for M_EACH(item, hash_Table, OLIGONUCLEOTIDE_OPLIST) {
-      elements_Array[i++]=((element){item->key, item->value});
-    }
-  // Clear the hash table
-  dict_oligonucleotide_clear(hash_Table);
-  
-  // Sort elements_Array.
-  qsort(elements_Array, elements_Array_Size, sizeof(element), element_Compare);
-  
-  // Print the frequencies for each oligonucleotide.
-  for(int output_Position=0, i=0; i<elements_Array_Size; i++) {
-    
-    // Convert the key for the oligonucleotide to a string.
-    char oligonucleotide[desired_Length_For_Oligonucleotides+1];
-    for(int j=desired_Length_For_Oligonucleotides-1; j>-1; j--){
-      oligonucleotide[j]=NUCLEOTIDE_FOR_CODE(elements_Array[i].key);
-      elements_Array[i].key>>=2;
-    }
-    oligonucleotide[desired_Length_For_Oligonucleotides]='\0';
-    
-    // Output the frequency for oligonucleotide to output.
-    output_Position+=snprintf(output+output_Position
-                              , MAXIMUM_OUTPUT_LENGTH-output_Position, "%s %.3f\n", oligonucleotide
-                              , 100.0f*elements_Array[i].value
-                              /(polynucleotide_Length-desired_Length_For_Oligonucleotides+1));
-  }
-  
-  free(elements_Array);
 }
 
+// Generate frequencies for all oligonucleotides in polynucleotide that are of
+// desiredLength and then save it to output.
+static void compute_freq(const polynucleotide_t p, const size_t desiredLength, char output[])
+{  
+  dict_oligonucleotide_t hash_Table;
+  dict_oligonucleotide_init(hash_Table);
+  compute_hash(hash_Table, p, desiredLength);
+
+  // Order key is the count frequency (inverse value & key)
+  tree_dict_oligonucleotide_t tree;
+  tree_dict_oligonucleotide_init(tree);
+  for M_EACH(item, hash_Table, dict_oligonucleotide_t) {
+    tree_dict_oligonucleotide_set_at (tree, item->value, item->key);
+  }
+
+  // Create Output
+  char oligonucleotide[64];
+  const double invSum = 100.0 / (double) (polynucleotide_size(p) - desiredLength + 1);
+  int output_Position = 0;
+  for M_EACH(item, tree, tree_dict_oligonucleotide_t) {
+    // Convert the key for the oligonucleotide to a string.
+    uint64_t freq  = *item->value_ptr;
+    uint32_t count = *item->key_ptr;
+    for(size_t j = desiredLength-1; j < desiredLength; j--){
+      oligonucleotide[j]=NUCLEOTIDE_FOR_CODE(freq);
+      freq>>=2;
+    }
+    oligonucleotide[desiredLength]='\0';
+    // Output the frequency for oligonucleotide to output.
+    output_Position+=snprintf(output+output_Position, MAXIMUM_OUTPUT_LENGTH-output_Position,
+                              "%s %.3f\n", oligonucleotide, invSum * count);
+  }
+  
+  dict_oligonucleotide_clear(hash_Table);
+  tree_dict_oligonucleotide_init(tree);
+}
 
 // Generate a count for the number of times oligonucleotide appears in
 // polynucleotide and then save it to output.
-static void
+static void compute_count(polynucleotide_t p, const char oligonucleotide[], char output[])
-generate_Count_For_Oligonucleotide(const char * const polynucleotide,
-                                   const int polynucleotide_Length,
-                                   const char * const oligonucleotide,
-                                   char * const output)
 {
-  const int oligonucleotide_Length=strlen(oligonucleotide);
+  const size_t oligonucleotide_Length = strlen(oligonucleotide);
-  assert (2*oligonucleotide_Length <= 64);
-  
   dict_oligonucleotide_t hash_Table;
   dict_oligonucleotide_init(hash_Table);
-  
+  compute_hash(hash_Table, p, oligonucleotide_Length);
-  uint64_t key = 0;
-  const uint64_t mask=(1ULL<<(2*oligonucleotide_Length))-1;
-  
-  // For the first several nucleotides we only need to append them to key in
-  // preparation for the insertion of complete oligonucleotides to hash_Table.
-  for(int i=0; i<oligonucleotide_Length-1; i++)
-    key=(key<<2 & mask) | polynucleotide[i];
-
-  // Add all the complete oligonucleotides of oligonucleotide_Length to
-  // hash_Table and update the count for each oligonucleotide.
-  for(int i=oligonucleotide_Length-1; i<polynucleotide_Length; i++){
-    key=(key<<2 & mask) | polynucleotide[i];
-    
-    uint32_t *p = dict_oligonucleotide_get(hash_Table, key);
-    if (p != NULL)
-      (*p) ++;
-    else {
-      dict_oligonucleotide_set_at(hash_Table, key, 1);
-    }
-  }
   
   // Generate the key for oligonucleotide.
-  key=0;
+  uint64_t key=0;
-  for(int i=0; i<oligonucleotide_Length; i++)
+  for(size_t i=0; i<oligonucleotide_Length; i++)
-    key=(key<<2) | CODE_FOR_NUCLEOTIDE(oligonucleotide[i]);
+    key = (key<<2) | CODE_FOR_NUCLEOTIDE(oligonucleotide[i]);
   
   // Output the count for oligonucleotide to output.
-  uint32_t *p = dict_oligonucleotide_get(hash_Table, key);
+  uint32_t *val = dict_oligonucleotide_get(hash_Table, key);
-  uintmax_t count= p== NULL ? 0 : *p;
+  uintmax_t count= val== NULL ? 0 : *val;
   snprintf(output, MAXIMUM_OUTPUT_LENGTH, "%ju\t%s", count, oligonucleotide);
   
   dict_oligonucleotide_clear(hash_Table);
@@ -199,94 +161,49 @@ generate_Count_For_Oligonucleotide(const char * const polynucleotide,
 
 int main()
 {
-  char buffer[MAXIMUM_OUTPUT_LENGTH];
-  worker_t worker;
-
   // Initialize the different worker threads by auto-detecting the number of cores of the system
+  worker_t worker;
   worker_init(worker, 0, 1, NULL);
+
+  // Read the input
+  polynucleotide_t p, *ptr = &p;
+  polynucleotide_init_stream(p, stdin);
   
-  // Find the start of the third polynucleotide.
+  // Prepare outputs
-  while(fgets(buffer, sizeof(buffer), stdin) && memcmp(">THREE", buffer
+  char outputBuffer[MAXMUM_NUMBER_OUTPUT][MAXIMUM_OUTPUT_LENGTH];
-                                                       , sizeof(">THREE")-1));
+  char (*output)[MAXIMUM_OUTPUT_LENGTH] = outputBuffer;
-  
-  // Start with 1 MB of storage for reading in the polynucleotide and grow
-  // geometrically.
-  int polynucleotide_Capacity=1048576;
-  int polynucleotide_Length=0;
-  char * polynucleotide=malloc(polynucleotide_Capacity);
-  if (!polynucleotide) abort();
-  
-  // Start reading and encoding the third polynucleotide.
-  while(fgets(buffer, sizeof(buffer), stdin) && buffer[0]!='>'){
-    for(int i=0; buffer[i]!='\0'; i++)
-      if(buffer[i]!='\n')
-        polynucleotide[polynucleotide_Length++]
-          =CODE_FOR_NUCLEOTIDE(buffer[i]);
-    
-    // Make sure we still have enough memory allocated for any potential
-    // nucleotides in the next line.
-    if(polynucleotide_Capacity-polynucleotide_Length<sizeof(buffer)) {
-      polynucleotide=realloc(polynucleotide, polynucleotide_Capacity*=2);
-      if (!polynucleotide) abort();
-    }
-  }
-  
-  char output_Buffer[7][MAXIMUM_OUTPUT_LENGTH];
-  char (*output)[MAXIMUM_OUTPUT_LENGTH] = output_Buffer;
   
   // Do the following functions in parallel.
-  // Define a synchronization block
+  worker_sync_t block;   // Define a synchronization block
-  worker_sync_t block;
   worker_start(block);
   
-  WORKER_SPAWN(worker, block, (polynucleotide, polynucleotide_Length, output), {
+  WORKER_SPAWN(worker, block, (ptr, output), {
-      generate_Count_For_Oligonucleotide(polynucleotide,
+      compute_count(*ptr, "GGTATTTTAATTTATAGT", output[6]);
-                                         polynucleotide_Length,
-                                         "GGTATTTTAATTTATAGT", output[6]);
     }, (/*no output*/));
+  WORKER_SPAWN(worker, block, (ptr, output), {
+      compute_count(*ptr, "GGTATTTTAATT", output[5]);
+    }, (/*no output*/));
+  WORKER_SPAWN(worker, block, (ptr, output), {
+      compute_count(*ptr, "GGTATT", output[4]);
+    }, (/*no output*/));
+  WORKER_SPAWN(worker, block, (ptr, output), {
+      compute_count(*ptr, "GGTA", output[3]);
+    }, (/*no output*/));
+  WORKER_SPAWN(worker, block, (ptr, output), {
+      compute_count(*ptr, "GGT", output[2]);  
+    }, (/*no output*/));
+  WORKER_SPAWN(worker, block, (ptr, output), {
+      compute_freq(*ptr, 2, output[1]);
+    }, (/*no output*/));
+  compute_freq(p, 1, output[0]);
 
-  WORKER_SPAWN(worker, block, (polynucleotide, polynucleotide_Length, output), {
-      generate_Count_For_Oligonucleotide(polynucleotide,
-                                         polynucleotide_Length,
-                                         "GGTATTTTAATT", output[5]);
-    }, (/*no output*/));
-  
-  WORKER_SPAWN(worker, block, (polynucleotide, polynucleotide_Length, output), {
-      generate_Count_For_Oligonucleotide(polynucleotide,
-                                         polynucleotide_Length,
-                                         "GGTATT", output[4]);
-    }, (/*no output*/));
-  
-  WORKER_SPAWN(worker, block, (polynucleotide, polynucleotide_Length, output), {
-      generate_Count_For_Oligonucleotide(polynucleotide,
-                                         polynucleotide_Length,
-                                         "GGTA", output[3]);
-    }, (/*no output*/));
-  
-  WORKER_SPAWN(worker, block, (polynucleotide, polynucleotide_Length, output), {
-      generate_Count_For_Oligonucleotide(polynucleotide,
-                                         polynucleotide_Length,
-                                         "GGT", output[2]);  
-    }, (/*no output*/));
-
-  WORKER_SPAWN(worker, block, (polynucleotide, polynucleotide_Length, output), {
-      generate_Frequencies_For_Desired_Length_Oligonucleotides(polynucleotide,
-                                                               polynucleotide_Length,
-                                                               2,
-                                                               output[1]);
-    }, (/*no output*/));
-
-  generate_Frequencies_For_Desired_Length_Oligonucleotides(polynucleotide,
-                                                           polynucleotide_Length,
-                                                           1,
-                                                           output[0]);
   worker_sync(block);
   
   // Output the results to stdout.
-  for(int i=0; i<7; i++)
+  for(int i=0; i<MAXMUM_NUMBER_OUTPUT; i++)
-    printf("%s\n", output_Buffer[i]);
+    printf("%s\n", outputBuffer[i]);
   
-  free(polynucleotide);
+  polynucleotide_clear(p);
   worker_clear(worker);
   
   return 0;