在C++中寻找一个经过网格/矩阵的成本最优路径

Question

在C++中寻找一个经过网格/矩阵的成本最优路径

c++algorithmgraphpermutationshortest-path

13

我遇到了一个难题，在线上没有找到太多有用的帮助。我需要找到由多个数字向量中组成的最小成本组合。所有向量的大小都相同。例如，考虑以下情况：

row [0]:  a  b  c  d   
row [1]:  e  f  g  h  
row [2]:  i  j  k  l

现在我需要从每一行即向量中取一个元素来找到数字的组合，例如：aei
接下来，我需要找到另外三个不相交的组合，例如：bfj，cgk，dhl。我基于这四个选定的组合计算成本。目标是找到给出最小成本的组合。另一个可能的组合可以是：afj，bei，chk，dgl。如果总列数为d，行数为k，则可能的总组合数为d^k。行被存储为向量。我现在遇到了困难，难以编写上述过程的算法。如果有人能帮忙，我将不胜感激。
谢谢。

// I am still working on the algorithm. I just have the vectors and the cost function.  

//Cost Function  , it also depends on the path chosen
float cost(int a, int b, PATH to_a) {  
float costValue;  
...  
...  
return costValue;  
}  

vector< vector < int > > row;  
//populate row  
...   
...
//Suppose  

//    row [0]:  a  b  c  d   
//    row [1]:  e  f  g  h  
//    row [2]:  i  j  k  l   

// If a is chosen from row[0] and e is chosen from row[1] then,  
float subCost1 = cost(a,e, path_to_a);  

// If i is chosen from row[2] ,  
float subCost2 = cost(e,i,path_to_e);  

// Cost for selecting aei combination is  
float cost1 = subCost1 + subCost2;  

//similarly other three costs need to be calculated by selecting other remaining elements  
//The elements should not intersect with each other eg. combinations aei and bej cannot exist on the same set.  

//Suppose the other combinations chosen are bfj with cost cost2, cgk with cost cost3   and dhl with cost cost4  
float totalCost = cost1 + cost2 + cost3 + cost4;   

//This is the cost got from one combination. All the other possible combinations should be enumerated to get the minimum cost combination.

- coolcav

1

如果您只是在寻找C++的答案，正如标题所示，为什么还要标记为C呢？ - Flexo

这是一道作业题吗？如果是的话，应该打上相应的标签。 - GWW

1

如果它们被存储为向量，则“C”标签可能不适用。 - Mooing Duck

@MStodd 我的意思是多个数字向量。我现在正在编辑它。 - coolcav

@coolcav：我看到你把成本函数变成了动态函数。我们可以假设两个给定节点之间的成本是恒定的吗（即成本函数是确定性的）？ - sehe

显示剩余8条评论

5个回答

6

首先，有一个稍微有帮助的观察结果。

我认为4！^3的结果没有体现出{ aei, bfj, cgk, dhl }和（例如）{ bfj, aei, cgk, dhl }具有相同的成本。

这意味着我们只需要考虑形式为的序列。

{ a??, b??, c??, d?? }

这种等价关系可以将不同情况的数量减少4倍！另一方面，@sehe认为3x4等于4!^3（我同意），因此类似地，6x48需要48!^6。其中“仅有”的是48!^5个不同情况。现在这个数字是2.95 × 10^305。使用3x4的例子，以下是一个算法的开始，可以给出某种答案。

Enumerate all the triplets and their costs. 
Pick the lowest cost triplet.
Remove all remaining triplets containing a letter from that triplet.
Now find the lowest cost triplet remaining.
And so on.

请注意，这不是完整的穷举搜索。

我从中发现这仍然需要大量的计算。第一次搜索仍需计算48^6（12,230,590,464）次成本。我想这可以做到，但需要付出很多努力。相比之下，后续的搜索成本将会更加便宜。

- Keith

你的想法和我一模一样。不过，我本来要“问”一下楼主成本函数是否只对单独的“垂直线”进行操作。在我之前发表了这个想法，给你点赞 :) - sehe

@sehe。我是根据原始评论工作的，“对于两种组合，路径cgk的成本将是相同的”。 - Keith

