1 //
2 // graph2.cc
3 //
4
5 #include "wx/string.h"
6
7 #include <list>
8 #include <map>
9 #include <queue>
10 #include <set>
11 #include <vector>
12 #include "edge.h"
13 #include "graph.h"
14 #include "main.h"
15 #include "math.h"
16 #include "matrix.h"
17 #include "polynomial.h"
18 #include "vertex.h"
19
20
is_bridge(Edge * e) const21 bool Graph::is_bridge (Edge *e) const
22 {
23 Graph *tmp = new Graph (*this);
24 bool res;
25
26 tmp->remove (tmp->find (tmp->find (e->v->label),
27 tmp->find (e->w->label), e->directed));
28 res = tmp->is_connected ();
29
30 delete tmp;
31 return !res;
32 }
33
is_undirected() const34 bool Graph::is_undirected () const
35 {
36 Graph::e_const_iterator eit;
37
38 for (eit = e_begin (); eit != e_end (); ++eit)
39 if ((*eit)->directed)
40 return false;
41 return true;
42 }
43
is_connected(bool dir,Vertex * start) const44 bool Graph::is_connected (bool dir, Vertex *start) const
45 {
46 bool done;
47 unsigned int marks;
48 Graph::e_const_iterator eit;
49 Graph::v_const_iterator vit;
50
51 if (V.empty ())
52 return true; // null graph is connected
53 if (num_edges () < (order () - 1))
54 return false; // not enough edges
55
56 // Unmark all edges and vertices initially
57 for (eit = e_begin (); eit != e_end (); ++eit)
58 (*eit)->mark = 0;
59 for (vit = v_begin (); vit != v_end (); ++vit)
60 (*vit)->mark = 0;
61
62 if (!start)
63 start = V[0];
64 start->mark = 1;
65
66 // Algorithm: Iterate through edge list, marking vertices adjacent
67 // to already-marked vertices, until no more can be marked
68 marks = 1;
69 do {
70 done = true;
71 for (eit = e_begin (); eit != e_end (); ++eit) {
72 Edge *e = *eit;
73
74 if (e->mark)
75 continue;
76 if (!e->v->mark && !e->w->mark)
77 continue;
78
79 if (!dir || !e->directed) {
80 if (e->v->mark != e->w->mark)
81 marks++;
82 e->v->mark = e->w->mark = e->mark = 1;
83 done = false;
84 } else {
85 if (!e->v->mark && e->w->mark)
86 continue;
87 marks++;
88 e->w->mark = e->mark = 1;
89 done = false;
90 }
91 }
92 } while (!done);
93
94 if (marks < order ())
95 return false;
96 return true;
97 }
98
99 // TODO: this is horribly slow (but at least correct!)
is_strongly_connected() const100 bool Graph::is_strongly_connected () const
101 {
102 Graph::v_const_iterator vit;
103
104 for (vit = v_begin (); vit != v_end (); ++vit) {
105 if (!is_connected (true, *vit))
106 return false;
107 }
108 return true;
109 }
110
eulericity(bool & euler,bool & semi,wxString & tour) const111 void Graph::eulericity (bool &euler, bool &semi, wxString &tour) const
112 {
113 int cnt;
114 Graph::v_const_iterator vit;
115 const Vertex *v, *w;
116 bool undir;
117
118 euler = false;
119 semi = false;
120
121 if (!is_connected ())
122 return;
123
124 undir = is_undirected ();
125 if (undir) {
126 cnt = 0;
127 v = w = 0;
128 for (vit = v_begin (); vit != v_end (); ++vit) {
129 int odd = (*vit)->degree () & 1;
130 if (!odd)
131 continue;
132 cnt += 1;
133 if (cnt == 1)
134 v = *vit;
135 else if (cnt == 2)
136 w = *vit;
137 if (cnt > 2)
138 break;
139 }
140 if (cnt == 0) // Eulerian
141 euler = true;
142 else if (cnt == 2) // Semi-Eulerian
143 semi = true;
144 } else {
145 cnt = 0;
146 v = w = 0;
147 for (vit = v_begin (); vit != v_end (); ++vit) {
148 int ind = (*vit)->indegree ();
149 int outd = (*vit)->outdegree ();
150 if (ind != outd) {
151 ++cnt;
152 if (outd == ind + 1) {
153 if (v)
154 break;
155 v = *vit;
156 } else if (ind == outd + 1) {
157 if (w)
158 break;
159 w = *vit;
160 } else
161 break;
162 }
163 }
164 if (cnt == 0)
165 euler = true;
166 else if ((cnt == 2) && v && w)
167 semi = true;
168 }
169
170 // Use Fleury's Algorithm
171 Graph *red = new Graph (*this);
172 std::vector<const Edge *> tour_e;
173 std::vector<const Vertex *> tour_v;
174 const Vertex *curr_this, *curr_red;
175
176 // We always start with vertex v when semi-Eulerian; otherwise pick any
177 curr_this = semi ? v : V[0];
178 curr_red = red->find (curr_this->label);
179 tour_v.push_back (curr_this);
180
181 while (1) {
182 if ((curr_red->degree () == 0) && (curr_red->outdegree () == 0))
183 break;
184 // pick an edge
185 Edge *sel, *fail = 0;
186 Graph::e_const_iterator eit, last = curr_red->e_end ();
187 for (eit = curr_red->e_begin (); eit != last; ++eit) {
188 if ((*eit)->directed && ((*eit)->v != curr_red))
189 continue;
190 if (red->is_bridge (*eit)) {
191 fail = *eit;
192 continue;
193 }
194 if (!(undir || (*eit)->directed)) {
195 fail = *eit;
196 continue;
197 }
198 break;
199 }
200 sel = (eit == last) ? fail : *eit;
201 Vertex *dst_red = curr_red->opposite (sel);
202 Edge *e = find (curr_this, find (dst_red->label), sel->directed);
203 curr_red = dst_red;
204 red->remove (sel);
205 curr_this = curr_this->opposite (e);
206 tour_e.push_back (e);
207 tour_v.push_back (curr_this);
208 }
209
210 delete red;
211
212 std::vector<const Vertex *>::const_iterator cvit;
213 int num = 0;
214 for (cvit = tour_v.begin (); cvit != tour_v.end (); ++cvit) {
215 if (num > 0)
216 tour += wxT(", ");
217 ++num;
218 tour += (*cvit)->label;
219 }
220 }
221
mark_shortest_path(Vertex * v1,Vertex * v2)222 void Graph::mark_shortest_path (Vertex *v1, Vertex *v2)
223 {
224 Graph::e_const_iterator eit;
225 Graph::v_iterator vit;
226
227 // Method: Dijkstra's Algorithm
228 // (vertex->mark is shortest distance from source + 1)
229 // (if vertex->mark is negative, it's an estimate)
230
231 // Mark all vertices as belonging to the V\S set
232 for (vit = v_begin (); vit != v_end (); ++vit)
233 (*vit)->mark = 0;
234 v1->mark = 1;
235
236 while (1) {
237 // If every vertex is marked, break
238 for (vit = v_begin (); vit != v_end (); ++vit)
239 if ((*vit)->mark < 1)
240 break;
241 if (vit == v_end ())
242 break; // all marked
243
244 // Revise estimates
245 for (vit = v_begin (); vit != v_end (); ++vit) {
246 Vertex *v = *vit;
247
248 if (v->mark > 0)
249 continue; // already marked
250 v->mark = 0;
251 eit = v->e_begin ();
252 for (; eit != v->e_end (); ++eit) {
253 Edge *e = *eit;
254 Vertex *alt = v->opposite (e);
255 if (e->directed && (v != e->w))
256 continue;
257
258 if (alt->mark < 1)
259 continue;
260 if (((alt->mark + e->weight) < -v->mark) ||
261 (v->mark == 0))
262 v->mark = -(alt->mark + e->weight);
263 }
264 }
265
266 // Find smallest estimate, and mark it permanently
267 Vertex *min = v1;
268 bool none_yet = true;
269 int min_val = 0;
270 for (vit = v_begin (); vit != v_end (); ++vit) {
271 if ((*vit == v1) || ((*vit)->mark >= 0))
272 continue;
273 if ((min_val > -(*vit)->mark) || none_yet) {
274 min = *vit;
275 min_val = -min->mark;
276 none_yet = false;
277 }
278 }
279 if (min == v1)
280 break; // filled all possible
281
282 // Mark this vertex permanently
283 min->mark = -min->mark;
284
285 if (min == v2)
286 break; // cut search short
287 }
288
289 // Adjust marks to be correct distances
290 for (vit = v_begin (); vit != v_end (); ++vit) {
291 Vertex *v = *vit;
292
293 if (v->mark > 0)
294 v->mark -= 1;
295 else
296 v->mark = -1;
297 }
298 }
299
300 // TODO: digraph fixes?
diameter(bool sel)301 int Graph::diameter (bool sel)
302 {
303 int best_dist;
304 Vertex *best_src, *best_dst;
305 Graph::v_iterator vit1, vit2;
306
307 if (order () < 2)
308 return 0;
309
310 best_dist = -1;
311 best_src = best_dst = 0;
312 for (vit1 = v_begin (); (vit1 + 1) != v_end (); ++vit1) {
313 mark_shortest_path (*vit1, 0);
314 for (vit2 = vit1 + 1; vit2 != v_end (); ++vit2) {
315 if (best_dist < (*vit2)->mark) {
316 best_src = *vit1;
317 best_dst = *vit2;
318 best_dist = best_dst->mark;
319 }
320 }
321 }
322
323 if (best_dist < 0)
324 return 0;
325
326 if (!sel)
327 return best_dist;
328
329 // Select longest shortest path marked
330 mark_shortest_path (best_src, best_dst);
331 unselect_all ();
332 select (best_dst);
333 while (best_dst != best_src) {
334 Graph::e_const_iterator eit;
335 Vertex *alt = 0;
336
337 eit = best_dst->e_begin ();
338 for (; eit != best_dst->e_end (); ++eit) {
339 alt = best_dst->opposite (*eit);
340 if (best_dst->mark == (alt->mark + (*eit)->weight))
341 break;
342 }
343 if (eit == best_dst->e_end ())
344 break; // Uh, oh!
345 select (*eit);
346 select (alt);
347 best_dst = alt;
348 }
349
350 return best_dist;
351 }
352
353 // TODO: digraph fixes?
radius(bool sel)354 int Graph::radius (bool sel)
355 {
356 int best_dist, sub_best_dist;
357 Vertex *best_src;
358 Graph::v_iterator vit1, vit2;
359
360 if (order () < 2)
361 return 0;
362
363 best_dist = -1;
364 best_src = 0;
365 for (vit1 = v_begin (); vit1 != v_end (); ++vit1) {
366 mark_shortest_path (*vit1, 0);
367 sub_best_dist = -1;
368 for (vit2 = v_begin (); vit2 != v_end (); ++vit2) {
369 if (sub_best_dist < (*vit2)->mark)
370 sub_best_dist = (*vit2)->mark;
371 }
372 if (sub_best_dist == -1)
373 continue;
374 if ((best_dist == -1) || (sub_best_dist < best_dist)) {
375 best_src = *vit1;
376 best_dist = sub_best_dist;
377 }
378 }
379
380 if (best_dist < 0)
381 return 0;
382
383 if (sel) {
384 // Select most central vertex
385 unselect_all ();
386 select (best_src);
387 }
388
389 return best_dist;
390 }
391
adjacency_matrix() const392 Matrix Graph::adjacency_matrix () const
393 {
394 Matrix ret (order (), order ());
395 unsigned int i, j;
396
397 for (j = 0; j < order (); ++j)
398 for (i = 0; i < order (); ++i) {
399 if (i == j) {
400 ret (i, j) = 0;
401 continue;
402 }
403 Edge *e = find (V[i], V[j], true);
404 if (!e) {
405 e = find (V[i], V[j], false);
406 if (e && e->directed && (e->v == V[j]))
407 e = NULL;
408 }
409 ret (i, j) = e ? e->weight : 0;
410 }
411
412 return ret;
413 }
414
bfs(Vertex * v,wxString & s)415 void Graph::bfs (Vertex *v, wxString &s)
416 {
417 Graph::e_const_iterator eit;
418 Graph::v_iterator vit;
419 std::queue<Vertex *> vq;
420 int cnt = 1;
421
422 // Set all marks to 0
423 for (vit = v_begin (); vit != v_end (); ++vit)
424 (*vit)->mark = 0;
425
426 traversal_visit (v, s, cnt);
427 vq.push (v);
428
429 while (!vq.empty ()) {
430 v = vq.front ();
431 vq.pop ();
432
433 for (eit = v->e_begin (); eit != v->e_end (); ++eit) {
434 if ((*eit)->directed && ((*eit)->v != v))
435 continue;
436 Vertex *opp = v->opposite (*eit);
437 if (opp->mark)
438 continue;
439 traversal_visit (opp, s, cnt);
440 vq.push (opp);
441 }
442 }
443 }
444
dfs(Vertex * v,wxString & s)445 void Graph::dfs (Vertex *v, wxString &s)
446 {
447 Graph::v_iterator vit;
448
449 // Set all marks to 0
450 for (vit = v_begin (); vit != v_end (); ++vit)
451 (*vit)->mark = 0;
452
453 dfs_do (v, s, 1);
454 }
455
lt_weight(const Edge * a,const Edge * b)456 static bool lt_weight (const Edge *a, const Edge *b)
457 {
458 return a->weight < b->weight;
459 }
460
461 struct lt_vertex
462 {
operator ()lt_vertex463 bool operator () (const Vertex *v1, const Vertex *v2) const
464 {
465 return ((v1->x < v2->x) ||
466 ((v1->x == v2->x) && (v1->y < v2->y)));
467 }
468 };
469
minimum_spanning_tree(std::set<Edge * > & result) const470 void Graph::minimum_spanning_tree (std::set<Edge *> &result) const
471 {
472 // Use Kruskal's Algorithm
473
474 result.clear ();
475 // set of sets, each set contains a connected component
476 std::set< std::set<Vertex *> > comps;
477 // list to contain the edges sorted by their weight
478 std::list<Edge *> w_edge;
479
480 Graph::v_const_iterator vit;
481 Graph::e_const_iterator eit;
482
483 // firstly, the vertices are in separated components
484 for (vit = v_begin (); vit != v_end (); ++vit) {
485 std::set<Vertex *> tmp;
486 tmp.insert (*vit);
487 comps.insert (tmp);
488 }
489
490 // add all the edges to the edge set
491 for (eit = e_begin (); eit != e_end (); ++eit)
492 w_edge.push_back (*eit);
493 w_edge.sort (lt_weight);
494
495 // loop while not all edges are processed
496 while (!w_edge.empty ()) {
497 // get edge with the smallest weight
498 Edge *current = w_edge.front ();
499 // does it straddle two components?
500 std::set< std::set<Vertex *> >::iterator sit = comps.begin (),
501 foundv = comps.end (), foundw = comps.end ();
502 for (; sit != comps.end (); ++sit) {
503 if ((*sit).find (current->v) != (*sit).end ())
504 foundv = sit;
505 if ((*sit).find (current->w) != (*sit).end ())
506 foundw = sit;
507 }
508
509 if (foundv != foundw) {
510 // in different components => add to MST
511 std::set<Vertex *> un;
512 result.insert (current);
513 un.insert (foundv->begin (), foundv->end ());
514 un.insert (foundw->begin (), foundw->end ());
515 comps.erase (foundv);
516 comps.erase (foundw);
517 comps.insert (un);
518 }
519
520 // done with this edge
521 w_edge.pop_front ();
522 }
523 }
524
525 //#define FF_DEBUG
526 // This implementation is based on what was found on p170 of:
527 // "Graphs, Algorithms, and Optimization" (Kocay and Kreher)
528 // It is an augmenting-path implementation of Ford-Fulkerson,
529 // using a BFS approach to finding shortest augmenting paths.
ford_fulkerson(Vertex * src,Vertex * dest)530 int Graph::ford_fulkerson (Vertex *src, Vertex *dest)
531 {
532 Graph *g;
533 std::map<Vertex *, Vertex *> PrevPt;
534 std::map<Vertex *, int> ResCap;
535
536 std::map<Vertex *, Vertex *>::iterator vvit;
537 Graph::v_iterator vit;
538 Graph::e_iterator eit;
539 int total_flow = 0;
540
541 // Handle undirected edges by making a new graph, replacing all the
542 // undirected edges with a pair of directed edges (one each way), of
543 // the same capacity
544 g = new Graph (*this);
545 src = g->find (src->label);
546 dest = g->find (dest->label);
547 std::list<Edge *> undir_e; // list of undirected edges
548 for (eit = g->e_begin (); eit != g->e_end (); ++eit) {
549 if (!(*eit)->directed)
550 undir_e.push_back (*eit);
551 }
552 std::list<Edge *>::iterator u_eit;
553 for (u_eit = undir_e.begin (); u_eit != undir_e.end (); ++u_eit) {
554 Edge *e = *u_eit;
555 Vertex *v = e->v, *w = e->w;
556 int weight = e->weight;
557 g->remove (e);
558 g->add (new Edge (v, w, true, weight));
559 g->add (new Edge (w, v, true, weight));
560 }
561
562 // Usage of vertex->mark: [bitmask]
563 // 1 = is in ScanQ (for speeding up lookup)
564
565 // Set all flow to 0
566 for (eit = g->e_begin (); eit != g->e_end (); ++eit) {
567 Edge *e = *eit;
568 e->flow = 0;
569 }
570 for (vit = g->v_begin (); vit != g->v_end (); ++vit) {
571 Vertex *v = *vit;
572 PrevPt[v] = 0;
573 v->mark = 0;
574 //ResCap[v] = 0;
575 }
576
577 // Compute residual capacity at source vertex
578 int max_cap = 0;
579 for (eit = src->e_begin (); eit != src->e_end (); ++eit)
580 max_cap += (*eit)->weight;
581
582
583 while (true) { // search for an augmenting path
584 std::queue<Vertex *> ScanQ;
585
586 for (vit = g->v_begin (); vit != g->v_end (); ++vit) {
587 (*vit)->mark &= ~1; // ScanQ is empty
588 ResCap[*vit] = 0; // XXX: check this
589 }
590 ResCap[src] = max_cap;
591 ScanQ.push (src);
592 src->mark |= 1;
593 #ifdef FF_DEBUG
594 std::cerr << "FF: Looking for augmenting path...\n";
595 #endif
596 do {
597 Vertex *u = ScanQ.front ();
598 ScanQ.pop ();
599 for (eit = u->e_begin (); eit != u->e_end (); ++eit) {
600 Edge *e = *eit;
601 Vertex *v = u->opposite (e);
602 if (v->mark & 1) // v is in ScanQ
603 continue;
604 // v is not in ScanQ
605
606 #ifdef FF_DEBUG
607 std::cerr << "\tTrying (" << u->label << "," <<
608 v->label << "), " << v->label <<
609 "->mark=" << v->mark << ", " <<
610 " e->flow=" << e->flow << "\n";
611 #endif
612
613 // Assume directed edge
614
615 if (e->v == u) { // forward edge
616 if (e->weight > e->flow) {
617 ScanQ.push (v);
618 v->mark |= 1;
619 PrevPt[v] = u;
620
621 int del = e->weight - e->flow;
622 if (ResCap[u] < del)
623 ResCap[v] = ResCap[u];
624 else
625 ResCap[v] = del;
626 if (v == dest)
627 goto ford_fulkerson_augment;
628 }
629 } else { // backward edge
630 if (e->flow > 0) {
631 ScanQ.push (v);
632 v->mark |= 1;
633 PrevPt[v] = u;
634
635 if (ResCap[u] < e->flow)
636 ResCap[v] = ResCap[u];
637 else
638 ResCap[v] = e->flow;
639 if (v == dest)
640 goto ford_fulkerson_augment;
641 }
642 }
643 }
644
645 } while (!ScanQ.empty ());
646
647 ford_fulkerson_all_done:
648 // All done!
649 for (eit = e_begin (); eit != e_end (); ++eit) {
650 Edge *e_orig = *eit;
651 Vertex *v_dup = g->find (e_orig->v->label),
652 *w_dup = g->find (e_orig->w->label);
653 if (e_orig->directed) {
654 Edge *e_dup = g->find (v_dup, w_dup, true);
655 e_orig->flow = e_dup->flow;
656 } else {
657 Edge *e1_dup = g->find (v_dup, w_dup, true),
658 *e2_dup = g->find (w_dup, v_dup, true);
659 int f1 = e1_dup->flow, f2 = e2_dup->flow;
660 if (f1 > f2)
661 e_orig->flow = f1 - f2;
662 else
663 e_orig->flow = f2 - f1;
664 }
665 }
666 return total_flow;
667
668 ford_fulkerson_augment:
669 // AUGMENTFLOW(t):
670 {
671 Vertex *v = dest;
672 Vertex *u = PrevPt[v];
673 int delta = ResCap[dest];
674 #ifdef FF_DEBUG
675 std::cerr << "\tAugmentFlow w/ delta=" << delta << "\n";
676
677 std::cerr << "\t" << v->label;
678 #endif
679 if (delta == 0)
680 goto ford_fulkerson_all_done;
681 while (u) {
682 #ifdef FF_DEBUG
683 std::cerr << " <- " << u->label;
684 #endif
685 Edge *e = g->find (u, v, true);
686 if (e) // forward edge
687 e->flow += delta;
688 else { // backward edge
689 e = g->find (v, u, true);
690 e->flow -= delta;
691 }
692 v = u;
693 u = PrevPt[v];
694 }
695 total_flow += delta;
696 #ifdef FF_DEBUG
697 std::cerr << ".\n";
698 #endif
699 }
700 for (vit = v_begin (); vit != v_end (); ++vit)
701 PrevPt[*vit] = 0;
702 }
703 }
704
chromatic_number() const705 int Graph::chromatic_number () const
706 {
707 switch (order ()) {
708 case 0: return -1; // undefined!
709 case 1: return 1;
710 }
711
712 if (num_edges () == 0)
713 return 1;
714
715 // try colouring with only a few colours
716 {
717 // TODO: not an ideal approach...
718 Graph g (*this);
719 if (g.try_colouring (2))
720 return 2;
721 if (g.try_colouring (3))
722 return 3;
723 }
724 // Note that the only value that try_colouring is guaranteed to succeed
725 // for (when a colouring is possible) is 2. Thus we can only go up to
726 // 3, since try_colouring(3) failing does not preclude the graph from
727 // being 3-colourable. Once we've retrofitted a better algorithm into
728 // try_colouring(), then we may safely increase this number.
729
730 Polynomial p = chromatic_polynomial ();
731 unsigned int i;
732
733 for (i = 1; i < order (); ++i)
734 if (p.eval (i) != 0)
735 return i;
736 return i;
737 }
738
traversal_visit(Vertex * v,wxString & s,int & cnt)739 void Graph::traversal_visit (Vertex *v, wxString &s, int &cnt)
740 {
741 v->mark = cnt++;
742 if (cnt > 2)
743 s += wxT(", ");
744 s += v->label;
745 }
746
dfs_do(Vertex * v,wxString & s,int cnt)747 int Graph::dfs_do (Vertex *v, wxString &s, int cnt)
748 {
749 Graph::e_const_iterator eit;
750
751 traversal_visit (v, s, cnt);
752
753 for (eit = v->e_begin (); eit != v->e_end (); ++eit) {
754 if ((*eit)->directed && ((*eit)->v != v))
755 continue;
756 Vertex *opp = v->opposite (*eit);
757 if (!opp->mark)
758 cnt = dfs_do (opp, s, cnt);
759 }
760
761 return cnt;
762 }
763
chromatic_polynomial() const764 Polynomial Graph::chromatic_polynomial () const
765 {
766 Graph *g = flattened ();
767 Polynomial p;
768
769 p = g->chromatic_poly (0.0, 1.0);
770 delete g;
771 return p;
772 }
773
chromatic_poly(double startP,double endP)774 Polynomial Graph::chromatic_poly (double startP, double endP)
775 {
776 Graph::e_const_iterator eit;
777 Graph::v_const_iterator vit;
778 unsigned int edge_num, vertex_num;
779 bool subtract;
780
781 edge_num = E.size ();
782 vertex_num = order ();
783
784 // Some quick, simple cases
785 if (vertex_num == 0)
786 return Polynomial (1);
787 else if (vertex_num == 1)
788 return Polynomial (1, 0); // x
789 else if (vertex_num == 2) {
790 if (edge_num == 0)
791 return Polynomial (1, 0, 0); // x^2
792 else
793 return Polynomial (1, -1, 0); // x^2 - x
794 } else if (vertex_num == 3) {
795 switch (edge_num) {
796 case 0: return Polynomial (1, 0, 0, 0);
797 case 1: return Polynomial (1, -1, 0, 0);
798 case 2: return Polynomial (1, -2, 1, 0);
799 case 3: return Polynomial (1, -3, 2, 0);
800 }
801 }
802 if (edge_num == 0) {
803 Polynomial p;
804 p[vertex_num] = 1;
805 return p;
806 }
807
808 // Remove degree 0 vertices
809 std::queue<Vertex *> vq;
810 int cnt = 0;
811 for (vit = v_begin (); vit != v_end (); ++vit)
812 if ((*vit)->degree () == 0)
813 vq.push (*vit);
814 while (!vq.empty ()) {
815 remove (vq.front ());
816 vq.pop ();
817 ++cnt;
818 }
819 if (cnt) {
820 Polynomial p;
821 p[cnt] = 1;
822 return chromatic_poly (startP, endP) * p;
823 }
824
825 // Remove degree 1 vertices
826 int tot = 0;
827 do {
828 std::queue<Vertex *> vq;
829 cnt = 0;
830 for (vit = v_begin (); vit != v_end (); ++vit)
831 if ((*vit)->degree () == 1)
832 vq.push (*vit);
833 while (!vq.empty ()) {
834 if (vq.front ()->degree () == 1) {
835 remove (vq.front ());
836 ++cnt;
837 }
838 vq.pop ();
839 }
840 tot += cnt;
841 } while (cnt > 0);
842 if (tot)
843 return chromatic_poly (startP, endP) // XXX: fine tune
844 * Polynomial::binomial (-1, tot);
845
846 // Assumptions: |V| >= 4, |E| >= 1
847
848 // Firstly, if the graph is not connected, calculate polynomials
849 // for each component, then compute their product
850 if (!is_connected ()) {
851 Graph gm; // marked subgraph
852 std::queue<Vertex *> vq; // marked vertices
853 std::queue<Edge *> eq; // marked edges
854
855 // Extract marked subgraph
856 for (vit = v_begin (); vit != v_end (); ++vit) {
857 if (!(*vit)->mark)
858 continue;
859 vq.push (*vit);
860 gm.add (new Vertex ((*vit)->label, 0, 0));
861 }
862 for (eit = e_begin (); eit != e_end (); ++eit) {
863 if (!(*eit)->mark)
864 continue;
865 eq.push (*eit);
866 gm.add (new Edge (gm.find ((*eit)->v->label),
867 gm.find ((*eit)->w->label)));
868 }
869
870 while (!eq.empty ()) {
871 remove (eq.front ());
872 eq.pop ();
873 }
874 while (!vq.empty ()) {
875 remove (vq.front ());
876 vq.pop ();
877 }
878
879 // XXX: fine-tune
880 double weight = (double) order () / (order () + gm.order ());
881 double midP = startP + weight * (endP - startP);
882
883 Polynomial p1 = chromatic_poly (startP, midP);
884 Polynomial p2 = gm.chromatic_poly (midP, endP);
885
886 return p1 * p2;
887 }
888
889 // Assumptions: |V| >= 4, |E| >= 1, connected
890
891 // Is it a tree?
892 if (edge_num == (vertex_num - 1)) {
893 // binomial expansion: x(x-1)^(v-1)
894 return Polynomial (1, 0) *
895 Polynomial::binomial (-1, edge_num);
896 }
897
898 // Assumptions: |E| >= |V| >= 4, connected, cyclic
899
900 // Special cases for common graphs
901 if ((vertex_num == 4) && (edge_num == 4)) {
902 // either C4, or irregular
903 if ((V[0]->degree () != 2) ||
904 (V[1]->degree () != 2) ||
905 (V[2]->degree () != 2))
906 return Polynomial (1, -4, 5, -2, 0); // irregular
907 else
908 return Polynomial (1, -4, 6, -3, 0); // C4
909 }
910 if ((vertex_num == 4) && (edge_num == 5))
911 return Polynomial (1, -5, 8, -4, 0); // K4 - 1 edge
912 if ((vertex_num == 4) && (edge_num == 6))
913 return Polynomial (1, -6, 11, -6, 0); // K4
914 if ((vertex_num == 5) && (edge_num == 10))
915 return Polynomial (1, -10, 35, -50, 24, 0); // K5
916
917 // Assumptions: |E| >= |V| >= 5, connected, cyclic
918
919 // Handle complete graphs
920 if ((2 * edge_num) == (vertex_num * (vertex_num - 1))) {
921 int i, sign = (vertex_num & 1) ? -1 : 1;
922 Polynomial p;
923
924 p[vertex_num] = 1;
925 for (i = 0; i < (signed) vertex_num; ++i, sign = -sign)
926 p[i] = sign * Math::stirling (vertex_num, i);
927
928 return p;
929 }
930
931 // Assumptions: |E| >= |V| >= 5, connected, cyclic, not complete
932
933 // Determine which reduction formula to use. For relatively sparse
934 // graphs (4*e < (v^2 + v)), use the subtractive, otherwise use
935 // the additive. This probably needs some fine tuning.
936 subtract = ((4 * edge_num) < (vertex_num * (vertex_num + 1)));
937
938 setProgress (startP);
939 double midP = (startP + endP) / 2;
940
941 if (subtract) {
942 // P(G,x) = P(G-e,x) - P(G|e,x)
943
944 Graph gm (*this);
945
946 Edge *gme = 0;
947 // Find a lightly connected edge
948 for (eit = gm.e_begin (); eit != gm.e_end (); ++eit) {
949 int wt;
950 wt = (*eit)->v->degree () + (*eit)->w->degree ();
951 if (!gme || (wt < gme->mark)) {
952 gme = *eit;
953 gme->mark = wt;
954 }
955 }
956
957 remove (find (find (gme->v->label), find (gme->w->label)));
958 gm.identify (gme->v, gme->w);
959
960 //return (chromatic_poly () - gm.chromatic_poly ());
961 Polynomial a = chromatic_poly (startP, midP);
962 setProgress (midP);
963 Polynomial b = gm.chromatic_poly (midP, endP);
964 setProgress (endP);
965 return a - b;
966 } else {
967 // P(G,x) = P(G+e,x) + P(G|e,x)
968
969 Graph gm (*this);
970
971 // Find two non-complete vertices to add an edge to
972 Vertex *v1 = 0, *v2 = 0;
973 for (vit = v_begin (); vit != v_end (); ++vit) {
974 if ((*vit)->degree () >= (vertex_num - 1))
975 continue;
976 if (!v1)
977 v1 = *vit;
978 else if (!are_adjacent (v1, *vit)) {
979 v2 = *vit;
980 break;
981 }
982 }
983
984 add (new Edge (v1, v2));
985 gm.identify (gm.find (v1->label), gm.find (v2->label));
986
987 //return (chromatic_poly () + gm.chromatic_poly ());
988 Polynomial a = chromatic_poly (startP, midP);
989 setProgress (midP);
990 Polynomial b = gm.chromatic_poly (midP, endP);
991 setProgress (endP);
992 return a + b;
993 }
994
995 // Never reached
996 }
997
try_colouring(unsigned int colours)998 bool Graph::try_colouring (unsigned int colours)
999 {
1000 Graph::v_iterator vit;
1001 Graph::e_iterator eit;
1002
1003 if (colours < 2)
1004 return false;
1005
1006 for (vit = v_begin (); vit != v_end (); ++vit)
1007 (*vit)->mark = -1;
1008
1009 // Special case: 2 colours
1010 if (colours == 2) {
1011 bool done = false, next_component = true;
1012 unsigned int marks_made = 0;
1013 while (!done) {
1014 done = true;
1015 for (vit = v_begin (); vit != v_end (); ++vit) {
1016 Vertex *v = *vit;
1017 if (v->mark >= 0)
1018 continue;
1019 if (next_component) {
1020 next_component = false;
1021 v->mark = 0;
1022 ++marks_made;
1023 continue;
1024 }
1025 for (eit = v->e_begin (); eit != v->e_end ();
1026 ++eit) {
1027 Vertex *v2 = v->opposite (*eit);
1028 if (v2->mark < 0)
1029 continue;
1030 v->mark = 1 - v2->mark;
1031 ++marks_made;
1032 done = false;
1033 break;
1034 }
1035 }
1036 if (done && (marks_made < order ())) {
1037 done = false;
1038 next_component = true;
1039 }
1040 }
1041 return check_colouring (colours);
1042 }
1043
1044 // Try the greedy algorithm, with some small modifications:
1045 // - the next vertex to colour will be the one with the least
1046 // number of uncoloured neighbours
1047 // TODO:
1048 // * don't assume a connected graph!
1049 bool done = false, next_component = true, failed = false;
1050 bool *can_use = new bool[colours];
1051 while (!done) {
1052 // Find vertex with the least number of uncoloured neighbours,
1053 // but with still 1 coloured neighbour
1054 Vertex *best_v = 0;
1055 unsigned int uncol_neigh = order ();
1056 for (vit = v_begin (); vit != v_end (); ++vit) {
1057 Vertex *v = *vit;
1058 unsigned int num_col = 0, num_uncol = 0;
1059 if (v->mark >= 0)
1060 continue;
1061 if (next_component) {
1062 best_v = v;
1063 next_component = false;
1064 break;
1065 }
1066 for (eit = v->e_begin (); eit != v->e_end (); ++eit) {
1067 Vertex *valt = v->opposite (*eit);
1068 if (valt->mark < 0)
1069 ++num_uncol;
1070 else
1071 ++num_col;
1072 }
1073 if (num_col < 1)
1074 continue;
1075 if (num_uncol < uncol_neigh) {
1076 best_v = v;
1077 uncol_neigh = num_uncol;
1078 }
1079 }
1080 if (!best_v) {
1081 if (next_component) {
1082 done = true;
1083 break; // no more components to colour
1084 }
1085 next_component = true;
1086 continue;
1087 }
1088 //std::cout << "Colouring " << best_v->label << "\n";
1089 unsigned int i;
1090 for (i = 0; i < colours; ++i)
1091 can_use[i] = true;
1092 for (eit = best_v->e_begin (); eit != best_v->e_end (); ++eit) {
1093 Vertex *valt = best_v->opposite (*eit);
1094 if (valt->mark >= 0)
1095 can_use[valt->mark] = false;
1096 }
1097 for (i = 0; i < colours; ++i)
1098 if (can_use[i]) {
1099 best_v->mark = i;
1100 break;
1101 }
1102 if (best_v->mark < 0) {
1103 // Uh, oh! Failed!
1104 failed = true;
1105 break;
1106 }
1107 }
1108 delete can_use;
1109
1110 if (failed)
1111 return false;
1112 return check_colouring (colours);
1113 }
1114
check_colouring(unsigned int colours) const1115 bool Graph::check_colouring (unsigned int colours) const
1116 {
1117 Graph::e_const_iterator eit;
1118
1119 for (eit = e_begin (); eit != e_end (); ++eit) {
1120 if ((*eit)->v->mark == (*eit)->w->mark)
1121 return false;
1122 }
1123 return true;
1124 }
1125