Skip to content

Commit

Permalink
Re upload of every file
Browse files Browse the repository at this point in the history
  • Loading branch information
@GambarottoFilippo
GambarottoFilippo authored Aug 28, 2017
1 parent cf07259 commit d316cd6
Show file tree
Hide file tree
Showing 10 changed files with 736 additions and 66 deletions.
10 changes: 4 additions & 6 deletions ParallelGeneralized/Makefile
View file
Original file line number Diff line number Diff line change
@@ -1,20 +1,18 @@
OBJS = SLT.c SLT_MAWs.c dbwt_queue.c indexed_DNA5_seq.c DNA5_tables.c dbwt.c dbwt_utils.c mt19937ar.c DNA5_Basic_BWT.c sais.c ../malloc_count-master/malloc_count.c ../malloc_count-master/stack_count.c
OBJNAIVE= naive_MAWs.c
OBJS = SLT.c SLT_MAWs.c SLT_single_string.c SLT_MAWs_single_string.c dbwt_queue.c indexed_DNA5_seq.c DNA5_tables.c dbwt.c dbwt_utils.c mt19937ar.c DNA5_Basic_BWT.c sais.c ../malloc_count-master/malloc_count.c ../malloc_count-master/stack_count.c naive_MAWs.c

HDRS = SLT.h dbwt_queue.h indexed_DNA5_seq.h dbwt.h dbwt_utils.h mt19937ar.h DNA5_Basic_BWT.h SLT_MAWs.h ../malloc_count-master/malloc_count.h ../malloc_count-master/stack_count.h
HDRS = SLT.h dbwt_queue.h indexed_DNA5_seq.h dbwt.h dbwt_utils.h mt19937ar.h DNA5_Basic_BWT.h SLT_MAWs.h ../malloc_count-master/malloc_count.h ../malloc_count-master/stack_count.h naive_MAWs.h

HDRNAIVE= naive_MAWs.h




LIBS = -ldl
LIBS = -ldl -lm
#CFLAGS = -g -Wall -O2 -fopenmp
CFLAGS = -Wall -O3 -fopenmp

CC = gcc

test_SLT_MAWs : $(OBJS) test_SLT_MAWs.o ../malloc_count-master/malloc_count.o ../malloc_count-master/stack_count.o $(HDRS)
test_SLT_MAWs : $(OBJS) test_SLT_MAWs.o ../malloc_count-master/malloc_count.o ../malloc_count-master/stack_count.o $(HDRS)
$(CC) $(CFLAGS) $(OBJS) test_SLT_MAWs.o -o test_SLT_MAWs $(LIBS)
test_SLT_MAWs.o: test_SLT_MAWs.c $(HDRS)
$(CC) $(CFLAGS) -c test_SLT_MAWs.c
Expand Down
217 changes: 187 additions & 30 deletions ParallelGeneralized/SLT_MAWs.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,8 @@
#define alloc_growth_denom 3
static unsigned int length1;
static unsigned int length2;
static unsigned int F1;
static unsigned int F2;


