Найти оптимизированный по стоимости путь через сетку/матрицу в С++

Ответ 1

Проводка кода служебной программы

см. github: https://gist.github.com/1233012#file_new.cpp

Это в основном гораздо лучший подход к созданию всех возможных перестановок, основанных на гораздо более простом подходе (поэтому у меня не было реальной причины опубликовать его раньше: поскольку он стоит прямо сейчас, он не делает ничего больше, чем код python).

Я решил поделиться им так или иначе, так как вы могли бы получить некоторую прибыль от этого в качестве основы для возможного решения.

Pro:

• намного быстрее
  • умный алгоритм (использует STL и математику:))
  • оптимизация команд
  • оптимизация хранилища
• общая модель проблемы
• модель и алгоритмические идеи могут быть использованы в качестве основы для правильного алгоритма
• для хорошей рассылки OpenMP (n-way, для n строк), разработанной (но не сглаженной)

Contra:

• Код намного эффективнее за счет гибкости: адаптация кода для построения логики об ограничениях и эвристике затрат будет намного проще с более пошаговым подходом к Python.

В целом я чувствую, что мой код на С++ может быть большим выигрышем IFF, оказывается, что Simulated Annealing подходит с учетом функции (-ов) затрат; Подход, принятый в коде, даст

• высокоэффективная модель хранения
• высокоэффективный способ генерации случайных/тесно связанных новых конфигураций сетки
• удобные функции отображения

---

Обязательная (abritrary...) контрольная точка данных (сравнение с версией python:)

  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

Вот что мы подняли до сих пор:

• Из описания я подбираю предположение, что у вас есть базовый граф вроде

• должен быть построен путь, который посещает все узлы сетки (гамильтонов цикл).

• Дополнительное ограничение состоит в том, что последующие узлы должны быть взяты из следующего ранга (ad, eh, il, являющегося тремя рангами, после того, как был посещен node из последнего ранга, путь должен продолжаться с любым невидимым node из первого ранга

• Края взвешиваются, поскольку они связаны со стоимостью. Однако функция веса не является традиционной для алгоритмов графа, поскольку стоимость зависит от полного пути, а не только от конечных точек каждого ребра.

В свете этого я считаю, что мы находимся в сфере проблем "полного покрытия" (требуя алгоритма A *, самого известного из бумаги Knuths Dancing Links).

В частности Без дополнительной информации (эквивалентность путей, конкретных свойств функции стоимости) наиболее известным алгоритмом для получения "дешевого" гамильтонова пути, который удовлетворяет ограничениям, будет

• генерировать все возможные пути
• рассчитать фактическую функцию стоимости для каждого
• выберите путь с минимальной стоимостью

Вот почему я отправился и написал действительно тупой генератор грубой силы, который генерирует все уникальные пути, возможные в общей сетке NxM.

Конец Вселенной

Выход для сетки с образцами 3 × 4 равен 4! ³= 13824 возможных путей... Экстраполируя это на 6 × 48 столбцов, приводит к 6! ⁴⁸= 1.4 × 10 ¹³⁷ возможности. Очень ясно , что без дальнейшей оптимизации проблема недостижима (NP Hard или что-то - я никогда не помню довольно тонких определений).

Взрыв времени исполнения оглушает:

• 3 × 4 (измеряется) занимает около 0,175 с
• 4 × 5 (измеренный) занял около 6 мкс (работает без вывода и под PyPy 1.6 на быстрой машине)
• 5 × 6 займет примерно 10 лет и 9 месяцев...

В 48x6 мы будем смотреть на... что... 8.3x10 ¹⁰⁷ lightyears (внимательно прочитайте)

Смотрите это в прямом эфире: 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)

Ответ 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). Я думаю, это можно сделать, но приложим немало усилий. Последующие проходы будут дешевы в сравнении.

Ответ 3

EDIT: ДОБАВЛЕННОЕ ПОЛНОЕ РЕШЕНИЕ

Как уже указывалось в других ответах, проблема слишком сложная, чтобы столкнуться с грубой силой. Отправной точкой такого рода проблем является Имитированный отжиг. Я создал небольшое приложение, реализующее алгоритм.

Еще один способ увидеть вашу проблему - свести к минимуму сложную функцию. Также у вас есть дополнительное ограничение на возможное решение. Я начинаю со случайной допустимой конфигурации (которая соответствует вашим ограничениям), затем я модифицировал это случайное решение, изменяющее элемент за раз. Я заставляю приложение выполнять только действительный переход. Код довольно ясный.

Я создал функцию шаблона, поэтому вам нужно просто предоставить необходимые функциональные объекты и структуру.

#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 * я * C, где I = количество итераций C = количество случайных конфигураций (во избежание локального минимума) E = расчет энергетической функции (или стоимости функции)

В этом случае E фактически является N * M, где N и M - размеры исходной матрицы.

Если вас не устраивает результат имитации отжига, вы можете попробовать генетические алгоритмы.

Ответ 4

Рекурсивно вы можете решить проблему.

Ввод метода - это индекс первого вектора для вычисления, а вектор - вне функции.

В случае оставшихся двух строк решение может быть вычислено с использованием обратного отслеживания. В этом случае вам нужно найти менее дорогостоящую пару.

В этом случае существует более двух строк, вы должны вызвать метод со следующим индексом, получить частичный результат, а затем вычислить минимальное значение, используя обратный трафик снова.

Когда поток возвращается к первому вектору, вы можете консолидировать результат в качестве конечного результата.

Ответ 5

Стоит отметить, что для некоторых интересных вариантов стоимости пути существует многопользовательский алгоритм, например, если стоимость пути является суммой затрат на грань, то оптимальное решение можно найти, выполнив для всех i, венгерский алгоритм на строках я и я + 1.