c++ - 如何优化这种组合算法？

Question

我正在编写一个分子动力学程序，需要获取分子中的原子并找到它们可能的键合方式。为此，我有一个 Atom 对象 vector ，并使用以下算法生成组合对:

    void CombinationKN(std::vector<std::vector<int>> &indices, int K, int N) {
        std::string bitmask(K, 1);
        bitmask.resize(N, 0);

        do {
            /* This loop takes forever with larger N values (approx. 3000) */
            std::vector<int> indexRow;

            for (int i = 0; i < N; i++)
            {
                if (bitmask[i]) indexRow.push_back(i);
            }

            indices.push_back(indexRow);
        } while (std::prev_permutation(bitmask.begin(), bitmask.end()));
    }

这是一个简单的 N 选择 K 算法(即返回的索引可以包含 (1, 2) 但不包含 (2, 1))，在我的例子中 N 是分子中的原子数，K 是 2。

然后我这样调用算法:

void CalculateBondGraph(const std::vector<Atom *> &atoms, std::map<int, 
    std::map<int, double>> &bondGraph, ForceField *forceField) {
    int natoms = atoms.size();

    std::vector<std::vector<int>> indices;

    utils::CombinationKN(indices, 2, natoms);

    for (auto &v : indices) {
        int i = v[0];
        int j = v[1];

        /*... Check if atoms i and j are bonded based on their coordinates.*/
    }
}

此算法的问题在于，对于具有 3000 多个原子的大分子，它需要很长时间才能完成。我考虑过将其并行化(特别是使用 OpenMP)，但即便如此，工作也必须在几个线程之间进行拆分，并且仍然需要很多时间才能完成。我需要一种方法来优化此算法，以便计算组合不会花费很长时间。感谢您的帮助。

谢谢，维卡斯

Answer 1

如果 K 远小于 N 并且如果 N 很大，当然 K 比 N 小得多，否则你会很快耗尽内存。

请注意，每个有效的 index_row 都是 K 小于 N 的整数的严格单调递增序列，反之亦然。直接生成这些很容易:

void CombinationKN(std::vector<std::vector<int>> &indices, int K, int N) {
    std::vector<int> index_row;
    // lexographically first valid row
    for (int i=0; i<K; ++i) {
        index_row.push_back(i);
    }

    for(;;) {
        // output current row
        indeces.push_back(index_row);

        // increment index_row the the lexically next valid sequence
        // find the right-most index we can increment
        // This loop does O(1) amortized iterations if K is not large.  O(K) worst case
        int inc_index=K-1;
        int index_limit=N-1;
        while(inc_index >= 0 && index_row[inc_index] >= index_limit) {
          --inc_index;
          --index_limit;
        }
        if (inc_index < 0) {
            break; //all done
        }
        // generate the lexically first valid row with matching prefix and
        // larger value at inc_index
        int val = index_row[inc_index]+1;
        for (;inc_index<K; ++inc_index, ++val) {
            index_row[inc_index] = val;
        }
    }
}

此外，如果您对这些组合所做的唯一事情就是遍历它们，那么就没有理由浪费存储它们的整个列表所需的(可能非常大量的)内存。上面的函数包含一个在您需要时从上一个生成下一个的过程。