+2 这个回答会得到我的悬赏投票：我还没有找时间去做这方面的调查研究（如果我确信这些假设是成立的，我可能会更有动力；我认为一个有前途的程序员在着手编码之前应该更密切地与问题所有者合作）。 - sehe

4

编辑：添加完整解决方案

正如其他答案已经指出的那样，你的问题太难以用蛮力解决了。这种问题的起点总是模拟退火算法。我创建了一个实现该算法的小应用程序。

看待你的问题的另一种方式是最小化一个复杂的函数。此外，你还有一个额外的限制条件。我从一个随机有效的配置（符合你的限制条件）开始，然后修改那个随机解决方案，每次更改一个元素。我强制应用程序执行有效的转换。代码非常清晰。

我创建了一个模板函数，所以你只需要提供必要的函数对象和结构即可。

#include <iostream>
#include <cmath>
#include <ctime>
#include <vector>
#include <algorithm>
#include <functional>

//row [0]:  c00  c01  c02  c03   
//row [1]:  c10  c11  c12  c13  
//row [2]:  c20  c21  c22  c23 


typedef std::pair<int,int> RowColIndex;
// the deeper vector has size 3 (aei for example)
// the outer vector has size 4
typedef std::vector<std::vector<RowColIndex> > Matrix;

size_t getRandomNumber(size_t up)
{
    return rand() % up;
}

struct Configuration
{
    Configuration(const Matrix& matrix) : matrix_(matrix){}
    Matrix matrix_;
};

std::ostream& operator<<(std::ostream& os,const Configuration& toPrint)
{
    for (size_t row = 0; row < toPrint.matrix_.at(0).size(); row++)
    {
        for (size_t col = 0; col < toPrint.matrix_.size(); col++)
        {
            os << toPrint.matrix_.at(col).at(row).first  << "," 
               << toPrint.matrix_.at(col).at(row).second << '\t';
        }
        os << '\n';
    }   
    return os;
}

struct Energy 
{ 
    double operator()(const Configuration& conf)
    {
        double result = 0;
        for (size_t col = 0; col < conf.matrix_.size(); col++)
        {
            for (size_t row =0; row < conf.matrix_.at(col).size(); row++)
            {
                result += pow(static_cast<double>(row) - static_cast<double>(conf.matrix_.at(col).at(row).first),2) +
                          pow(static_cast<double>(col) - static_cast<double>(conf.matrix_.at(col).at(row).second),2);
            }
        }
        return result;
    }
};

size_t calculateNewColumn(std::vector<int>& isAlreadyUse)
{
    size_t random;
    do
    {
        random = getRandomNumber(isAlreadyUse.size());
    }
    while (isAlreadyUse.at(random) != 0);

    isAlreadyUse.at(random) = 1;
    return random;
}

Configuration createConfiguration(size_t numberOfRow,size_t numberOfColumn)
{
    //create suitable matrix
    Matrix matrix;
    //add empty column vector
    for (size_t col = 0; col < numberOfColumn; col++)
        matrix.push_back(std::vector<RowColIndex>());

    //loop over all the element
    for (size_t row = 0; row < numberOfRow; row++)
    {
        std::vector<int> isAlreadyUse(numberOfColumn);
        for (size_t col = 0; col < numberOfColumn; col++)
        {
            size_t newCol = calculateNewColumn(isAlreadyUse);
            matrix.at(newCol).push_back(std::make_pair(row,col));
        }
    }   

    return Configuration(matrix);
}


struct CreateNewConfiguration
{
    Configuration operator()(const Configuration& conf)
    {
        Configuration result(conf);

        size_t fromRow = getRandomNumber(result.matrix_.at(0).size());

        size_t fromCol = getRandomNumber(result.matrix_.size());
        size_t toCol = getRandomNumber(result.matrix_.size());

        result.matrix_.at(fromCol).at(fromRow) = conf.matrix_.at(toCol).at(fromRow);
        result.matrix_.at(toCol).at(fromRow) = conf.matrix_.at(fromCol).at(fromRow);

        return result;
    }
};

template<typename Conf,typename CalcEnergy,typename CreateRandomConf>
Conf Annealing(const Conf& start,CalcEnergy energy,CreateRandomConf createNewConfiguration,
               int maxIter = 100000,double minimumEnergy = 1.0e-005)
{
    Conf first(start);
    int iter = 0;
    while (iter < maxIter && energy(first) > minimumEnergy )
    {
        Configuration newConf(createNewConfiguration(first));
        if( energy(first) > energy(newConf))
        {
            first = newConf;
        }
        iter++;
    }
    return first;
}

