[ create a new paste ] login | about

Link: http://codepad.org/HCeGa7st    [ raw code | fork ]

Evetro - C++, pasted on May 21:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
/**
 * @file Graph5.cpp
 * @brief matchings finding
 * @author Dmitri Kuvshinov
 */

#include "Graph4.h" // maxflow algorithms

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

/// fill capacities-matrix A for maxflow algorithm for bipartite graph adjacency matrix B of size X x Y
/**
 * 0-1 matrix is assumed.
 * Source will have index (X + Y), sink will have index (X + Y + 1), hence A must be
 * of size (X + Y + 2) x (X + Y + 2).
 */
template<class OutWeightMatrix, class InWeightMatrix>
void prepareBipartiteForMaxFlow
  (
    OutWeightMatrix &A      // матрица пропускных способностей для алгоритма максимального потока
  , const InWeightMatrix &B // матрица смежности двудольного графа
  , const unsigned X        // размер первой доли, количество строк B
  , const unsigned Y        // размер второй доли, количество столбцов B
  )
{
  typedef typename OutWeightMatrix::value_type Value;
  const unsigned source = X + Y, sink = source + 1, N = sink + 1;

  // из X можно придти только в Y
  for (unsigned i = 0; i < X; ++i)
  {
    for (unsigned j = 0; j < X; ++j)
      A(i, j) = Value(0);

    for (unsigned j = 0; j < Y; ++j)
      A(i, X + j) = B(i, j);

    A(i, source) = Value(0);
    A(i, sink) = Value(0);
  }

  // из каждого элемента Y можно придти только в сток
  for (unsigned i = X; i < source; ++i)
  {
    for (unsigned j = 0; j < sink; ++j)
      A(i, j) = Value(0);

    A(i, sink) = Value(1);
  }

  // из источника можно придти в каждый элемент X
  for (unsigned i = 0; i < X; ++i)
    A(source, i) = Value(1);
  for (unsigned i = X; i < N; ++i)
    A(source, i) = Value(0);

  // из стока никуда нельзя придти
  for (unsigned i = 0; i < N; ++i)
    A(sink, i) = Value(0);
}


/*
Вспомогательная функция, заполняющая вектор индексов 
найденным паросочетанием по матрице потока.
*/

/// extract matching as index vector from a 0-1-flow matrix of bipartite X-Y-graph
/**
 * Flow matrix is of size (X+Y+2) x (X+Y+2) as in prepareBipartiteForMaxFlow function.
 * IndexVector matching must be of size X. 
 * When there is no matching for a given index in X-part then it's given value Y
 * (while correct indices for Y-part begin with 0 and end up with Y-1.
 *
 * @return true if the matching is semiperfect.
 */
template<class IndexVector, class WeightMatrix>
bool matchingFromFlow
  (
    IndexVector &matching    // паросочетание в виде отношения X -> Y
  , const WeightMatrix &flow // матрица пропускных способностей
  , const unsigned X         // размер первой доли, длина matching
  , const unsigned Y         // размер второй доли
  )
{
  typedef typename WeightMatrix::value_type Value;
  
  unsigned matched = 0;
  for (unsigned i = 0; i < X; ++i)
  {
    matching[i] = Y;

    for (unsigned j = 0; j < Y; ++j)
    {
      if (flow(i, X + j) > Value(0))
      {
        ++matched;
        matching[i] = j;
        break;
      }
    }
  }

  return matched == X;
}


/*
Вспомогательный класс-обёртка. Нужен, чтобы указать способ построения паросочетания (MaxFlowWrapper).
По умолчанию считается, что этот способ задаёт алгоритм построения максимального потока.
Однако, можно определить специализацию по любому другому типу, что используется для
"вписывания" алгоритма Хопкрофта-Карпа в ту же схему.
*/

template<class MaxFlowWrapper>
struct FullMatching
{
/*
Функция, находящая (по возможности полное) паросочетание в двудольном графе,
заданном X x Y матрицей смежности B. Сохраняет паросочетание в вектор индексов matching.
По умолчанию использует алгоритм построения максимального потока, передаваемый с помощью
"обёрточного" класса MaxFlow в виде статической функции-члена maxFlow.
*/
  /// find (possibly perfect) matching for a bipartite graph given as X x Y adjacency matrix B
  template<class IndexVector, class BipartiteWeightMatrix>
  static bool find
    (
      IndexVector &matching          // паросочетание в виде отношения X -> Y
    , const BipartiteWeightMatrix &B // матрица смежности двудольного графа X x Y
    , const unsigned X               // размер первой доли, длина matching, количество строк в B
    , const unsigned Y               // размер второй доли, количество столбцов B
    )
  {
    // получить матрицу пропускных способностей из матрицы двудольного графа B
    const unsigned N = X + Y + 2;
    NaiveMatrix<int> A(N, N);
    prepareBipartiteForMaxFlow(A, B, X, Y);

    // рассчитать максимальный поток
    NaiveMatrix<int> F(N, N);
    MaxFlowWrapper::maxFlow(A, F, N - 2, N - 1, N);

    // получить паросочетание из потока
    return matchingFromFlow(matching, F, X, Y);
  }
};


