优化或提出c++、c#代码使用omp查找所有相似的k个motifs
Posted
技术标签:
【中文标题】优化或提出c++、c#代码使用omp查找所有相似的k个motifs【英文标题】：Optimize or propose c++, c# code using omp to find all similar k motifs 【发布时间】：2015-09-30 21:14:56 【问题描述】：给定一个正整数k (k≤50)，一个最长为s 的DNA 串5,000 代表一个基序，一个最长为t 的DNA 串50,000 代表一个基因组。
问题在于返回 t 的所有子字符串 t′ 使得编辑距离 d_E(s,t′)
两个字符串之间的编辑距离是基本的最小数量 操作（插入、删除和替换）来转换一个 串入另一个，例如
s = 'TGCATAT'和t' = 'ATCCGAT'here is a c++ implementation byuser1131146 账号被封 也许使用它会更好？？
小于或等于k。每个子字符串应由一对包含其在t 中的位置和其长度的对进行编码。
例如
2
ACGTAG
ACGGATCGGCATCGT
应该输出
1 4
1 5
1 6
对于这个例子k=2 结果意味着：
对于索引1 to 4，d_E(s,t′)=2（添加T THEN 一个A BEFORE LAST t's G）
s = ACGTAG
t' = ACGG
对于索引1 to 5d_E(s,t′)=2（将G添加到t'的末尾，然后将索引4处的t's G替换为T）
s = ACGTAG
t' = ACGGA
对于索引 1 to 6 d_E(s,t′)=2（替换最后一个 t's T BY GTHEN 将索引 4 处的 t's G 替换为 T）
s = ACGTAG
t' = ACGGAT
拥有获得all substrings of a genome that are within a certain fixed distance of the desired motif 的解决方案，使用omp 并行化解决方案的最佳方法是什么。随着字符串变得越长，程序需要太多时间。
我已经使用 omp #pragma omp parallel for 进行了测试，然后在写入文件部分使用了 lock 以及 #pragma omp critical 但是我不知道我是否正确并行化它。
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
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
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;
我无法完成或优化它。有没有更好的办法？
【问题讨论】：
在内循环中做file io不是个好主意 您无法按照您的方式并行化a[0][j]=a[0][j-1]-1。如果我理解正确，您可以将其重写为a[0][j] = a[0][0]-j，这样可以很好地并行化。例如 a[] = 9,3,4,5 在循环之后是 9,8,7,6。
您的代码还有其他问题。在您的函数alignment 中，您写入每个线程共享的a。这将导致竞争条件。不仅如此，您写入它的方式还依赖于i 和j 上的先前迭代。您需要删除依赖项，或者尽可能使用前缀和之类的东西。
那么更好的方法是什么？
【参考方案1】：
如您帖子的评论中所述，为了并行化对alignment 的调用，每个线程都需要拥有自己的a 副本。您可以使用firstprivate OpenMP 子句做到这一点：
#pragma omp parallel for firstprivate(a)
在alignment 本身中，您在循环中重复计算，优化器可能不会消除这些计算。在循环条件下对a.size 和使用min 的调用可以计算一次并存储在局部变量中。不断计算a[i]和a[i-1]也可以提升出内循环。
int asize = int(a.size());
for(i = 1; i < asize; i++)
int jend = min(ylen, i + k);
vi &cura = a[i];
vi &preva = a[i-1];
for(j = max(1, i - 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));
Windows 标头定义 max 和 min 宏。如果您使用那些（而不是 STL 中的内联函数）也可能会不必要地重复代码。
另一个瓶颈可能是匹配的输出，这取决于找到匹配的频率。改善这一点的一种方法是将i-1,j 对存储到一个新的局部变量中（也将其包含在private 子句中），然后使用reduction 子句来组合结果（尽管我不确定如何适用于容器）。完成 for 循环后，您可以输出结果，可能先对它们进行排序。
尝试并行化初始化a[0] 的j 循环可能不值得。代码需要修复才能使其正常工作（也在评论中提到），如果启动线程的开销过多，或者如果多个线程尝试在线程之间存在缓存行争用，则多个线程可能会导致其运行速度变慢写入内存中几乎相邻的值。如果您的代码中的示例 s 是典型的，我会在一个线程上运行它。
当你构造你的a向量时，你可以在构造函数中包含-999的初始值，如果包含向量：
vector<vi> a(slen + 1, vi(slen + k + 1, -999));
那么它下面的初始化循环只需要设置每个包含向量中第一个元素的值。
首先应该这样做。注意：代码示例是建议，未经编译或测试。
编辑我已经解决了这个问题，结果不是你所希望的。
最初我只是运行您的代码（以获取我的基准性能数据）。使用VC2010，最大优化，64位编译。
cl /nologo /W4 /MD /EHsc /Ox /favor:INTEL64 sequencing.cpp
由于我没有你的数据文件，我随机生成了一个包含 50,000 个 [AGCT] 字符的 t，s 的 5,000 个字符，并使用了一个 50 个的 k。这花了 135 秒 没有命中输出。当我将t 设置为s 加上45,000 个随机字符时，我得到了很多命中，而对执行时间没有明显影响。
我使用 omp 进行了这项工作（使用 firstprivate 而不是 private 来获取复制的 a 的初始值）并且它崩溃了。仔细观察，我意识到对alignment 的调用取决于上一次调用的结果。所以这不能在多核上执行。
我确实重写了 alignment 以删除所有冗余计算，并减少了大约 25% 的执行时间（103 秒）：
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; i < asize; i++)
auto xi = x[i - 1];
int jend = min(ylen, i + k);
vi &cura = a[i];
const vi &preva = a[i-1];
j = max(1, i - 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;
我做的最后一件事是使用 VS2015 编译它。这将执行时间又减少了 28% 至大约 74 秒。
可以在main 中进行类似的调整和调整，但对性能没有明显影响，因为大部分时间都花在了alignment 中。
使用 32 位二进制文件的执行时间相似。
我确实想到，由于 alignment 处理大量数据，可能可以在两个（或更多）线程上运行它，并重叠线程使用的区域（以便第二个线程仅在计算第一个线程结果的元素后才访问它们，但在第一个线程完成对整个数组的处理之前。但是，这需要创建自己的线程并在线程之间进行一些非常仔细的同步以确保在正确的位置提供正确的数据。我没有尝试完成这项工作。
EDIT 2 优化版 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;
alignment 的代码与我上面列出的相同。
【讨论】：
如果您可以测试程序以减少其执行时间，那么赏金就是您的 @cMinor 这个周末我去看看，昨天没空。 改进aligment的好方法，但是，你认为有没有更好的方法来加速主函数，特别是调用对齐的部分？ 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++) ....我正在用这个文件测试程序，你能告诉你在你的电脑上完成需要多少时间吗？dropbox.com/s/vb9wjho6mp8lvgf/test_k.txt?dl=0。最后，VS 2015 对代码有这么大的改进吗？为什么？
@cMinor VS2010，94 秒。 VS2015，67 秒。 5318 场比赛。查看两个编译器的反汇编，看起来 VS2015 在指令调度方面做得更好，尽管指令选择也略有不同。对于主函数，您可以将与j 进行比较的min 分配给局部变量，并将a[slen][j] 更改为指针或迭代器，就像我在alignment 中所做的那样。不过，这不会产生太大影响，因为大部分时间都花在alignment。
好的，您能否分享您使用的所有代码（主要和对齐功能），以便我可以奖励您？以上是关于优化或提出c++、c#代码使用omp查找所有相似的k个motifs的主要内容，如果未能解决你的问题，请参考以下文章