int main(int argc,char* argv[])
{

    size_t nRows = 25;
    size_t nCols = 25;
    std::vector<Configuration> res;
    for (int i =0; i < 10; i++)
    {
        std::cout << "Configuration #" << i << std::endl;
        Configuration c = createConfiguration(nRows,nCols);
        res.push_back(Annealing(c,Energy(),CreateNewConfiguration()));
    }

    std::vector<Configuration>::iterator it = res.begin();


    std::vector<Configuration>::iterator lowest = it;
    while (++it != res.end())
    {
        if (Energy()(*it) < Energy()(*lowest))
            lowest = it;
    }

    std::cout << Energy()(*lowest) << std::endl;

    std::cout << std::endl;

    std::cout << *lowest << std::endl;


    std::cin.get();
    return 0;
}

当然，你不能保证这个解决方案是最好的（这是一种启发式方法）。但它是一个很好的起点。

由于您并没有提供完整的函数成本，所以我实现了自己的函数成本，这使我可以轻松检查最终结果。你只需要提供函数成本即可完成工作。

你可能会让程序更加高效，因为有很大的改进空间，但是逻辑已经存在，你可以很容易地实现你的函数。

复杂度

该算法的复杂度为E*I*C，其中 I = 迭代次数 C = 随机配置次数（避免局部最小值） E = 能量函数（或函数成本）的计算

在这种情况下，E实际上是N*M，其中N和M是初始矩阵的尺寸。

如果您对模拟退火的结果不满意，可以尝试遗传算法。

- Alessandro Teruzzi

+1，绝对要使用模拟退火或遗传算法来解决这个问题。 - TiansHUo

目前我认为这个样本既不简单也不完整。成本函数尚未确定。我不能假设它适用于SA。如果是这样，那么这是一个非常好的选择（实现可能更有效率）。然而，现在请参阅WP：问题在于真的没有任何逻辑的“邻居”生成器（能量函数可能是离散的，甚至在随机选择的转换处是奇异的）。同样，如果已知问题域中的启发式方法，则可以完全弥补这些差距。 - sehe

我尽力根据现有的信息来实现代码，虽然我知道这个实现并不高效，但我认为清晰比快速更重要。代码的目的是展示一个算法。 - Alessandro Teruzzi

3

你可以采用递归的方式解决这个问题。

该方法的输入是要计算的第一个向量的索引，而向量在函数外共享。

对于只剩下两行的情况，可以使用回溯法来计算。在这种情况下，您只需要找到成本更低的配对即可。

对于有多于两行的情况，您应该调用下一个索引的方法，获取部分结果，然后再次使用回溯法计算最小值。

当流程回到第一个向量时，您可以将结果合并为最终结果。

- robermorales

1

值得注意的是，对于一些有趣的路径成本选择，存在一个多项式时间算法，例如，如果路径成本是边缘成本的总和，则可以通过在第i行和第i+1行上运行匈牙利算法来找到最优解。

- the guy formerly known as d

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- sehe · Accepted Answer

Posting more utility code

see github: https://gist.github.com/1233012#file_new.cpp

This is basically an much better approach to generating all possible permutations based on much simpler approach (therefore I had no real reason to post it before: As it stands right now, it doesn't do anything more than the python code).

I decided to share it anyways, as you might be able to get some profit out of this as the basis for an eventual solution.

Pro:

much faster

smarter algorithm (leverages STL and maths :))

instruction optimization

storage optimization

generic problem model

model and algorithmic ideas can be used as basis for proper algorithm

basis for a good OpenMP parallelization (n-way, for n rows) designed-in (but not fleshed out)

Contra:

The code is much more efficient at the cost of flexibility: adapting the code to build in logic about the constraints and cost heuristics would be much easier with the more step-by-step Python approach

All in all I feel that my C++ code could be a big win IFF it turns out that Simulated Annealing is appropriate given the cost function(s); The approach taken in the code would give

a highly efficient storage model

a highly efficient way to generate random / closely related new grid configurations

