优化或建议使用omp的c++、c#代码找到所有相似的k motifs

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

给定一个正整数k (k≤50)，一个长度至多5,000的DNA串表示一个基序，一个长度至多50,000的DNA串表示一个基因组。

问题在于返回 t 的所有子串 t′ 使得编辑距离 d_E(s,t′)

The edit distance between two strings is the minimum number of elementary operations (insertions, deletions, and substitutions) to transform one string into the other, for instamce s = 'TGCATAT' and t' = 'ATCCGAT' here is a c++ implementation by user1131146 account aban Maybe it would be better to use that??

小于等于k。每个子字符串应由一对包含其在 t 中的位置后跟其长度的对进行编码。

例如

2
ACGTAG
ACGGATCGGCATCGT

应该输出

1 4
1 5
1 6

对于此示例 k=2 结果表示：
对于索引 1 to 4、d_E(s,t′)=2（ADD T 然后一个 A 在最后一个 t's G 之前）

s  = ACGTAG
t' = ACGG

对于索引 1 to 5 d_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 替换为 G 然后将索引 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;
}

我无法完成或优化它。有没有更好的办法？

如您的 post 评论中所述，为了并行调用 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 headers 定义了 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

因为我没有你的数据文件，我 运行domly 生成了一个 t 的 50,000 个 [AGCT] 字符，s 的 5,000 个，并使用了 k 共 50 次。这花费了 135 秒 ，没有命中输出。当我将 t 设置为 s 加上 45,000 个 运行dom 字符时，我得到了很多命中，对执行时间没有明显影响。

我使用 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 需要处理大量数据，因此 可能 可以 运行 这两个（或更多） ) 线程，并与线程工作的区域重叠（以便第二个线程仅在计算第一个线程结果的元素后，但在第一个线程完成对整个数组的处理之前访问它们。但是，这需要创建你自己的线程和线程之间的一些非常小心的同步，以确保正确的数据在正确的地方可用。我没有尝试让这个工作。

编辑 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 的代码与我上面列出的代码没有变化。