diff --git a/bench/nucleotide/nucleotide.c b/bench/nucleotide/nucleotide.c old mode 100644 new mode 100755 index 69edd183..273b6df1 --- a/bench/nucleotide/nucleotide.c +++ b/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 #include #include #include -#include -#include #include #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. -static inline bool oor_equal_p(uint64_t k, unsigned char n) +// Number of outputs to display +#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 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) +#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) -typedef struct { - uint64_t key; - uint32_t value; -} element; +// Define a B*TREE-21 of uint32_t --> uint64_t and register its oplist +#define TREE_CMP(a,b) ((a) > (b) ? -1 : (a) < (b)) +BPTREE_DEF2(tree_dict_oligonucleotide, 21, uint32_t, M_OPEXTEND(M_DEFAULT_OPLIST, CMP(TREE_CMP)), uint64_t, M_DEFAULT_OPLIST) +#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]) - -// 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) +static void polynucleotide_init_stream(polynucleotide_t p, FILE *stream) { - const element * const left_Element = a; - const element * const right_Element = b; - // Sort based on element values. - if(left_Element->value < right_Element->value) return 1; - if(left_Element->value > right_Element->value) return -1; - - // If we got here then both items have the same value so then sort based on - // key. - return left_Element->key > right_Element->key ? 1 : -1; + char buffer[MAXIMUM_OUTPUT_LENGTH], c; + + polynucleotide_init(p); + polynucleotide_reserve(p, MAXIMUM_OUTPUT_LENGTH); + // Find the start of the third polynucleotide. + while(fgets(buffer, sizeof(buffer), stream) && memcmp(">THREE", buffer, 6)); + // Start reading and encoding the third polynucleotide. + while(fgets(buffer, sizeof(buffer), stream) && buffer[0]!='>'){ + 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; ikey, 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-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 -generate_Count_For_Oligonucleotide(const char * const polynucleotide, - const int polynucleotide_Length, - const char * const oligonucleotide, - char * const output) +static void compute_count(polynucleotide_t p, const char oligonucleotide[], char output[]) { - const int oligonucleotide_Length=strlen(oligonucleotide); - assert (2*oligonucleotide_Length <= 64); - + const size_t oligonucleotide_Length = strlen(oligonucleotide); dict_oligonucleotide_t hash_Table; dict_oligonucleotide_init(hash_Table); - - 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; iTHREE", buffer - , sizeof(">THREE")-1)); - - // 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