convenient display functions

Mandatory (abritrary...) benchmark data point (comparison to the python version:)
  a  b  c  d e
  f  g  h  i j
  k  l  m  n o
  p  q  r  s t

Result: 207360000

real  0m13.016s
user  0m13.000s
sys   0m0.010s

到目前为止，我们得到了以下内容：

从描述中我了解到你有一个基本的图形，如
必须构建一条访问网格中所有节点的路径 (哈密顿回路)。
额外的约束是后续节点必须从下一个rank中取出 (a-d、e-h、i-l 是三个rank；一旦访问了最后一个rank的节点，路径必须继续与第一个rank中的任何未访问节点连接)
边缘具有权重，因为它们具有相关的成本。然而，权重函数对于图算法来说不是传统的，因为成本取决于整个路径，而不仅仅是每个边缘的端点。

根据这些信息，我认为我们处于“全覆盖”问题的领域中（需要A*算法，最著名的是Knuths Dancing Links论文）。

具体而言，如果没有更多的信息（路径的等价性、成本函数的特定属性），则获得满足约束条件的“最便宜”的哈密顿路径的最佳已知算法将是：

生成所有可能的这样的路径
计算每个实际成本函数
选择最小成本路径

因此，我已经开始编写了一个非常愚蠢的暴力生成器，可以在一个通用的NxM网格中生成所有可能的唯一路径。

宇宙的终结

对于3×4的示例网格，输出结果为4！³=13824种可能的路径...将其推广到6×48列，得到6！⁴⁸=1.4×10¹³⁷种可能性。很明显，如果没有进一步优化，这个问题是难以处理的（NP难或其他什么——我从来没有记住过那些微妙的定义）。

运行时间的爆炸声震耳欲聋：

3×4（测量）需要大约0.175秒
4×5（测量）需要大约6分5秒（在快速机器上无输出运行并在PyPy 1.6下运行）
5×6需要大约10年9个月...

在48x6的情况下，我们将会看到...什么...8.3x10¹⁰⁷ 年（请仔细阅读）

在线查看：http://ideone.com/YsVRE

无论如何，这是Python代码（所有预设为2×3网格）

#!/usr/bin/python
ROWS = 2
COLS = 3

## different cell representations
def cell(r,c): 
    ## exercise for the reader: _gues_ which of the following is the fastest
    ## ...
    ## then profile it :)
    index = COLS*(r) + c
    # return [ r,c ]
    # return ( r,c )
    # return index
    # return "(%i,%i)" % (r,c)

    def baseN(num,b,numerals="abcdefghijklmnopqrstuvwxyz"):
        return ((num == 0) and numerals[0]) or (baseN(num // b, b, numerals).lstrip(numerals[0]) + numerals[num % b])

    return baseN(index, 26)

ORIGIN = cell(0,0)

def debug(t): pass; #print t
def dump(grid): print("\n".join(map(str, grid)))

def print_path(path):
    ## Note: to 'normalize' to start at (1,1) node:
    # while ORIGIN != path[0]: path = path[1:] + path[:1] 
    print " -> ".join(map(str, path))

def bruteforce_hamiltonians(grid, whenfound):
    def inner(grid, whenfound, partial):

        cols = len(grid[-1]) # number of columns remaining in last rank
        if cols<1:
            # assert 1 == len(set([ len(r) for r in grid ])) # for debug only
            whenfound(partial)                             # disable when benchmarking
            pass
        else:
            #debug(" ------ cols: %i ------- " % cols)

            for i,rank in enumerate(grid):
                if len(rank)<cols: continue
                #debug("debug: %i, %s (partial: %s%s)" % (i,rank, "... " if len(partial)>3 else "", partial[-3:]))
                for ci,cell in enumerate(rank):
                    partial.append(cell)
                    grid[i] = rank[:ci]+rank[ci+1:] # modify grid in-place, keeps rank

                    inner(grid, whenfound, partial)

                    grid[i] = rank # restore in-place
                    partial.pop()
                break
        pass
    # start of recursion
    inner(grid, whenfound, [])

grid = [ [ cell(c,r) for r in range(COLS) ] for c in range(ROWS) ]

dump(grid)

bruteforce_hamiltonians(grid, print_path)