在 c++ 中通过网格/矩阵找到成本优化的路径

Question

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

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. 

Answer 1

发布更多实用程序代码

参见github：https://gist.github.com/1233012#file_new.cpp< /a>

这基本上是一种更好的方法，可以基于更简单的方法生成所有可能的排列（因此我之前没有真正的理由发布它：就目前而言，它不做除了 python 代码之外的任何内容）。

无论如何，我决定分享它，因为您也许可以从中获得一些利润，作为最终解决方案的基础。

优点：

• 快得多
  • 更智能的算法（利用 STL 和数学:)）
  • 指令优化
  • 存储优化
• 通用问题模型
• 模型和算法思想可以作为正确算法的基础
• 设计良好 OpenMP 并行化（n 路，n 行）的基础（但未充实）

魂斗罗：

• 代码的效率要高得多，但代价是灵活性：使用更分步的 Python 方法，调整代码以构建有关约束和成本启发式的逻辑会容易得多

总而言之，我觉得我的 C++ 代码可能是一个大胜利 IFF 事实证明，考虑到成本函数，模拟退火是合适的；代码中采用的方法将给出

• 高效的存储模型
• 一种生成随机/密切相关的新网格配置的高效方法
• 方便的显示功能

<小时>

强制（abritary...）基准数据点（与python版本比较：）

<前><代码> abcde
弗吉吉
克勒姆诺
普奎斯特

结果：207360000

真实0米13.016秒
用户 0m13.000s
系统0m0.010s

这是我们到目前为止所得到的：

• 从描述中我收集到的建议是您有一个像

• < p>必须构建一条访问网格中所有节点的路径（Hamiltonian循环）。

• 额外的约束是后续节点必须从下一个等级中获取（ad、eh、il 是三个等级；一旦访问了最后一个rank中的节点，该路径必须从第一个rank

  中的任何未访问节点继续。

• 边是加权的，因为它们具有关联的成本。但是，权重函数对于图算法来说不是传统的，因为成本取决于。完整的路径，而不仅仅是每个路径的终点边。

鉴于此，我相信我们正处于“Full Cover”问题的领域（需要 A* 算法，最著名的是 Knuths Dancing Links 论文）。

具体如果没有进一步的信息（路径的等价性、成本函数的特定属性），获得满足约束的“最便宜”哈密尔顿路径的最著名算法将是生成

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

这就是为什么我开始并编写一个非常愚蠢的强力生成器，它生成 NxM 通用网格中可能的所有唯一路径。

3×4 样本网格的宇宙尽头

输出为 4!³ = 13824 条可能的路径...推断为 6×48 列，得到 6!⁴⁸ = 1.4×10¹³⁷ 种可能性。很明显如果不进一步优化，问题将难以解决（NP Hard 之类的东西——我从来不记得那些微妙的定义）。

运行时间的爆炸性增长是震耳欲聋的：

• 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)

Answer 2

首先，一个观察结果有一点帮助。

我认为 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) 成本。我想这是可以做到的，但需要付出很多努力。相比之下，后续的通行证会便宜一些。

Answer 3

编辑：添加了完整的解决方案

正如其他答案已经指出的那样，你的问题太难了，无法用暴力来面对。这类问题的出发点始终是模拟退火。我创建了一个实现该算法的小应用程序。

查看问题的另一种方法是最小化复杂函数。此外，您对可能的解决方案还有额外的限制。我从随机有效配置（满足您的约束）开始，然后修改每次更改一个元素的随机解决方案。我强制应用程序仅执行有效的转换。代码非常清楚。

我创建了一个模板函数，因此您只需提供必要的函数对象和结构。

#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 是初始矩阵的维度。

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

Answer 4

您可以递归地解决该问题。

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

对于剩余两行的情况，可以使用回溯来计算解决方案。在这种情况下，您只需找到更便宜的配对即可。

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

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

Answer 5

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