From f9ddc433493af65f608b0f4b683fc9b4ccc702b4 Mon Sep 17 00:00:00 2001 From: Patrick Pelissier Date: Tue, 9 Jan 2018 18:44:11 +0100 Subject: [PATCH] Add bench k-nucleotide as per http://benchmarksgame.alioth.debian.org/u64q/knucleotide.html --- bench/nucleotide/Makefile | 22 +++ bench/nucleotide/fasta.c | 344 ++++++++++++++++++++++++++++++++++ bench/nucleotide/nucleotide.c | 292 +++++++++++++++++++++++++++++ 3 files changed, 658 insertions(+) create mode 100644 bench/nucleotide/Makefile create mode 100644 bench/nucleotide/fasta.c create mode 100644 bench/nucleotide/nucleotide.c diff --git a/bench/nucleotide/Makefile b/bench/nucleotide/Makefile new file mode 100644 index 00000000..d906758f --- /dev/null +++ b/bench/nucleotide/Makefile @@ -0,0 +1,22 @@ +CC=cc -std=c11 +CXX=c++ -std=c++11 +CFLAGS=-O2 -march=native -Wall +CPPFLAGS=-I../.. + +all: k-nucleotide + +k-nucleotide: k-nucleotide.exe data.txt + time ./k-nucleotide.exe < data.txt + +k-nucleotide.exe: nucleotide.c + $(CC) $(CFLAGS) $(CPPFLAGS) nucleotide.c -o k-nucleotide.exe -lpthread + +data.txt: fasta.exe + ./fasta.exe 25000000 > data.txt + +fasta.exe: fasta.c + $(CC) $(CFLAGS) $(CPPFLAGS) -fopenmp fasta.c -o fasta.exe + +clean: + rm -f *.exe *~ *.o data.txt + diff --git a/bench/nucleotide/fasta.c b/bench/nucleotide/fasta.c new file mode 100644 index 00000000..d27ce3ba --- /dev/null +++ b/bench/nucleotide/fasta.c @@ -0,0 +1,344 @@ +// The Computer Language Benchmarks Game +// http://benchmarksgame.alioth.debian.org/ +// +// contributed by Jeremy Zerfas +// rewritten by Аноним Легионов, inspired by fasta Rust #2 program +// use two OpenMP locks instead of one critical section +// decouples IO activity from random number generation +// +// modified by Josh Goldfoot, adding use of a buffer for fasta_repeat + + +// This controls the width of lines that are output by this program. +#define MAXIMUM_LINE_WIDTH 60 + +// This program will generate the random nucleotide sequences in parallel which +// are worked on in blocks of lines. The number of lines in those blocks is +// controlled by this setting. +#define LINES_PER_BLOCK 1024 + +#define CHARACTERS_PER_BLOCK (MAXIMUM_LINE_WIDTH*LINES_PER_BLOCK) + +#define THREADS_TO_USE 4 + +#include +#include +#include +#include + +#ifdef _OPENMP +#include +#endif + +// intptr_t should be the native integer type on most sane systems. +typedef intptr_t intnative_t; + +typedef struct +{ + char letter; + float probability; +} nucleotide_info; + + +// Repeatedly print string_To_Repeat until it has printed +// number_Of_Characters_To_Create. The output is also wrapped to +// MAXIMUM_LINE_WIDTH columns. +static char* make_Sequence_Buffer(const char string_To_Repeat[]) +{ + const intnative_t string_To_Repeat_Length = strlen(string_To_Repeat); + + // JG: Changed to that this writes a sequence to a buffer, which is used + intnative_t number_Of_Characters_To_Create = string_To_Repeat_Length * MAXIMUM_LINE_WIDTH; + char* buffer = (char*) malloc(number_Of_Characters_To_Create + + number_Of_Characters_To_Create / MAXIMUM_LINE_WIDTH + 1); + if (buffer == NULL) + exit(-1); + char* bufferOffset = buffer; + + // Create an extended_String_To_Repeat which is a copy of string_To_Repeat + // but extended with another copy of the first MAXIMUM_LINE_WIDTH characters + // of string_To_Repeat appended to the end. Later on this allows us to + // generate a line of output just by doing simple memory copies using an + // appropriate offset into extended_String_To_Repeat. + char extended_String_To_Repeat[string_To_Repeat_Length + MAXIMUM_LINE_WIDTH]; + for (intnative_t column = 0; column0;) + { + // Figure out the length of the line we need to write. If it's less than + // MAXIMUM_LINE_WIDTH then we also need to add a line feed in the right + // spot too. + intnative_t line_Length = MAXIMUM_LINE_WIDTH; + if (current_Number_Of_Characters_To_Createstring_To_Repeat_Length) + offset -= string_To_Repeat_Length; + + // "Output" that line to our buffer and update the + // current_Number_Of_Characters_To_Create. + memcpy(bufferOffset, line, line_Length + 1); // JG: used to be fwrite(line, line_Length + 1, 1, stdout); + bufferOffset += line_Length + 1; + current_Number_Of_Characters_To_Create -= line_Length; + } + *bufferOffset = 0; + return buffer; +} + +void repeat_And_Wrap_String(const char string_To_Repeat[], intnative_t number_Of_Characters_To_Create) +{ + /* JG: fasta_repeat repeats every len(alu) * line-length = 287 * 61 = 17507 characters. + So, calculate this once, then just print that buffer over and over. */ + char* sequence = make_Sequence_Buffer(string_To_Repeat); + intnative_t sequenceLen = (intnative_t) strlen(sequence); + intnative_t outputBytes = number_Of_Characters_To_Create + number_Of_Characters_To_Create / 60; + while (outputBytes >= sequenceLen) { + fwrite(sequence, sequenceLen, 1, stdout); + outputBytes -= sequenceLen; + } + if (outputBytes > 0) { + fwrite(sequence, outputBytes, 1, stdout); + printf("\n"); + } + free(sequence); +} + + +// Generate a pseudorandom number from 0 to max using a linear +// congruential generator. +#define IM 139968 +#define IA 3877 +#define IC 29573 +uint32_t seed = 42; + +int rng_tid; //Thread ID +int rng_tnum = 1; //Thread number +intnative_t rng_cnt = 0; + +#ifdef _OPENMP +omp_lock_t rng_lock; +#define RNG_LOCK_INIT() omp_init_lock(&rng_lock) +#define RNG_LOCK() omp_set_lock(&rng_lock) +#define RNG_FREE() omp_unset_lock(&rng_lock) +#else +#define RNG_LOCK_INIT() do{}while(0) +#define RNG_LOCK() do{}while(0) +#define RNG_FREE() do{}while(0) +#endif + +static void rng_init(void) +{ + RNG_LOCK_INIT(); + rng_tid = 0; +} + +static intnative_t rng_gen_blk(uint32_t * buf, intnative_t len, int curr_tid) +{ + intnative_t gen_cnt = -1;//Error by default + RNG_LOCK(); + + if (rng_tid == curr_tid) + { + if (++rng_tid >= rng_tnum) + { + rng_tid = 0; + } + + gen_cnt = (len= out_tnum) + { + out_tid = 0; + } + wr_cnt = fwrite(buf, len, 1, stdout); + } + + OUT_FREE(); + return wr_cnt; //-1 - thread error, 0 - IO error, 1 - ОК +} + +static void generate_And_Wrap_Pseudorandom_DNA_Sequence( + const nucleotide_info nucl_info[], + const intnative_t nucl_num, + const intnative_t char_num) +{ + uint32_t cumul_p[nucl_num]; + float cumul_acc = 0.0; + + for (intnative_t i = 0; iTHREADS_TO_USE) tnum = THREADS_TO_USE; + omp_set_num_threads(tnum); + rng_tnum = tnum; + out_tnum = tnum; +#endif + rng_tid = 0; + out_tid = 0; + rng_cnt = char_num; + +#pragma omp parallel + { + char block[CHARACTERS_PER_BLOCK + LINES_PER_BLOCK]; + char * line; + uint32_t rnd[CHARACTERS_PER_BLOCK], r; + intnative_t cnt, col, prid, nid, ncnt; + int cur_tid; +#ifdef _OPENMP + cur_tid = omp_get_thread_num(); +#else + cur_tid = 0; +#endif + while (1) + { + do + { + cnt = rng_gen_blk(rnd, CHARACTERS_PER_BLOCK, cur_tid); + } while (-1 == cnt); + + if (0 == cnt) + { + break;//Work finished! + } + + line = block; + + for (col = 0, prid = 0; prid < cnt; prid++) + { + r = rnd[prid]; + + ncnt = 0; + for (nid = 0; nid < nucl_num; nid++) + { + if (cumul_p[nid] <= r) + { + ncnt++; + } + } + + *line++ = nucl_info[ncnt].letter; + if (++col >= MAXIMUM_LINE_WIDTH) + { + col = 0; + *line++ = '\n'; + } + } + //Check if we need to end the line + if (0 != col) + { + //Last iteration didn't end the line, so finish the job. + *line++ = '\n'; + } + //Print results + do + { + cnt = out_write(block, line - block, cur_tid); + } while (-1 == cnt); + //Check fot IO error + if (0 == cnt) + { + exit(1); + } + } + } +} + + +int main(int argc, char ** argv) +{ + const intnative_t n = atoi(argv[1]); + + fputs(">ONE Homo sapiens alu\n", stdout); + const char homo_Sapiens_Alu[] = + "GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGGGAGGCCGAGGCGGGCGGATCACCTGAGGTC" + "AGGAGTTCGAGACCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAATACAAAAATTAGCCGGGCG" + "TGGTGGCGCGCGCCTGTAATCCCAGCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGGAGGCGG" + "AGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCCAGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA"; + repeat_And_Wrap_String(homo_Sapiens_Alu, 2 * n); + + rng_init(); + out_init(); + + fputs(">TWO IUB ambiguity codes\n", stdout); + nucleotide_info iub_Nucleotides_Information[] = + { + { 'a', 0.27 },{ 'c', 0.12 },{ 'g', 0.12 },{ 't', 0.27 },{ 'B', 0.02 }, + { 'D', 0.02 },{ 'H', 0.02 },{ 'K', 0.02 },{ 'M', 0.02 },{ 'N', 0.02 }, + { 'R', 0.02 },{ 'S', 0.02 },{ 'V', 0.02 },{ 'W', 0.02 },{ 'Y', 0.02 } + }; + generate_And_Wrap_Pseudorandom_DNA_Sequence(iub_Nucleotides_Information, + sizeof(iub_Nucleotides_Information) / sizeof(nucleotide_info), 3 * n); + + fputs(">THREE Homo sapiens frequency\n", stdout); + nucleotide_info homo_Sapien_Nucleotides_Information[] = + { + { 'a', 0.3029549426680 },{ 'c', 0.1979883004921 }, + { 'g', 0.1975473066391 },{ 't', 0.3015094502008 } + }; + generate_And_Wrap_Pseudorandom_DNA_Sequence(homo_Sapien_Nucleotides_Information, + sizeof(homo_Sapien_Nucleotides_Information) / sizeof(nucleotide_info), 5 * n); + + return 0; +} diff --git a/bench/nucleotide/nucleotide.c b/bench/nucleotide/nucleotide.c new file mode 100644 index 00000000..515def60 --- /dev/null +++ b/bench/nucleotide/nucleotide.c @@ -0,0 +1,292 @@ +// The Computer Language Benchmarks Game +// http://benchmarksgame.alioth.debian.org/ +// +// 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 + +// This controls the maximum length for each set of oligonucleotide frequencies +// and each oligonucleotide count output by this program. +#define NDEBUG + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "m-worker.h" +#include "m-dict.h" + +#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) +{ + return k == ~(uint64_t) n; +} +static inline void 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 dictionnary +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)), + 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) + +typedef struct { + uint64_t key; + uint32_t value; +} element; + +// 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 +// need the three least significant bits to differentiate the letters 'A', 'C', +// 'G', and 'T'. Spaces in this array/string will never be used as long as +// 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) +{ + 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; +} + + +// 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}); + } + 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 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) +{ + const int oligonucleotide_Length=strlen(oligonucleotide); + assert (2*oligonucleotide_Length <= 64); + + 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