typedef struct
Expand Down Expand Up @@ -81,8 +83,10 @@ void SLT_callback_MAWs(const SLT_joint_params_t * SLT_joint_params,void * intern
double x=(double)1/((SLT_joint_params->string_depth+2)*(SLT_joint_params->string_depth+2));
char_mask1=1;
for(i=1;i<5;i++) {
char_mask1<<=1;
char_mask2=1;
for(j=1;j<5;j++) {
char_mask2<<=1;
if((SLT_joint_params->right_extension_bitmap1&char_mask2)
&& (SLT_joint_params->left_extension_bitmap1&char_mask1)
&& SLT_joint_params->left_right_extension_freqs1[i][j]==0) {
Expand All @@ -106,6 +110,87 @@ void SLT_callback_MAWs(const SLT_joint_params_t * SLT_joint_params,void * intern
// We have a MAW. Write it to the output
state->nMAWs++;
state->LW-=2*x;
if(mem) {
if(state->MAW_buffer_idx+SLT_joint_params->string_depth+3>state->nMAW_capacity) {
#pragma omp critical
fwrite(state->MAW_buffer, state->MAW_buffer_idx, sizeof(char), state->file);
state->MAW_buffer_idx=0;
}
state->MAW_buffer[state->MAW_buffer_idx++]=alpha4_to_ACGT[i-1];
for(k=0;k<SLT_joint_params->string_depth;k++){
if (state->char_stack[SLT_joint_params->string_depth-k-1]-1 > 3)
printf("%d ", state->char_stack[SLT_joint_params->string_depth-k-1]-1);
state->MAW_buffer[state->MAW_buffer_idx++]=
alpha4_to_ACGT[state->char_stack[SLT_joint_params->string_depth-k-1]-1];
}
state->MAW_buffer[state->MAW_buffer_idx++]=alpha4_to_ACGT[j-1];
state->MAW_buffer[state->MAW_buffer_idx++]='\n';
}
}
}
}

}
void SLT_callback_RWs(const SLT_joint_params_t * SLT_joint_params,void * intern_state, unsigned int mem)
{
MAWs_callback_state_t * state= (MAWs_callback_state_t*)(intern_state);
unsigned int char_mask1;
unsigned int char_mask2;
unsigned int i,j,k,h;

if(state->nMAW_capacity==0 && mem) {
state->nMAW_capacity=1<<16;
state->MAW_buffer=(unsigned char *) malloc(state->nMAW_capacity);
}
if(SLT_joint_params->string_depth!=0)
{
if(SLT_joint_params->string_depth>state->char_stack_capacity)
{
state->char_stack_capacity=(state->char_stack_capacity+1)*alloc_growth_num/alloc_growth_denom;
state->char_stack=(unsigned char *)realloc(state->char_stack,state->char_stack_capacity);
};
state->char_stack[SLT_joint_params->string_depth-1]=SLT_joint_params->WL_char;
}

// Check that we are at a maximal repeat of length at least minlen-2
if(SLT_joint_params->string_depth + 2 < state->minlen ||
((SLT_joint_params->nleft_extensions1 < 2 || SLT_joint_params->nright_extensions1 < 2) &&
(SLT_joint_params->nleft_extensions2 < 2 || SLT_joint_params->nright_extensions2 < 2)))
return;

char_mask1=1;
for(i=1;i<5;i++) {
char_mask1<<=1;
char_mask2=1;
for(j=1;j<5;j++) {
char_mask2<<=1;
unsigned int freqLeft1= 0;
unsigned int freqRight1= 0;
unsigned int freqLeft2= 0;
unsigned int freqRight2= 0;
for(k=1; k<5; k++) {
freqLeft1+= SLT_joint_params->left_right_extension_freqs1[i][k];
freqLeft2+= SLT_joint_params->left_right_extension_freqs2[i][k];
freqRight1+= SLT_joint_params->left_right_extension_freqs1[k][j];
freqRight2+= SLT_joint_params->left_right_extension_freqs2[k][j];
}

if(freqLeft1 >= F2 && freqRight1 >= F2
&& SLT_joint_params->left_right_extension_freqs1[i][j]<=F1) {
// We have a MAW. Write it to the output
state->nMAWs1++;
}
if(freqLeft2 >= F2 && freqRight2 >= F2
&& SLT_joint_params->left_right_extension_freqs2[i][j]<=F1) {
// We have a MAW. Write it to the output
state->nMAWs2++;
}
if(freqLeft1 >= F2 && freqRight1 >= F2
&& freqLeft2 >= F2 && freqRight2 >= F2
&& SLT_joint_params->left_right_extension_freqs1[i][j]<=F1
&& SLT_joint_params->left_right_extension_freqs2[i][j]<=F1) {
// We have a MAW. Write it to the output
state->nMAWs++;
if(mem) {
if(state->MAW_buffer_idx+SLT_joint_params->string_depth+3>state->nMAW_capacity) {
#pragma omp critical
Expand All @@ -120,9 +205,7 @@ void SLT_callback_MAWs(const SLT_joint_params_t * SLT_joint_params,void * intern
state->MAW_buffer[state->MAW_buffer_idx++]='\n';
}
}
char_mask2<<=1;
}
char_mask1<<=1;
}

}
Expand Down Expand Up @@ -194,8 +277,6 @@ void SLT_callback_kernel(const SLT_joint_params_t * SLT_joint_params,void * inte
int N=0;
int d1=0;
int d2=0;
mem=1;
unsigned int maw_pres=0; //if maw in t1 and present in t2

if(state->nMAW_capacity==0 && mem) {
state->nMAW_capacity=1<<16;
Expand Down Expand Up @@ -421,6 +502,32 @@ unsigned int SLT_find_MAWs(Basic_BWT_t * BBWT1,Basic_BWT_t * BBWT2,
f=fopen("output.txt", "a");
state.file=f;
}
length1= BBWT1->textlen+2;
length2= BBWT2->textlen+2;
state.N=0;
state.D1=0;
state.D2=0;

//Initializing D1 and D2
double prefix_sum= 0;
for(i=1; i<=BBWT1->textlen+2;i++) {
prefix_sum+= (g1(i)-1)*(g1(i)-1);
state.D1+=prefix_sum;
}
prefix_sum= 0;
for(i=1; i<=BBWT2->textlen+2;i++){
prefix_sum+= (g2(i)-1)*(g2(i)-1);
state.D2+=prefix_sum;
}

//Prefix sum allocation and initialization
state.prefix_capacity=4;
state.prefix_sum1=(double *)malloc(state.prefix_capacity*sizeof(double));
state.prefix_sum2=(double *)malloc(state.prefix_capacity*sizeof(double));
state.prefix_sumN=(double *)malloc(state.prefix_capacity*sizeof(double));
state.prefix_sum1[1]=(g1(1)-1)*(g1(1)-1);
state.prefix_sum2[1]=(g2(1)-1)*(g2(1)-1);
state.prefix_sumN[1]=(g1(1)-1)*(g2(1)-1);
//resize!
/*if(d> state.bitvec_sizes[i]){
state.bitvec_sizes[i]*=2;
Expand Down Expand Up @@ -459,32 +566,6 @@ unsigned int SLT_find_MAWs(Basic_BWT_t * BBWT1,Basic_BWT_t * BBWT2,
*_output_result= state.LW;
break;
case 6:
length1= BBWT1->textlen+2;
length2= BBWT2->textlen+2;
state.N=0;
state.D1=0;
state.D2=0;

//Initializing D1 and D2
double prefix_sum= 0;
for(i=1; i<=BBWT1->textlen+2;i++) {
prefix_sum+= (g1(i)-1)*(g1(i)-1);
state.D1+=prefix_sum;
}
prefix_sum= 0;
for(i=1; i<=BBWT2->textlen+2;i++){
prefix_sum+= (g2(i)-1)*(g2(i)-1);
state.D2+=prefix_sum;
}

//Prefix sum allocation and initialization
state.prefix_capacity=4;
state.prefix_sum1=(double *)malloc(state.prefix_capacity*sizeof(double));
state.prefix_sum2=(double *)malloc(state.prefix_capacity*sizeof(double));
state.prefix_sumN=(double *)malloc(state.prefix_capacity*sizeof(double));
state.prefix_sum1[1]=(g1(1)-1)*(g1(1)-1);
state.prefix_sum2[1]=(g2(1)-1)*(g2(1)-1);
state.prefix_sumN[1]=(g1(1)-1)*(g2(1)-1);

SLT_iterator=new_SLT_joint_iterator(SLT_callback_kernel,SLT_cloner, SLT_combiner,SLT_free,&state,BBWT1,BBWT2,SLT_stack_trick, mem, cores);
SLT_joint_execute_iterator(SLT_iterator);
Expand All @@ -502,6 +583,82 @@ unsigned int SLT_find_MAWs(Basic_BWT_t * BBWT1,Basic_BWT_t * BBWT2,
return state.nMAWs;
};


unsigned int SLT_find_RWs(Basic_BWT_t * BBWT1,Basic_BWT_t * BBWT2,
unsigned int minlen, unsigned int * _nMAWs1,
unsigned int * _nMAWs2, double * _output_result, unsigned int mem, unsigned int cores, unsigned int result, unsigned int f1, unsigned int f2)
{
SLT_joint_iterator_t * SLT_iterator;
MAWs_callback_state_t state;

unsigned int i;
state.nMAWs=0;
state.nMAWs1=0;
state.nMAWs2=0;
state.MAW_buffer=0;
state.MAW_buffer_idx=0;
state.nMAW_capacity=0;
state.minlen=minlen;
state.LW=0;
state.char_stack_capacity=4;
state.char_stack=(unsigned char *) malloc(state.char_stack_capacity);
FILE *f;
if(mem) {
f=fopen("output.txt", "a");
state.file=f;
}
length1= BBWT1->textlen+2;
length2= BBWT2->textlen+2;
F1= f1;
F2= f2;
state.N=0;
state.D1=0;
state.D2=0;

//Initializing D1 and D2
double prefix_sum= 0;
for(i=1; i<=BBWT1->textlen+2;i++) {
prefix_sum+= (g1(i)-1)*(g1(i)-1);
state.D1+=prefix_sum;
}
prefix_sum= 0;
for(i=1; i<=BBWT2->textlen+2;i++){
prefix_sum+= (g2(i)-1)*(g2(i)-1);
state.D2+=prefix_sum;
}

//Prefix sum allocation and initialization
state.prefix_capacity=4;
state.prefix_sum1=(double *)malloc(state.prefix_capacity*sizeof(double));
state.prefix_sum2=(double *)malloc(state.prefix_capacity*sizeof(double));
state.prefix_sumN=(double *)malloc(state.prefix_capacity*sizeof(double));
state.prefix_sum1[1]=(g1(1)-1)*(g1(1)-1);
state.prefix_sum2[1]=(g2(1)-1)*(g2(1)-1);
state.prefix_sumN[1]=(g1(1)-1)*(g2(1)-1);
//resize!
/*if(d> state.bitvec_sizes[i]){
state.bitvec_sizes[i]*=2;
state.bitvec[i]=realloc(state.bitvec[i],state.bitvec_sizes[i]/8);
memset(state.bitvec[i]+state.bitvec_sizes[i]/16,0,state.bitvec_sizes[i]/16);
state.bitvec_capacity=(unsigned int *) malloc(36*sizeof(unsigned int));
}
for(i=0;i<36;i++) state.bitvec_capacity[i]=32;
state.bitvec1= (unsigned int **)malloc(36*sizeof(int));
for(i=0;i<36;i++) {
for(i=0;i<36;i++) state.bitvec1[i]=(unsigned int *) calloc(1,sizeof(unsigned int));
}
state.KL_capacity= BBWT1->textlen+1;
state.KL= (double *) malloc(state.KL_capacity*sizeof(double));
*/

SLT_iterator=new_SLT_joint_iterator(SLT_callback_RWs,SLT_cloner, SLT_combiner,SLT_free,&state,BBWT1,BBWT2,SLT_stack_trick, mem, cores);
SLT_joint_execute_iterator(SLT_iterator);

*_nMAWs1=state.nMAWs1;
*_nMAWs2=state.nMAWs2;
return state.nMAWs;
};

void convert_MAWs_to_ACGT(unsigned char ** MAW_ptr,unsigned int nMAWs)
{
unsigned int i;
Expand Down
2 changes: 2 additions & 0 deletions ParallelGeneralized/SLT_MAWs.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,9 +2,11 @@
#define SLT_MAWs_h
#include"SLT.h"

void SLT_callback_MAWs_single_string(const SLT_joint_params_t * SLT_params,void * intern_state, unsigned int memory);
void SLT_callback_MAWs(const SLT_joint_params_t * SLT_params,void * intern_state, unsigned int memory);
void SLT_callback_MAWs_present(const SLT_joint_params_t * SLT_params,void * intern_state, unsigned int memory);
void SLT_callback_MAWs_kernel(const SLT_joint_params_t * SLT_params,void * intern_state, unsigned int memory);
void SLT_callback_RWs(const SLT_joint_params_t * SLT_params,void * intern_state, unsigned int memory);
void* SLT_cloner(void* p, unsigned int t);
void SLT_combiner(void** intern_state, void* state, unsigned int t, unsigned int mem);
void SLT_free(void* intern_state, unsigned int mem);
Expand Down
Loading

0 comments on commit d316cd6

Please sign in to comment.