/*
Обёртки трёх алгоритмов построения максимального потока 
для передачи этих алгоритмов функции fullMatching.
*/

/// Edmonds-Karp maximal flow algorithm (O(N5)) wrapper for fullMatching
struct EdmondsKarpMaxFlow
{
  template<class WeightMatrix>
  static typename WeightMatrix::value_type maxFlow
    (
      const WeightMatrix &A 
    , WeightMatrix &F       
    , const unsigned s      
    , const unsigned t
    , const unsigned N
    )
  {
    return maxFlowFordFulkerson<EdmondsKarp>(A, F, s, t, N);
  }
};

/// basic preflow pushing maximal flow algorithm (O(N4)) wrapper for fullMatching
struct PreflowPushingMaxFlow
{
  template<class WeightMatrix>
  static typename WeightMatrix::value_type maxFlow
    (
      const WeightMatrix &A 
    , WeightMatrix &F       
    , const unsigned s      
    , const unsigned t
    , const unsigned N
    )
  {
    return maxFlowPushing<HeightsBFS>(A, F, s, t, N);
  }
};

/// push relabel maximal flow algorithm wrapper (O(N3)) for fullMatching
struct PushRelabelMaxFlow
{
  template<class WeightMatrix>
  static typename WeightMatrix::value_type maxFlow
    (
      const WeightMatrix &A 
    , WeightMatrix &F       
    , const unsigned s      
    , const unsigned t
    , const unsigned N
    )
  {
    return maxFlowPushRelabel<HeightsBFS>(A, F, s, t, N);
  }
};


/*
Реализация алгоритма Хопкрофта-Карпа.
Класс используется для:
  а) удобства реализации в виде набора функций с разделяемыми данными (поля класса);
  б) для использования в качестве параметра MaxFlowWrapper шаблона FullMatching.
*/

class HopcroftKarp
{
  typedef std::vector<unsigned> Adjacency;    // вектор соседей вершины (список смежности)
  typedef std::vector<Adjacency> Adjacencies; // вектор списков смежности (представление графа)

  const unsigned X, Y, N; // размеры доль и индекс дополнительной вершины (N == X + Y)
  Adjacency pair; // текущее паросочетание
  Adjacencies G;  // вспомогательный граф, построенный на основе исходного двудольного графа

  typedef std::vector<unsigned> Distances; // вектор расстояний (метки фронтов)
  Distances dist; // текущие метки фронтов вершин G (расстояния от свободных вершин из X)  

  /*
  Поиск в ширину для вычисления меток фронтов.
  */
  bool bfs()
  {
    std::queue<unsigned> Q;
    for (unsigned v = 0; v < X; ++v)
    {
      if (pair[v] == N)
      {
        dist[v] = 0;
        Q.push(v);
      }
      else
        dist[v] = Max<unsigned>::value();
    }

    dist[N] = Max<unsigned>::value();
    while (!Q.empty())
    {
      const unsigned v = Q.front();
      Q.pop();

      if (v != N)
      {
        for (Adjacency::const_iterator p = G[v].begin(), pend = G[v].end(); p != pend; ++p)
        {
          const unsigned upair = pair[*p];
          if (dist[upair] == Max<unsigned>::value())
          {
            dist[upair] = dist[v] + 1;
            Q.push(upair);
          }
        }
      }
    }

    return dist[N] != Max<unsigned>::value();
  }


  /*
  Поиск в глубину, строящий и применяющий (если такая будет построена) чередующуюся цепь.
  */
  bool dfs(const unsigned v)
  {
    if (v != N)
    {
      for (Adjacency::const_iterator p = G[v].begin(), pend = G[v].end(); p != pend; ++p)
      {
        const unsigned upair = pair[*p];
        if (dist[upair] == dist[v] + 1)
        {
          if (dfs(upair))
          {
            pair[*p] = v;
            pair[v] = *p;
            return true;
          }
        }
      }

      dist[v] = Max<unsigned>::value();
      return false;
    }

    return true;
  }

public:

