Add bench k-nucleotide as per

http://benchmarksgame.alioth.debian.org/u64q/knucleotide.html
This commit is contained in:
Patrick Pelissier committed 2018-01-09 18:44:11 +01:00
1 parent 1c9a232a9d
commit f9ddc43349
3 files changed
+658

No files matched your search

+22
View File
@@ -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
+344
View File
@@ -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 <stdint.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#ifdef _OPENMP
#include <omp.h>
#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; column<string_To_Repeat_Length + MAXIMUM_LINE_WIDTH;
column++)
extended_String_To_Repeat[column] =
string_To_Repeat[column%string_To_Repeat_Length];
intnative_t offset = 0;
char line[MAXIMUM_LINE_WIDTH + 1];
line[MAXIMUM_LINE_WIDTH] = '\n';
for (intnative_t current_Number_Of_Characters_To_Create =
number_Of_Characters_To_Create;
current_Number_Of_Characters_To_Create>0;)
{
// 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_Create<MAXIMUM_LINE_WIDTH)
{
line_Length = current_Number_Of_Characters_To_Create;
line[line_Length] = '\n';
}
memcpy(line, extended_String_To_Repeat + offset, line_Length);
// Update the offset, reducing it by string_To_Repeat_Length if
// necessary.
offset += line_Length;
if (offset>string_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<rng_cnt) ? len : rng_cnt;
rng_cnt -= gen_cnt;
len = gen_cnt;
while (0 != len--)
{
seed = (seed*IA + IC) % IM;
*(buf++) = seed;//This is stupid actually!
}
}
RNG_FREE();
return gen_cnt;
}
int out_tid; //Thread ID
int out_tnum = 1; //Thread number
#ifdef _OPENMP
omp_lock_t out_lock;
#define OUT_LOCK_INIT() omp_init_lock(&out_lock)
#define OUT_LOCK() omp_set_lock(&out_lock)
#define OUT_FREE() omp_unset_lock(&out_lock)
#else
#define OUT_LOCK_INIT() do{}while(0)
#define OUT_LOCK() do{}while(0)
#define OUT_FREE() do{}while(0)
#endif
static void out_init(void)
{
OUT_LOCK_INIT();
rng_tid = 0;
}
static intnative_t out_write(char * buf, intnative_t len, int curr_tid)
{
intnative_t wr_cnt = -1;//Error by default
OUT_LOCK();
if (out_tid == curr_tid)
{
if (++out_tid >= 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; i<nucl_num; i++)
{
cumul_acc += nucl_info[i].probability;
cumul_p[i] = 1ul + (uint32_t)(cumul_acc*(float)IM); //Compensate rounding errors on test file
}
#ifdef _OPENMP
intnative_t tnum = omp_get_num_procs();
if (tnum>THREADS_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;
}
+292
View File
@@ -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 <stdint.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <time.h>
#include <sys/time.h>
#include <assert.h>
#include <pthread.h>
#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; i<desired_Length_For_Oligonucleotides-1; i++)
key=(key<<2 & mask) | polynucleotide[i];
// Add all the complete oligonucleotides of
// desired_Length_For_Oligonucleotides to hash_Table and update the count
// for each oligonucleotide.
for(int i=desired_Length_For_Oligonucleotides-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);
}
}
// 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});
}
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 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; 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;
for(int i=0; i<oligonucleotide_Length; 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);
uintmax_t count= p== NULL ? 0 : *p;
snprintf(output, MAXIMUM_OUTPUT_LENGTH, "%ju\t%s", count, oligonucleotide);
dict_oligonucleotide_clear(hash_Table);
}
int main()
{
char buffer[MAXIMUM_OUTPUT_LENGTH];
worker_t worker;
worker_init(worker, 0, 1, NULL);
// Find the start of the third polynucleotide.
while(fgets(buffer, sizeof(buffer), stdin) && memcmp(">THREE", 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<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.
worker_sync_t block;
worker_start(block);
WORKER_SPAWN(worker, block, (polynucleotide, polynucleotide_Length, output), {
generate_Count_For_Oligonucleotide(polynucleotide,
polynucleotide_Length,
"GGTATTTTAATTTATAGT", output[6]);
}, ());
WORKER_SPAWN(worker, block, (polynucleotide, polynucleotide_Length, output), {
generate_Count_For_Oligonucleotide(polynucleotide,
polynucleotide_Length,
"GGTATTTTAATT", output[5]);
}, ());
WORKER_SPAWN(worker, block, (polynucleotide, polynucleotide_Length, output), {
generate_Count_For_Oligonucleotide(polynucleotide,
polynucleotide_Length,
"GGTATT", output[4]);
}, ());
WORKER_SPAWN(worker, block, (polynucleotide, polynucleotide_Length, output), {
generate_Count_For_Oligonucleotide(polynucleotide,
polynucleotide_Length,
"GGTA", output[3]);
}, ());
WORKER_SPAWN(worker, block, (polynucleotide, polynucleotide_Length, output), {
generate_Count_For_Oligonucleotide(polynucleotide,
polynucleotide_Length,
"GGT", output[2]);
}, ());
WORKER_SPAWN(worker, block, (polynucleotide, polynucleotide_Length, output), {
generate_Frequencies_For_Desired_Length_Oligonucleotides(polynucleotide,
polynucleotide_Length,
2,
output[1]);
}, ());
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++)
printf("%s\n", output_Buffer[i]);
free(polynucleotide);
worker_clear(worker);
return 0;
}