/**
* @file Graph4.cpp
* @brief finding maximal flow in a single source/sink network
* @author Dmitri Kuvshinov
*/
#include <limits>
#include <vector>
#include <list>
#include <queue>
#include <algorithm>
#include "Graph2.h"
// вспомогательная структура для выбора максимальных и минимальных значений
// некоторого числового типа (T)
template<class T>
struct MaxMinUtil
{
template<bool UseInf> // использовать значение "бесконечность"
struct Internal // <true> variant
{
static T maximum()
{
return std::numeric_limits<T>::infinity();
}
static T minimum()
{
return -std::numeric_limits<T>::infinity();
}
};
template<>
struct Internal<false> // не использовать значение "бесконечность"
{
static T maximum()
{
return (std::numeric_limits<T>::max)();
}
static T minimum()
{
return (std::numeric_limits<T>::min)();
}
};
};
/// provides maximal value of scalar type T
template<class T>
struct Max
{
static T value()
{ // использовать "бесконечность", если тип удовлетворяет IEC-559 == IEEE-754
return MaxMinUtil<T>::Internal<std::numeric_limits<T>::is_iec559>::maximum();
}
};
/// provides minimal value of scalar type T
template<class T>
struct Min
{
static T value()
{ // использовать "бесконечность", если тип удовлетворяет IEC-559 == IEEE-754
return MaxMinUtil<T>::Internal<std::numeric_limits<T>::is_iec559>::minimum();
}
};
/// simplifies augment function choice
struct EdmondsKarp
{
/*
Вспомогательная функция "помечивания", используемая реализацией алгоритма Форда-Фалкерсона.
Ищет дополняющую (s, t)-цепь. Если таковая существует (h[t] < бесконечности), то
она может быть восстановлена по массиву предшественников Pr.
*/
/// builds an augmenting (s, t)-chain
template<class WeightMatrix, class WeightVector, class ChoiceVector, class Previous>
static void augment
(
/* вход */
const WeightMatrix &A // матрица пропускных способностей дуг
, const WeightMatrix &F // матрица потоков
, unsigned s // индекс вершины-источника
, unsigned t // индекс вершины-стока
, unsigned N // количество вершин в графе
/* выход */
, WeightVector &h // метки (прирост/уменьшение потока при добавлении этой цепи)
, ChoiceVector &choice // признак "прямая" дуга (true) или "обратная" дуга (false)
, Previous &Pr // предшественники в дополняющей цепи
)
{
// очистим метки (бесконечность == признак "вершина не помечена")
for (unsigned v = 0; v < N; ++v)
h[v] = Max<typename WeightMatrix::value_type>::value(); //std::numeric_limits<double>::infinity();
// начинаем поиск в ширину с вершины s
std::queue<unsigned> Q;
Q.push(s);
Pr[s] = s; // условность
// пока не пометили вершину t и ещё есть непросмотренные вершины
while (h[t] == Max<typename WeightMatrix::value_type>::value() //std::numeric_limits<double>::infinity()
&& !Q.empty())
{
// берём следующую вершину из очереди
const unsigned w = Q.front();
Q.pop();
// проверяем дуги, выходящие из w
for (unsigned v = 0; v < N; ++v)
{
// если вершина v ещё не помечена
if (h[v] == Max<typename WeightMatrix::value_type>::value()) //std::numeric_limits<double>::infinity())
{
// и есть ещё куда увеличить поток через дугу w-v (на delta)
const typename WeightMatrix::value_type delta = A(w, v) - F(w, v);
if (delta > 0)
{
// добавим в цепь с соответствующей меткой
h[v] = std::min(h[w], delta);
Pr[v] = w;
Q.push(v);
choice[v] = true; // прямая дуга
}
}
}
// проверяем дуги, входящие в w
for (unsigned v = 0; v < N; ++v)
{
// если вершина не помечена и есть ненулевой поток из v в w
if (v != s && h[v] == Max<typename WeightMatrix::value_type>::value() //std::numeric_limits<double>::infinity()
&& F(v, w) > 0)
{
// добавим в цепь с соответствующей меткой
h[v] = (std::min)(h[w], F(v, w));
Pr[w] = v;
Q.push(v);
choice[v] = false; // обратная дуга
}
}
}
}
};
/*
Построение максимального потока через сеть алгоритмом Форда-Фалкерсона.
Тип Augmentor должен предоставлять статическую функцию augment, строящую дополняющую цепь.
Тип WeightMatrix должен предоставлять operator()(unsigned, unsigned) для доступа к элементам матрицы.
*/
/// maximal flow constuction by Ford-Fulkerson algorithm
/**
* @param A input matrix of edge capacities
* @param F output matrix of flows
* @param s source vertex index
* @param t sink vertex index
* @param N amount of vertices
* @return maximal flow from s to t through the network A
*/
template<class Augmentor, class WeightMatrix>
typename WeightMatrix::value_type maxFlowFordFulkerson
(
const WeightMatrix &A // матрица пропускных способностей рёбер (вход)
, WeightMatrix &F // матрица потоков (выход)
, unsigned s // индекс вершины-источника
, unsigned t // индекс вершины-стока
, unsigned N // количество вершин в графе
)
{
// сначала потоков нет
for (unsigned i = 0; i < N; ++i)
for (unsigned j = 0; j < N; ++j)
F(i, j) = 0;
// результат функции -- суммарный максимальный поток через цепь
typename WeightMatrix::value_type maxFlow = 0;
// вспомогательные массивы для augment
std::vector<typename WeightMatrix::value_type> h(N);
std::vector<bool> choice(N);
std::vector<unsigned> Pr(N);
for (;;)
{
// искать дополняющую (s, t)-цепь
Augmentor::augment(A, F, s, t, N, h, choice, Pr);
if (h[t] == Max<typename WeightMatrix::value_type>::value()) //std::numeric_limits<double>::infinity())
return maxFlow; // вернуть текущее значение потока, если цепь не найдена
// цепь найдена, увеличить поток, "раскручивая" цепь с конца по предшественникам
maxFlow += h[t];
for (unsigned v = t; v != s; )
{
const unsigned w = Pr[v];
if (choice[v]) // прямая дуга, увеличить поток
F(w, v) += h[t];
else // обратная дуга, уменьшить поток
F(v, w) -= h[t];
v = w;
}
}
}
/// calculate initial heights of vertices using reverse bfs from the sink
struct HeightsBFS
{
/*
Функция осуществляет заполнение вектора h высотами вершин, вычисляемыми
как расстояние от этой вершины до стока t.
Расстояние вычисляется поиском в глубину против направления рёбер, начиная
со стока.
*/
template<class WeightMatrix, class IntVector>
static void heights(const WeightMatrix &A, IntVector &h, unsigned t, unsigned N)
{
h[t] = 0;
std::queue<unsigned> Q;
std::vector<bool> visited(N, false);
Q.push(t);
visited[t] = true;
while (!Q.empty())
{
const unsigned v = Q.front();
const unsigned nh = h[v] + 1;
Q.pop();
for (unsigned w = 0; w < N; ++w)
{
if (!visited[w] && A(w, v) > 0)
{
h[w] = nh;
Q.push(w);
visited[w] = true;
}
}
}
}
};
/*
Очередь беззнаковых целых без дублирования элементов.
Для предотвращения дублирования сохраняется признак "уже поставлено в очередь"
в векторе inQ.
*/
/// a queue with unique unsigned values inside
class UniqueQueue
{
std::queue<unsigned> Q; // собственно очередь
std::vector<bool> inQ; // признак "поставлено в очередь"
public:
// конструктору нужно знать максимально возможное значение,
// чтобы подготовить вектор inQ
/// provide maximal possible value
explicit UniqueQueue(unsigned N)
: inQ(N, false)
{
}
// проверка на пустоту
bool empty() const
{
return Q.empty();
}
// получить доступ к началу очереди
unsigned front() const
{
return Q.front();
}
// убрать элемент из начала очереди
void pop()
{
inQ[Q.front()] = false;
Q.pop();
}
// добавить элемент к очереди, el < N, переданного в конструктор
void push(unsigned el)
{
if (!inQ[el])
{
inQ[el] = true;
Q.push(el);
}
}
};
/*
Базовый вариант алгоритма проталкивания предпотока.
Ставит активные вершины в очередь типа UniqueQueue.
*/
/// naive preflow pushing algorithm implementation
template<class Heights, class WeightMatrix>
static typename WeightMatrix::value_type maxFlowPushing
(
const WeightMatrix &A // матрица пропускных способностей рёбер (вход)
, WeightMatrix &F // матрица потоков (выход)
, unsigned s // индекс вершины-источника
, unsigned t // индекс вершины-стока
, unsigned N // количество вершин в графе
)
{
// вспомогательные определения типов
typedef typename WeightMatrix::value_type ValueType; // тип элементов матрицы
typedef std::vector<ValueType> WeightVector; // избытки
typedef std::vector<unsigned> HeightVector; // высоты
// максимально заполнить рёбра, исходящие из истока
for (unsigned i = 0; i < N; ++i)
F(i, s) = -(F(s, i) = A(s, i));
// инициализировать высоты вершин
HeightVector H(N);
Heights::heights(A, H, t, N);
H[s] = N;
UniqueQueue Q(N); // очередь активных вершин
WeightVector E(N, 0); // избытки потока в вершинах
for (unsigned i = 0; i < N; ++i)
{
E[i] = F(s, i); // получить поток из заполненных рёбер
if (E[i] > 0 && i != s && i != t)
Q.push(i); // если избыток > 0 => вершина активна, поместим её в очередь
}
// пока есть активные вершины
while (!Q.empty())
{
// вынем следующую вершину из очереди
unsigned i = Q.front();
Q.pop();
bool needLift = true; // признак "нужно поднять вершину"
// попытаемся распределить избыток по соседям
for (unsigned j = 0; j < N; ++j)
{
const ValueType delta = A(i, j) - F(i, j); // остаточная пропускная способность
if (delta > 0 && (H[i] == H[j] + 1)) // можно протолкнуть поток из i в j
{
needLift = false; // не надо поднимать
// протолкнуть поток, сколько получится
ValueType d = (std::min)(E[i], delta);
F(j, i) = -(F(i, j) += d);
E[i] -= d;
E[j] += d;
// если вершина, куда протолкнули поток, получила избыток,
// поместить её в очередь активных вершин
if (E[j] > 0 && j != s && j != t)
Q.push(j);
// если избытка уже нет, то можно прекратить цикл
if (E[i] == 0)
break;
}
}
// надо поднимать?
if (needLift)
{ // подъём вершины i
// ищем минимальную высоту соседей вершины i по остаточному графу
unsigned d = Max<unsigned>::value();
for (unsigned j = 0; j < N; ++j)
if (d > H[j] && A(i, j) > F(i, j))
d = H[j];
// если нашли, устанавливаем новую высоту на 1 больше минимальной
if (d < Max<unsigned>::value())
H[i] = d + 1;
}
// если вершина всё ещё имеет избыток, помещаем её в очередь активных
if (E[i] > 0 && i != s && i != t)
Q.push(i);
}
// посчитаем значение максимального потока по рёбрам, исходящим из истока
ValueType flow = 0;
for (unsigned i = 0; i < N; ++i)
if (A(s, i) > 0)
flow += F(s, i);
return flow;
}
/*
Вариант алгоритма проталкивания предпотока с выбором наивысшей вершины.
Перебирает вершины по списку, если вершина была поднята, переносит её в начало списка.
*/
/// "push-relabel to start" algorithm implementation
template<class Heights, class WeightMatrix>
static typename WeightMatrix::value_type maxFlowPushRelabel
(
const WeightMatrix &A // матрица пропускных способностей рёбер (вход)
, WeightMatrix &F // матрица потоков (выход)
, unsigned s // индекс вершины-источника
, unsigned t // индекс вершины-стока
, unsigned N // количество вершин в графе
)
{
// вспомогательные определения типов
typedef typename WeightMatrix::value_type ValueType; // тип элементов матрицы
typedef std::vector<ValueType> WeightVector; // избытки
typedef std::vector<unsigned> HeightVector; // высоты
// максимально заполнить рёбра, исходящие из истока
for (unsigned i = 0; i < N; ++i)
F(i, s) = -(F(s, i) = A(s, i));
// инициализировать высоты вершин
HeightVector H(N, 0);
Heights::heights(A, H, t, N);
H[s] = N;
std::list<unsigned> Q; // список обрабатываемых вершин
WeightVector E(N, 0); // избытки потока в вершинах
for (unsigned i = 0; i < N; ++i)
{
E[i] = F(s, i); // получить поток из заполненных рёбер
if (i != s && i != t)
Q.push_back(i); // поместить все вершины, кроме s и t, в список на обработку
}
// идти по списку, пока он не кончится
for (std::list<unsigned>::iterator q = Q.begin(); q != Q.end(); ++q)
{
// разрядка вершины с индексом *q
bool lifted = false; // была ли поднята в процессе разрядки
unsigned i = *q, j = 0;
while (E[i] > 0) // пока есть избыток, расталкивать поток по соседям
{
// проверили всех соседей?
if (j == N)
{
// поднять вершину
lifted = true;
// ищем минимальную высоту соседей вершины i по остаточному графу
unsigned d = Max<unsigned>::value() - 1;
for (unsigned j = 0; j < N; ++j)
if (d > H[j] && A(i, j) > F(i, j))
d = H[j];
H[i] = d + 1; // устанавливаем высоту на 1 больше минимальной среди соседей
j = 0; // снова будем обходить соседей
}
const ValueType delta = A(i, j) - F(i, j); // остаточная пропускная способность
if (delta > 0 && (H[i] == H[j] + 1)) // можно протолкнуть поток из i в j
{
// протолкнуть поток, сколько получится
ValueType d = (std::min)(E[i], delta);
F(j, i) = -(F(i, j) += d);
E[i] -= d;
E[j] += d;
}
++j;
}
// переставить вершину в начало списка, если она была поднята
if (lifted)
{
const unsigned el = *q;
Q.erase(q);
Q.push_front(el);
q = Q.begin();
// следующей вершиной будет та, которая стояла в начале списка
// до перестановки поднятой вершины
}
}
// посчитаем значение максимального потока по рёбрам, исходящим из истока
ValueType flow = 0;
for (unsigned i = 0; i < N; ++i)
if (A(s, i) > 0)
flow += F(s, i);
return flow;
}
///////////////////////////////////////////////////////////////////////////////
#include <iostream>
#include "Graph2.h"
template<class WeightMatrix>
void randomLayers(WeightMatrix &A, unsigned s, unsigned t, unsigned L, unsigned k, unsigned N)
{
// занулить
for (unsigned i = 0; i < N; ++i)
for (unsigned j = 0; j < N; ++j)
A(i, j) = 0;
// заполнять "по уровням": связи внутри уровня и между уровнями
const unsigned Lsz = N / L;
for (unsigned Ls = 0; Ls < N; Ls += Lsz)
{
const unsigned Lse = (std::min)(Ls + Lsz, N);
for (unsigned i = Ls; i < Lse; ++i)
{
// связи внутри уровня
for (unsigned round = (std::min)(k / 2, Lse - Ls); round > 0; )
{
const unsigned j = Ls + std::rand() % (Lse - Ls);
if (A(i, j) == 0)
{
--round;
A(i, j) = 1 + std::rand();
}
}
// связи со следующим уровнем
const unsigned L2se = (std::min)(Lse + Lsz, N);
for (unsigned j = (std::min)(k - k / 2, L2se - Lse); j > 0; --j)
A(i, Lse + std::rand() % (L2se - Lse)) += std::rand();
}
}
}
// сделать случайную сеть
template<class WeightMatrix>
void randomCapacities(WeightMatrix &A, unsigned s, unsigned t, unsigned k, unsigned N)
{
randomLayers(A, s, t, N / k, k, N); //randomMatrix(A, k, N);
for (unsigned i = 0; i < N; ++i)
for (unsigned j = i + 1; j < N; ++j)
{
const typename WeightMatrix::value_type Aij = A(i, j), Aji = A(j, i);
if (Aij > 0 && Aji > 0)
{
A(i, j) = Aij + Aji;
A(j, i) = 0;
}
}
for (unsigned i = 0; i < N; ++i)
{
A(t, i) = 0;
A(i, s) = 0;
}
A(s, t) = 1;
}
// тесты
template<class FlowValue>
class TestCase4
{
static const unsigned N = 1000; // количество вершин
NaiveMatrix<FlowValue> capacities; // тестовая матрица пропускных способностей
public:
TestCase4()
: capacities(N, N, 0)
{
}
// тест с простым графом на совпадение результата алгоритма с известным решением
bool testSimple() const
{
NaiveMatrix<> A(4, 4, 0.0), F(4, 4);
A(0, 1) = 20;
A(0, 2) = 10;
A(2, 1) = 2;
A(1, 3) = 60;
A(2, 3) = 5;
double f = maxFlowFordFulkerson<EdmondsKarp>(A, F, 0, 3, 4);
return f == 27
&& F(0, 1) == 20
&& F(0, 2) == 7
&& F(2, 1) == 2
&& F(1, 3) == 22
&& F(2, 3) == 5;
}
// тест с большим случайным графом на время работы двух алгоритмов
void testTime()
{
// vertex 0 -- source, vertex N-1 -- sink
randomCapacities(capacities, 0, N - 1, N / 4, N);
double ff, pushing;
{
NaiveMatrix<FlowValue> F(N, N);
Timer t;
ff = maxFlowFordFulkerson<EdmondsKarp>(capacities, F, 0, N - 1, N);
std::cout << ff;
std::cout << "\ntime = " << t() << std::endl;
}
{
NaiveMatrix<FlowValue> F(N, N);
Timer t;
pushing = maxFlowPushing<HeightsBFS>(capacities, F, 0, N - 1, N);
std::cout << pushing;
std::cout << "\ntime = " << t() << std::endl;
}
{
NaiveMatrix<FlowValue> F(N, N);
Timer t;
pushing = maxFlowPushRelabel<HeightsBFS>(capacities, F, 0, N - 1, N);
std::cout << pushing;
std::cout << "\ntime = " << t() << std::endl;
}
std::cout << (ff - pushing) << std::endl;
}
};
int main()
{
using namespace std;
srand(time(NULL));
cout << "for double flow-type:\n";
{
TestCase4<double> tc;
cout << tc.testSimple() << endl;
tc.testTime();
}
cout << "\nfor int flow-type:\n";
{
TestCase4<int> tc;
tc.testTime();
}
cin.get();
}