  /*
  Конструктор подготавливает разделяемые данные.
  */
  template<class BipartiteWeightMatrix>
  HopcroftKarp
    (
      const BipartiteWeightMatrix &B
    , const unsigned X
    , const unsigned Y
    )
      : X(X), Y(Y), N(X + Y)
      , pair(N + 1, N), dist(N + 1), G(N + 1)
  {
    Adjacency buffer;
    buffer.reserve( (std::max)(X, Y) );

    for (unsigned i = 0; i < X; ++i)
    {
      for (unsigned j = 0; j < Y; ++j)
        if (B(i, j) > 0)
        {
          const unsigned y = X + j;
          buffer.push_back(y);
          G[y].push_back(i); // prone of cache miss and array reallocation, yet simple and mixed up
        }

      G[i] = buffer;
      buffer.clear();
    }
  }

  /*
  Алгоритм Хопкрофта-Карпа.
  */
  template<class IndexVector>
  unsigned matching(IndexVector &match)
  {
    unsigned result = 0;
    while (bfs())
    {
      for (unsigned v = 0; v < X; ++v)
      {
        if (pair[v] == N && dfs(v))
          ++result;
      }
    }

    for (unsigned i = 0; i < X; ++i)
      match[i] = pair[i] - X;

    return result;
  }
};


/*
Специализация шаблона FullMatching для HopcroftKarp.
*/

template<>
struct FullMatching<HopcroftKarp>
{
  template<class IndexVector, class BipartiteWeightMatrix>
  static bool find
    (
      IndexVector &matching
    , const BipartiteWeightMatrix &B
    , const unsigned X
    , const unsigned Y
    )
  {
    HopcroftKarp HK(B, X, Y);
    return HK.matching(matching) == (std::min)(X, Y);
  }
};


/*
Функция для удобного применения класса FullMatching.
Требует явного указания способа построения паросочетания MaxFlowWrapper.
*/

template<class MaxFlowWrapper, class IndexVector, class BipartiteWeightMatrix>
bool fullMatching
  (
    IndexVector &matching
  , const BipartiteWeightMatrix &B
  , const unsigned X
  , const unsigned Y
  )
{
  return FullMatching<MaxFlowWrapper>::find(matching, B, X, Y); 
};


///////////////////////////////////////////////////////////////////////////////
// тесты
#include <iostream>
#include <iterator>

/*
Сгенерировать случайную матрицу смежности двудольного графа, в которой есть полное по X
паросочетание и до M дополнительных рёбер.
*/
void randomWithFullMatching(NaiveMatrix<int> &B, unsigned M)
{
  const unsigned X = B.size1(), Y = B.size2(), N = (std::min)(X, Y);
  // prepare a semiperfect matching
  std::vector<unsigned> pair(N);
  for (unsigned i = 0; i < N; ++i)
    pair[i] = i;
  
  std::random_shuffle(pair.begin(), pair.end());
  for (unsigned i = 0; i < N; ++i)
    B(i, pair[i]) = 1;

  // additional M edges
  M = (std::min)(M, X * Y - N);
  unsigned x = std::rand() % X;
  while (M--)
  {
    // try some random offset until match
    for (unsigned bound = 0; bound < 4u * Y; ++bound) // to prevent eternal loop
    {
      const unsigned y = (x + std::rand() - (unsigned)(RAND_MAX / 2)) % Y;
      if (B(x, y) == 0)
      {
        B(x, y) = 1;
        break;
      }
    }

    if (++x == X)
      x = 0;
  }
}


class TestCase5
{
  typedef std::vector<int> Match;
  static void print(const Match &M)
  {
    std::copy(M.begin(), M.end(), std::ostream_iterator<int>(std::cout, " "));
    std::cout << std::endl;
  }

  NaiveMatrix<int> B;

public:
  
  void prepare(const unsigned N, const unsigned M)
  {
    NaiveMatrix<int> B0(N, N);
    randomWithFullMatching(B0, M);

    if (N < 9)
    {
      std::cout << "\n";
      B.print();
      std::cout << "\n";
    }

    B = B0;
  };

  template<class MaxFlowWrapper>
  bool test1()
  {
    Match match(B.size1());
    bool res = fullMatching<MaxFlowWrapper>(match, B, B.size1(), B.size2());
    std::cout << "\n";
    print(match);
    return res;
  }
};


int main()
{
  using namespace std;
  srand(static_cast<unsigned>(time(NULL)));
  TestCase5 tc;
  do
  {
    tc.prepare(5, 10);
    cout << "Edmonds-Karp " << tc.test1<EdmondsKarpMaxFlow>();
    cout << "Push-Relabel " << tc.test1<PreflowPushingMaxFlow>();
    cout << "Relabel-to-Start " << tc.test1<PushRelabelMaxFlow>();
    cout << "Hopcroft-Karp " << tc.test1<HopcroftKarp>();
  }
  while (cin.get() != ' ');
}


Create a new paste based on this one


Comments: