Optimize or propose c++, c# code using omp to find all similar k motifs -

given positive integer k (k≤50), dna string s of length @ 5,000 representing motif, , dna string t of length @ 50,000 representing genome.

the problem consists on returning substrings t′ of t such edit distance d_e(s,t′)

the edit distance between 2 strings minimum number of elementary operations (insertions, deletions, , substitutions) transform 1 string other, instamce s = 'tgcatat' , t' = 'atccgat' here c++ implementation by user1131146 account aban maybe better use that??

is less or equal k. each substring should encoded pair containing location in t followed length.

for instance

2 acgtag acggatcggcatcgt 

should output

1 4 1 5 1 6 

for example k=2 results mean:
indices 1 4, d_e(s,t′)=2 (add t 1 a before last t's g)

s  = acgtag t' = acgg 

for indices 1 5 d_e(s,t′)=2 (add g end of t' replace t's g @ index 4 t)

s  = acgtag t' = acgga 

for indices 1 6 d_e(s,t′)=2 (replace last t's t gthen replace t's g @ index 4 t)

s  = acgtag t' = acggat 

having solution all substrings of genome within fixed distance of desired motif,what best way parallelize solution using omp. longer strings become program takes time.

i have tested using omp #pragma omp parallel for using lock in write file section, , #pragma omp critical not know if paralelizing correctly.

void alignment(vector<vi>&a, string &x, string y, int k){     string tx,ty;     int i,j;     int ylen=a[0].size();     for(i=1;i<a.size();i++){         for(j=max(1,i-k);j<=min(ylen,i+k);j++){              a[i][j] = max(x[i-1] == y[j-1]?a[i-1][j-1] : (a[i-1][j-1]-1), max(a[i-1][j]-1,a[i][j-1]-1));         }     } }  int main() {     int k = 23;     string s = "aattagctaaggtgtacgatgtcccattgtgtaagattaggaactccatttaggttacctccgtcttaagtgatggaccgtgggtagctgcgtccgatggactcatgcagcgcccggatacctgcagtatttattatataggtctgcgcaccaaacgatttctttcgtggtcgggattccggggtcctccgctattcagagagctaaata";     string t = "acaatgcagcaatccagcgccggaatttaagaataggtcagtttgtaaggcactgttcccgttattcgtaatgcagtattaacgttaatgctcgagaccatattggacgtcagtatgcagacctgtgctagggtggtctatttcaagatcaccgagctaggcgcgtgagctaacaggccgtaatggtggcgcccgctcccataatcacttcacgaagcattaggtagactaccatttaggaagccctctcgcccgcgtactggttacagcccactacaatggatactccttacttcggtgcaggcaagacttctacaaagaagcgtccaagaagttgtcgtagctcgttcttaccccacctgtataaaattgatccagtcgtacatatgacgatgctgagcctcggactggtaaatacaagtcaaaggaccaacccattacagtatgaactaccggtgg";     time_t start = time(null);     std::ofstream out("output.txt");     ifstream somestream( "data.txt" );     string line;     getline( somestream, line );    int k = atoi(line.c_str() );     getline( somestream, line );    string s =line;     getline( somestream, line );    string t= line;     int slen=s.length(), tlen=t.length();     vector<vi>a( slen+1, vi(slen+k+1));      int i,j;     for(i=1;i<a.size();i++)         fill(a[i].begin(),a[i].end(),-999),a[i][0]=a[i-1][0]-1;     #pragma omp parallel     {         for(j=1;j<a[0].size();j++)         {             a[0][j]=a[0][j-1]-1;         }     }     //cout << "init";     time_t endinit = time(null);     cout<<"execution init time: "<< (double)(endinit-start)<<" seconds"<<std::endl;     //omp_lock_t writelock;     //omp_init_lock(&writelock);      #pragma omp parallel     {     for(i=0;i<=tlen-slen+k;i++)     {         alignment(a,s,t.substr(i,slen+k),k);          for(j=max(0,slen-k);j<=min(slen+k,tlen-i);j++)         {             if(a[slen][j]>=-k)             {                  //omp_set_lock(&writelock);                     //cout<<(i+1)<<' '<<j<<endl;                 #pragma omp critical                  {                     out <<(i+1)<<' '<<j<<endl;                 }                 //omp_unset_lock(&writelock);             }         }     }     }     //omp_destroy_lock(&writelock);     time_t end = time(null);     cout<<"execution time: "<< (double)(end-start)<<" seconds"<<std::endl;       out.close();     return 0; } 

i have not been able complete or optimize it. there better way?

as mentioned in comment post, in order parallelize calls alignment, each thread needs have own copy of a. can firstprivate openmp clause:

#pragma omp parallel firstprivate(a) 

in alignment repeat calculations in loops optimizer might not eliminating. calls a.size , using min in conditions of loops computed once , stored in local variable. calculating a[i] , a[i-1] can lifted out of inner loop.

int asize = int(a.size()); for(i = 1; < asize; i++) {     int jend = min(ylen, + k);     vi &cura = a[i];     vi &preva = a[i-1];     for(j = max(1, - k);j <= jend; j++)         cura[j] = max(x[i-1] == y[j-1]?preva[j-1] : (preva[j-1]-1), max(preva[j]-1,cura[j-1]-1)); } 

the windows headers define max , min macros. if you're using (and not inline functions in stl) repeat code unnecessarily.

another bottleneck output of matches, depending on how match found. 1 way improve store i-1,j pairs new local variable (also include in private clause), use reduction clause combine results (although i'm not sure how works containers). once you're done loop can output results, possibly sorting them first.

trying parallelize j loop initializes a[0] not worthwhile. code needs fixed work (also mentioned in comment), , multiple threads cause run slower if overhead starting threads much, or if there cache line contention among threads if multiple threads try write adjacent values in memory. if sample s in code typical i'd run on 1 thread.

when construct a vector, can include -999 initial value in constructor if contained vector:

vector<vi> a(slen + 1, vi(slen + k + 1, -999)); 

then initialization loop right below have set value of first element in each contained vector.

that should start. note: code samples suggestions, , have not been compiled or tested.

edit i've worked thru this, , results not hoping for.

initially ran code (to baseline performance numbers). using vc2010, maximum optimization, 64 bit compile.

cl /nologo /w4 /md /ehsc /ox /favor:intel64 sequencing.cpp 

since don't have data files, randomly generated t of 50,000 [agct] characters, s of 5,000, , used k of 50. took 135 seconds no hits output. when set t s plus 45,000 random characters, got lot of hits no noticeable impact on execution time.

i got working omp (using firstprivate instead of private initial value of a copied in) , crashed. looking in realized calls alignment depend on result of previous call. cannot executed on multiple cores.

i did rewrite alignment remove redundant computations , got reduction of around 25% in execution time (to 103 seconds):

void alignment(vector<vi>&a, const string &x, const string y, int k){     string tx,ty;     int i,j;     int ylen=int(a[0].size());     int asize = int(a.size());     for(i = 1; < asize; i++) {         auto xi = x[i - 1];         int jend = min(ylen, + k);         vi &cura = a[i];         const vi &preva = a[i-1];         j = max(1, - k);         auto caj = cura[j - 1];         const auto *yp = &y[j - 1];         auto *pca = &cura[j];         auto *ppj = &preva[j - 1];         for(;j <= jend; j++) {             caj = *pca = max(*ppj - (xi != *yp), max(ppj[1] - 1, caj - 1));             ++yp;             ++pca;             ++ppj;         }     } } 

the last thing did compile using vs2015. reduced execution time 28% around 74 seconds.

similar tweaks , adjustments can made in main, no observable effect on performance majority of time spent in alignment.

execution times using 32 bit binary similar.

it did occur me since alignment works lot of data, might possible run on 2 (or more) threads, , overlap area threads work (so second thread accessing elements of first thread's results once computed, before first thread has completed working entire array. however, require creating own thread , careful synchronization between thread ensure right data available in right place. didn't try make work.

edit 2 source optimized version of main:

int main() { //  int k = 50; //  string s = "aattagctaaggtgtacgatgtcccattgtgtaagattaggaactccatttaggttacctccgtcttaagtgatggaccgtgggtagctgcgtccgatggactcatgcagcgcccggatacctgcagtatttattatataggtctgcgcaccaaacgatttctttcgtggtcgggattccggggtcctccgctattcagagagctaaata"; //  string t = "acaatgcagcaatccagcgccggaatttaagaataggtcagtttgtaaggcactgttcccgttattcgtaatgcagtattaacgttaatgctcgagaccatattggacgtcagtatgcagacctgtgctagggtggtctatttcaagatcaccgagctaggcgcgtgagctaacaggccgtaatggtggcgcccgctcccataatcacttcacgaagcattaggtagactaccatttaggaagccctctcgcccgcgtactggttacagcccactacaatggatactccttacttcggtgcaggcaagacttctacaaagaagcgtccaagaagttgtcgtagctcgttcttaccccacctgtataaaattgatccagtcgtacatatgacgatgctgagcctcggactggtaaatacaagtcaaaggaccaacccattacagtatgaactaccggtgg";     time_t start = time(null);     std::ofstream out("output.txt");     ifstream somestream( "data.txt" );     string line;     getline( somestream, line );    int k = atoi(line.c_str() );     getline( somestream, line );    string s =line;     getline( somestream, line );    string t= line;     int slen=int(s.length()), tlen=int(t.length());      int i,j;     vector<vi> a(slen + 1, vi(slen + k + 1, -999));     a[0][0]=0;     for(i=1;i<=slen;i++)         a[i][0]=-i;      {         int ej=int(a[0].size());         for(j=1;j<ej;j++)             a[0][j] = -j;     }     //cout << "init";     time_t endinit = time(null);     cout<<"execution init time: "<< (double)(endinit-start)<<" seconds"<<std::endl;      {         int endi=tlen-slen+k;         for(i=0;i<=endi;i++)         {             alignment(a,s,t.substr(i,slen+k),k);             int ej=min(slen+k,tlen-i);             j=max(0,slen-k);             const auto *aj = &a[slen][j];             for(;j<=ej;j++,++aj)             {                 if(*aj>=-k)                 {                     //cout<<(i+1)<<' '<<j<<endl;                     out <<(i+1)<<' '<<j<<endl;                 }             }         }     }     time_t end = time(null);     cout<<"execution time: "<< (double)(end-start)<<" seconds"<<std::endl;       out.close();     return 0; } 

the code alignment unchanged listed above.