SHOGUN  6.1.3
GraphCut.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2014 Jiaolong Xu
8  * Copyright (C) 2014 Jiaolong Xu
9  */
10 
12 #include <shogun/io/SGIO.h>
13 
14 using namespace shogun;
15 
17  : CMAPInferImpl()
18 {
19  SG_UNSTABLE("CGraphCut::CGraphCut()", "\n");
20 
21  init();
22 }
23 
25  : CMAPInferImpl(fg)
26 {
27  ASSERT(m_fg != NULL);
28 
29  init();
30 }
31 
32 CGraphCut::CGraphCut(int32_t num_nodes, int32_t num_edges)
33  : CMAPInferImpl()
34 {
35  init();
36 
37  m_num_nodes = num_nodes;
38  // build s-t graph
39  build_st_graph(m_num_nodes, num_edges);
40 }
41 
43 {
44  if (m_nodes!=NULL)
45  SG_FREE(m_nodes);
46 
47  if (m_edges!=NULL)
48  SG_FREE(m_edges);
49 }
50 
51 void CGraphCut::init()
52 {
53  m_nodes = NULL;
54  m_edges = NULL;
55  m_edges_last = NULL;
56  m_num_nodes = 0;
57  m_num_edges = 0;
58 
59  m_active_first[0] = NULL;
60  m_active_last[0] = NULL;
61  m_active_first[1] = NULL;
62  m_active_last[1] = NULL;
63  m_orphan_first = NULL;
64  m_orphan_last = NULL;
65 
66  m_timestamp = 0;
67  m_flow = 0;
68  m_map_energy = 0;
69 
70  if (m_fg == NULL)
71  return;
72 
74 
76 
77  for (int32_t i = 0; i < cards.size(); i++)
78  {
79  if (cards[i] != 2)
80  {
81  SG_ERROR("This implementation of the graph cut optimizer supports only binary variables.");
82  }
83  }
84 
85  m_num_factors_at_order = SGVector<int32_t> (4);
86  m_num_factors_at_order.zero();
87 
88  for (int32_t i = 0; i < facs->get_num_elements(); i++)
89  {
90  CFactor* fac = dynamic_cast<CFactor*>(facs->get_element(i));
91 
92  int32_t num_vars = fac->get_num_vars();
93 
94  SG_UNREF(fac);
95 
96  if (num_vars > 3)
97  {
98  SG_ERROR("This implementation of the graph cut optimizer supports only factors of order <= 3.");
99  }
100 
101  ++m_num_factors_at_order[num_vars];
102 
103  }
104 
105  m_num_variables = m_fg->get_num_vars();
106  int32_t max_num_edges = m_num_factors_at_order[2] + 3 * m_num_factors_at_order[3];
107  m_num_nodes = m_num_variables + m_num_factors_at_order[3];
108 
109  // build s-t graph
110  build_st_graph(m_num_nodes, max_num_edges);
111 
112  for (int32_t j = 0; j < m_fg->get_num_factors(); j++)
113  {
114  CFactor* fac = dynamic_cast<CFactor*>(facs->get_element(j));
115  add_factor(fac);
116  SG_UNREF(fac);
117  }
118 
119  SG_UNREF(facs);
120 }
121 
122 void CGraphCut::build_st_graph(int32_t num_nodes, int32_t num_edges)
123 {
124  m_num_nodes = num_nodes;
125 
126  // allocate s-t graph
127  m_nodes = SG_MALLOC(GCNode, m_num_nodes);
128  m_edges = SG_MALLOC(GCEdge, 2 * num_edges);
129  m_edges_last = m_edges;
130 
131  for (int32_t i = 0; i < m_num_nodes; i++)
132  {
133  m_nodes[i].id = i;
134  m_nodes[i].tree_cap = 0;
135  m_nodes[i].first = NULL;
136  }
137 
138  m_num_edges = 0; // m_num_edges will be counted in add_edge()
139  m_flow = 0;
140 
141  m_active_first[0] = NULL;
142  m_active_last[0] = NULL;
143  m_active_first[1] = NULL;
144  m_active_last[1] = NULL;
145  m_orphan_first = NULL;
146  m_orphan_last = NULL;
147 
148  m_timestamp = 0;
149 }
150 
152 {
153  GCNode* node_i;
154 
155  m_active_first[0] = NULL;
156  m_active_last[0] = NULL;
157  m_active_first[1] = NULL;
158  m_active_last[1] = NULL;
159  m_orphan_first = NULL;
160  m_orphan_last = NULL;
161 
162  m_timestamp = 0;
163 
164  for (int32_t i = 0; i < m_num_nodes; i++)
165  {
166  node_i = m_nodes + i;
167  node_i->next = NULL;
168  node_i->timestamp = m_timestamp;
169 
170  if (node_i->tree_cap > 0)
171  {
172  // i is connected to the source
173  node_i->type_tree = SOURCE;
174  node_i->parent = TERMINAL_EDGE;
175  set_active(node_i);
176  node_i->dist_terminal = 1;
177  }
178  else if (node_i->tree_cap < 0)
179  {
180  // i is connected to the sink
181  node_i->type_tree = SINK;
182  node_i->parent = TERMINAL_EDGE;
183  set_active(node_i);
184  node_i->dist_terminal = 1;
185  }
186  else
187  {
188  node_i->parent = NULL;
189  }
190  }
191 }
192 
194 {
195  REQUIRE(assignment.size() == m_fg->get_cardinalities().size(),
196  "%s::inference(): the output assignment should be prepared as"
197  "the same size as variables!\n", get_name());
198 
199  // compute max flow
200  init_maxflow();
201  compute_maxflow();
202 
203  for (int32_t vi = 0; vi < assignment.size(); vi++)
204  {
205  assignment[vi] = get_assignment(vi) == SOURCE ? 0 : 1;
206  }
207 
208  m_map_energy = m_fg->evaluate_energy(assignment);
209  SG_DEBUG("fg.evaluate_energy(assignment) = %f\n", m_fg->evaluate_energy(assignment));
210  SG_DEBUG("minimized energy = %f\n", m_map_energy);
211 
212  return m_map_energy;
213 }
214 
215 void CGraphCut::add_factor(CFactor* factor)
216 {
217  SGVector<int32_t> fcards = factor->get_cardinalities();
218 
219  for (int32_t i = 0; i < fcards.size(); i++)
220  {
221  ASSERT(fcards[i] == 2);
222  }
223 
224  int32_t f_order = factor->get_num_vars();
225 
226  switch (f_order)
227  {
228  case 0:
229  break;
230  case 1:
231  {
232  SGVector<int32_t> fvars = factor->get_variables();
233  SGVector<float64_t> fenrgs = factor->get_energies();
234  ASSERT(fenrgs.size() == 2);
235  int32_t var = fvars[0];
236  float64_t v0 = fenrgs[0];
237  float64_t v1 = fenrgs[1];
238 
239  if (v0 < v1)
240  {
241  add_tweights(var, v1 - v0, 0);
242  }
243  else
244  {
245  add_tweights(var, 0, v0 - v1);
246  }
247  }
248  break;
249  case 2:
250  {
251  SGVector<int32_t> fvars = factor->get_variables();
252  SGVector<float64_t> fenrgs = factor->get_energies();
253  int32_t var0 = fvars[0];
254  int32_t var1 = fvars[1];
255  float64_t A = fenrgs[0]; //E{0,0} = {y_var0, y_var1}
256  float64_t B = fenrgs[2]; //E{0,1}
257  float64_t C = fenrgs[1]; //E{1,0}
258  float64_t D = fenrgs[3]; //E{1,1}
259 
260  // Added "truncation" code below to ensure regularity / submodularity
261  if (A + D > C + B)
262  {
263  SG_DEBUG("Truncation is applied to ensure regularity / submodularity.");
264 
265  float64_t delta = A + D - C - B;
266  float64_t subtrA = delta / 3;
267  A = A - subtrA;
268  C = C + subtrA;
269  B = B + (delta - subtrA * 2) + 0.0001; // for numeric issue
270  }
271 
272  // first variabe
273  if (C > A)
274  {
275  add_tweights(var0, C - A, 0);
276  }
277  else
278  {
279  add_tweights(var0, 0, A - C);
280  }
281  // second varibale
282  if (D > C)
283  {
284  add_tweights(var1, D - C, 0);
285  }
286  else
287  {
288  add_tweights(var1, 0, C - D);
289  }
290 
291  // submodular term
292  float64_t term = B + C - A - D;
293 
294  // term >= 0 is the regularity condition.
295  // It is the sufficient and necessary condition for any function to be graph-representable
296  if (term < 0)
297  {
298  SG_ERROR("\nRegularity condition is not satisfied\n");
299  }
300 
301  add_edge(var0, var1, term, 0);
302  }
303  break;
304  case 3:
305  {
306  SGVector<int32_t> fvars = factor->get_variables();
307  SGVector<float64_t> fenrgs = factor->get_energies();
308  int32_t var0 = fvars[0];
309  int32_t var1 = fvars[1];
310  int32_t var2 = fvars[2];
311  float64_t A = fenrgs[0]; //{0,0,0}
312  float64_t E = fenrgs[1]; //{1,0,0}
313  float64_t C = fenrgs[2]; //{0,1,0}
314  float64_t G = fenrgs[3]; //{1,1,0}
315  float64_t B = fenrgs[4]; //{0,0,1}
316  float64_t F = fenrgs[5]; //{1,0,1}
317  float64_t D = fenrgs[6]; //{0,1,1}
318  float64_t H = fenrgs[7]; //{1,1,1}
319 
320  int32_t id = get_tripleId(fvars);
321  float64_t P = (A + D + F + G) - (B + C + E + H);
322 
323  if (P >= 0.0)
324  {
325  if (F - B >= 0)
326  {
327  add_tweights(var0, F - B, 0);
328  }
329  else
330  {
331  add_tweights(var0, 0, B - F);
332  }
333 
334  if (G - E >= 0)
335  {
336  add_tweights(var1, G - E, 0);
337  }
338  else
339  {
340  add_tweights(var1, 0, E - G);
341  }
342 
343  if (D - C >= 0)
344  {
345  add_tweights(var2, D - C, 0);
346  }
347  else
348  {
349  add_tweights(var2, 0, C - D);
350  }
351 
352  add_edge(var1, var2, B + C - A - D, 0);
353  add_edge(var2, var0, B + E - A - F, 0);
354  add_edge(var0, var1, C + E - A - G, 0);
355 
356  add_edge(var0, id, P, 0);
357  add_edge(var1, id, P, 0);
358  add_edge(var2, id, P, 0);
359  add_edge(id, 1, P, 0);
360  }
361  else
362  {
363  if (C - G >= 0)
364  {
365  add_tweights(var0, 0, C - G);
366  }
367  else
368  {
369  add_tweights(var0, G - C, 0);
370  }
371 
372  if (B - D >= 0)
373  {
374  add_tweights(var1, 0, B - D);
375  }
376  else
377  {
378  add_tweights(var1, D - B, 0);
379  }
380 
381  if (E - F >= 0)
382  {
383  add_tweights(var2, 0, E - F);
384  }
385  else
386  {
387  add_tweights(var2, F - E, 0);
388  }
389 
390  add_edge(var2, var1, F + G - E - H, 0);
391  add_edge(var0, var2, D + G - C - H, 0);
392  add_edge(var1, var0, D + F - B - H, 0);
393 
394  add_edge(id, var0, -P, 0);
395  add_edge(id, var1, -P, 0);
396  add_edge(id, var2, -P, 0);
397  add_tweights(id, -P, 0);
398  }
399  }
400  break;
401  default:
402  SG_ERROR("This implementation of the graph cut optimizer does not support factors of order > 3.");
403  break;
404  }
405 }
406 
407 int32_t CGraphCut::get_tripleId(SGVector<int32_t> triple)
408 {
409  // search for triple in list
410  int32_t counter = m_num_variables;
411 
412  for (int32_t i = 0; i < m_triple_list.get_num_elements(); i++)
413  {
414  SGVector<int32_t> vec = m_triple_list[i];
415 
416  if (triple[0] == vec[0] && triple[1] == vec[1] && triple[2] == vec[2])
417  {
418  return counter;
419  }
420 
421  m_num_variables++;
422  }
423  // add triple to list
424  m_triple_list.push_back(triple);
425 
426  ASSERT(counter - m_num_variables < m_num_factors_at_order[3]);
427 
428  return counter;
429 }
430 
431 void CGraphCut::add_tweights(int32_t i, float64_t cap_source, float64_t cap_sink)
432 {
433  ASSERT(i >= 0 && i < m_num_nodes);
434 
435  float64_t delta = m_nodes[i].tree_cap;
436 
437  if (delta > 0)
438  {
439  cap_source += delta;
440  }
441  else
442  {
443  cap_sink -= delta;
444  }
445 
446  m_flow += (cap_source < cap_sink) ? cap_source : cap_sink;
447 
448  m_nodes[i].tree_cap = cap_source - cap_sink;
449 }
450 
451 void CGraphCut::add_edge(int32_t i, int32_t j, float64_t capacity, float64_t reverse_capacity)
452 {
453  ASSERT(i >= 0 && i < m_num_nodes);
454  ASSERT(j >= 0 && j < m_num_nodes);
455  ASSERT(i != j);
456  ASSERT(capacity >= 0);
457  ASSERT(reverse_capacity >= 0);
458 
459  GCEdge* e = m_edges_last++;
460  e->id = m_num_edges++;
461  GCEdge* e_rev = m_edges_last++;
462  e_rev->id = m_num_edges++;
463 
464  GCNode* node_i = m_nodes + i;
465  GCNode* node_j = m_nodes + j;
466 
467  e->reverse = e_rev;
468  e_rev->reverse = e;
469  e->next = node_i->first;
470  node_i->first = e;
471  e_rev->next = node_j->first;
472  node_j->first = e_rev;
473  e->head = node_j;
474  e_rev->head = node_i;
475  e->residual_capacity = capacity;
476  e_rev->residual_capacity = reverse_capacity;
477 }
478 
479 void CGraphCut::set_active(GCNode* node_i)
480 {
481  if (node_i->next == NULL)
482  {
483  // it's not in the list yet
484  if (m_active_last[1])
485  {
486  m_active_last[1]->next = node_i;
487  }
488  else
489  {
490  m_active_first[1] = node_i;
491  }
492 
493  m_active_last[1] = node_i;
494  node_i->next = node_i;
495  }
496 }
497 
498 GCNode* CGraphCut::next_active()
499 {
500  // Returns the next active node. If it is connected to the sink,
501  // it stays in the list, otherwise it is removed from the list.
502  GCNode* node_i;
503 
504  while (true)
505  {
506  if ((node_i = m_active_first[0]) == NULL)
507  {
508  m_active_first[0] = node_i = m_active_first[1];
509  m_active_last[0] = m_active_last[1];
510  m_active_first[1] = NULL;
511  m_active_last[1] = NULL;
512 
513  if (node_i == NULL)
514  {
515  return NULL;
516  }
517  }
518 
519  // remove it from the active list
520  if (node_i->next == node_i)
521  {
522  m_active_first[0] = NULL;
523  m_active_last[0] = NULL;
524  }
525  else
526  {
527  m_active_first[0] = node_i->next;
528  }
529 
530  node_i->next = NULL;
531 
532  // a node in the list is active iff it has a parent
533  if (node_i->parent != NULL)
534  {
535  return node_i;
536  }
537  }
538 }
539 
541 {
542  GCNode* current_node = NULL;
543  bool active_set_found = true;
544 
545  // start the main loop
546  while (true)
547  {
548  if (sg_io->get_loglevel() == MSG_DEBUG)
549  test_consistency(current_node);
550 
551  GCEdge* connecting_edge;
552 
553  // find a path from source to sink
554  active_set_found = grow(connecting_edge, current_node);
555 
556  if (!active_set_found)
557  {
558  break;
559  }
560 
561  if (connecting_edge == NULL)
562  {
563  continue;
564  }
565 
566  m_timestamp++;
567 
568  // augment that path
569  augment_path(connecting_edge);
570 
571  // adopt orphans, rebuild the search tree structure
572  adopt();
573  }
574 
575  if (sg_io->get_loglevel() == MSG_DEBUG)
576  test_consistency();
577 
578  return m_flow;
579 }
580 
581 bool CGraphCut::grow(GCEdge* &edge, GCNode* &current_node)
582 {
583  GCNode* node_i, *node_j;
584 
585  if ((node_i = current_node) != NULL)
586  {
587  node_i->next = NULL; // remove active flag
588 
589  if (node_i->parent == NULL)
590  {
591  node_i = NULL;
592  }
593  }
594 
595  if (node_i == NULL && (node_i = next_active()) == NULL)
596  {
597  return false;
598  }
599 
600  if (node_i->type_tree == SOURCE)
601  {
602  // grow source tree
603  for (edge = node_i->first; edge != NULL; edge = edge->next)
604  {
605  if (edge->residual_capacity)
606  {
607  node_j = edge->head;
608 
609  if (node_j->parent == NULL)
610  {
611  node_j->type_tree = SOURCE;
612  node_j->parent = edge->reverse;
613  node_j->timestamp = node_i->timestamp;
614  node_j->dist_terminal = node_i->dist_terminal + 1;
615  set_active(node_j);
616  }
617  else if (node_j->type_tree == SINK)
618  {
619  break;
620  }
621  else if (node_j->timestamp <= node_i->timestamp && node_j->dist_terminal > node_i->dist_terminal)
622  {
623  // heuristic - trying to make the distance from j to the source shorter
624  node_j->parent = edge->reverse;
625  node_j->timestamp = node_i->timestamp;
626  node_j->dist_terminal = node_i->dist_terminal + 1;
627  }
628  }
629  }
630  }
631  else
632  {
633  // grow sink tree
634  for (edge = node_i->first; edge != NULL; edge = edge->next)
635  {
636  if (edge->reverse->residual_capacity)
637  {
638  node_j = edge->head;
639 
640  if (node_j->parent == NULL)
641  {
642  node_j->type_tree = SINK;
643  node_j->parent = edge->reverse;
644  node_j->timestamp = node_i->timestamp;
645  node_j->dist_terminal = node_i->dist_terminal + 1;
646  set_active(node_j);
647  }
648  else if (node_j->type_tree == SOURCE)
649  {
650  edge = edge->reverse;
651  break;
652  }
653  else if (node_j->timestamp <= node_i->timestamp && node_j->dist_terminal > node_i->dist_terminal)
654  {
655  // heuristic - trying to make the distance from j to the sink shorter
656  node_j->parent = edge->reverse;
657  node_j->timestamp = node_i->timestamp;
658  node_j->dist_terminal = node_i->dist_terminal + 1;
659  }
660  }
661  }
662  } // grow sink tree
663 
664  if (edge != NULL)
665  {
666  node_i->next = node_i; // set active flag
667  current_node = node_i;
668  }
669  else
670  {
671  current_node = NULL;
672  }
673 
674  return true;
675 }
676 
677 void CGraphCut::augment_path(GCEdge* connecting_edge)
678 {
679  GCNode* node_i;
680  GCEdge* edge;
681  float64_t bottleneck;
682 
683  // 1. Finding bottleneck capacity
684  // 1a the source tree
685  bottleneck = connecting_edge->residual_capacity;
686 
687  for (node_i = connecting_edge->reverse->head; ; node_i = edge->head)
688  {
689  edge = node_i->parent;
690 
691  if (edge == TERMINAL_EDGE)
692  {
693  break;
694  }
695 
696  if (bottleneck > edge->reverse->residual_capacity)
697  {
698  bottleneck = edge->reverse->residual_capacity;
699  }
700  }
701 
702  if (bottleneck > node_i->tree_cap)
703  {
704  bottleneck = node_i->tree_cap;
705  }
706 
707  // 1b the sink tree
708  for (node_i = connecting_edge->head; ; node_i = edge->head)
709  {
710  edge = node_i->parent;
711 
712  if (edge == TERMINAL_EDGE)
713  {
714  break;
715  }
716 
717  if (bottleneck > edge->residual_capacity)
718  {
719  bottleneck = edge->residual_capacity;
720  }
721  }
722 
723  if (bottleneck > - node_i->tree_cap)
724  {
725  bottleneck = - node_i->tree_cap;
726  }
727 
728 
729  // 2. Augmenting
730  // 2a the source tree
731  connecting_edge->reverse->residual_capacity += bottleneck;
732  connecting_edge->residual_capacity -= bottleneck;
733 
734  for (node_i = connecting_edge->reverse->head; ; node_i = edge->head)
735  {
736  edge = node_i->parent;
737 
738  if (edge == TERMINAL_EDGE)
739  {
740  break;
741  }
742 
743  edge->residual_capacity += bottleneck;
744  edge->reverse->residual_capacity -= bottleneck;
745 
746  if (edge->reverse->residual_capacity == 0)
747  {
748  set_orphan_front(node_i); // add node_i to the beginning of the adoptation list
749  }
750  }
751 
752  node_i->tree_cap -= bottleneck;
753 
754  if (node_i->tree_cap == 0)
755  {
756  set_orphan_front(node_i); // add node_i to the beginning of the adoptation list
757  }
758 
759  // 2b the sink tree
760  for (node_i = connecting_edge->head; ; node_i = edge->head)
761  {
762  edge = node_i->parent;
763 
764  if (edge == TERMINAL_EDGE)
765  {
766  break;
767  }
768 
769  edge->reverse->residual_capacity += bottleneck;
770  edge->residual_capacity -= bottleneck;
771 
772  if (edge->residual_capacity == 0)
773  {
774  set_orphan_front(node_i);
775  }
776  }
777 
778  node_i->tree_cap += bottleneck;
779 
780  if (node_i->tree_cap == 0)
781  {
782  set_orphan_front(node_i);
783  }
784 
785  m_flow += bottleneck;
786 }
787 
788 void CGraphCut::adopt()
789 {
790  GCNodePtr* np, *np_next;
791  GCNode* node_i;
792 
793  while ((np = m_orphan_first) != NULL)
794  {
795  np_next = np->next;
796  np->next = NULL;
797 
798  while ((np = m_orphan_first) != NULL)
799  {
800  m_orphan_first = np->next;
801  node_i = np->ptr;
802  SG_FREE(np);
803 
804  if (m_orphan_first == NULL)
805  {
806  m_orphan_last = NULL;
807  }
808 
809  process_orphan(node_i, node_i->type_tree);
810  }
811 
812  m_orphan_first = np_next;
813  }
814 }
815 
816 void CGraphCut::set_orphan_front(GCNode* node_i)
817 {
818  GCNodePtr* np;
819  node_i->parent = ORPHAN_EDGE;
820  np = SG_MALLOC(GCNodePtr, 1);
821  np->ptr = node_i;
822  np->next = m_orphan_first;
823  m_orphan_first = np;
824 }
825 
826 void CGraphCut::set_orphan_rear(GCNode* node_i)
827 {
828  GCNodePtr* np;
829  node_i->parent = ORPHAN_EDGE;
830  np = SG_MALLOC(GCNodePtr, 1);
831  np->ptr = node_i;
832 
833  if (m_orphan_last != NULL)
834  {
835  m_orphan_last->next = np;
836  }
837  else
838  {
839  m_orphan_first = np;
840  }
841 
842  m_orphan_last = np;
843  np->next = NULL;
844 }
845 
846 void CGraphCut::process_orphan(GCNode* node_i, ETerminalType terminalType_tree)
847 {
848  GCNode* node_j;
849  GCEdge* edge0;
850  GCEdge* edge0_min = NULL;
851  GCEdge* edge;
852  int32_t d;
853  int32_t d_min = INFINITE_D;
854 
855  // trying to find a new parent
856  for (edge0 = node_i->first; edge0 != NULL; edge0 = edge0->next)
857  {
858  if ((terminalType_tree == SOURCE && edge0->reverse->residual_capacity) ||
859  (terminalType_tree == SINK && edge0->residual_capacity))
860  {
861  node_j = edge0->head;
862 
863  if (node_j->type_tree == terminalType_tree && (edge = node_j->parent) != NULL)
864  {
865  // check the origin of node_j
866  d = 0;
867  while (1)
868  {
869  if (node_j->timestamp == m_timestamp)
870  {
871  d += node_j->dist_terminal;
872  break;
873  }
874 
875  edge = node_j->parent;
876  d++;
877 
878  if (edge == TERMINAL_EDGE)
879  {
880  node_j->timestamp = m_timestamp;
881  node_j->dist_terminal = 1;
882  break;
883  }
884 
885  if (edge == ORPHAN_EDGE)
886  {
887  d = INFINITE_D;
888  break;
889  }
890 
891  node_j = edge->head;
892  } // while
893 
894  if (d < INFINITE_D) // node_j originates from the source, done
895  {
896  if (d < d_min)
897  {
898  edge0_min = edge0;
899  d_min = d;
900  }
901  // set marks along the path
902  for (node_j = edge0->head; node_j->timestamp != m_timestamp; node_j = node_j->parent->head)
903  {
904  node_j->timestamp = m_timestamp;
905  node_j->dist_terminal = d--;
906  }
907  }
908 
909  } // if node_j->type_tree
910  } // if(edge0->reverse->residual_capacity)
911  } // for edge0 = node_i->first
912 
913  if ((node_i->parent = edge0_min) != NULL)
914  {
915  node_i->timestamp = m_timestamp;
916  node_i->dist_terminal = d_min + 1;
917  }
918  else
919  {
920  // no parent is found, process neighbors
921  for (edge0 = node_i->first; edge0 != NULL; edge0 = edge0->next)
922  {
923  node_j = edge0->head;
924 
925  if (node_j->type_tree == terminalType_tree && (edge = node_j->parent) != NULL)
926  {
927  bool is_active_source = (terminalType_tree == SOURCE && edge0->reverse->residual_capacity);
928  bool is_active_sink = (terminalType_tree == SINK && edge0->residual_capacity);
929 
930  if (is_active_source || is_active_sink)
931  {
932  set_active(node_j);
933  }
934 
935  if (edge != TERMINAL_EDGE && edge != ORPHAN_EDGE && edge->head == node_i)
936  {
937  set_orphan_rear(node_j); // add node_j to the end of the adoptation list
938  }
939  }
940  } // for edge0 = node_i->first
941  }
942 }
943 
945 {
946  if (m_nodes[i].parent != NULL)
947  {
948  return m_nodes[i].type_tree;
949  }
950  else
951  {
952  return default_terminal;
953  }
954 }
955 
957 {
958  // print SOURCE-node_i and node_i->SINK edges
959  for (int32_t i = 0; i < m_num_nodes; i++)
960  {
961  GCNode* node_i = m_nodes + i;
962  if (node_i->parent == TERMINAL_EDGE)
963  {
964  if (node_i->type_tree == SOURCE)
965  {
966  SG_SPRINT("\n s -> %d, cost = %f", node_i->id, node_i->tree_cap);
967  }
968  else
969  {
970  SG_SPRINT("\n %d -> t, cost = %f", node_i->id, node_i->tree_cap);
971  }
972  }
973  }
974 
975  // print node_i->node_j edges
976  for (int32_t i = 0; i < m_num_edges; i++)
977  {
978  GCEdge* edge = m_edges + i;
979  SG_SPRINT("\n %d -> %d, cost = %f", edge->reverse->head->id, edge->head->id, edge->residual_capacity);
980  }
981 
982 }
983 
985 {
986  for (int32_t i = 0; i < m_num_nodes; i++)
987  {
988  GCNode* node_i = m_nodes + i;
989 
990  if (get_assignment(i) == SOURCE)
991  {
992  SG_SPRINT("\nGCNode %2d: S", node_i->id);
993  }
994  else
995  {
996  SG_SPRINT("\nGCNode %2d: T", node_i->id);
997  }
998  }
999 }
1000 
1001 void CGraphCut::test_consistency(GCNode* current_node)
1002 {
1003  GCNode* node_i;
1004  GCEdge* edge;
1005  int32_t num1 = 0;
1006  int32_t num2 = 0;
1007 
1008  // test whether all nodes i with i->next!=NULL are indeed in the queue
1009  for (int32_t i = 0; i < m_num_nodes; i++)
1010  {
1011  node_i = m_nodes + i;
1012  if (node_i->next || node_i == current_node)
1013  {
1014  num1++;
1015  }
1016  }
1017 
1018  for (int32_t r = 0; r < 3; r++)
1019  {
1020  node_i = (r == 2) ? current_node : m_active_first[r];
1021 
1022  if (node_i)
1023  {
1024  for (; ; node_i = node_i->next)
1025  {
1026  num2++;
1027 
1028  if (node_i->next == node_i)
1029  {
1030  if (r < 2)
1031  ASSERT(node_i == m_active_last[r])
1032  else
1033  ASSERT(node_i == current_node)
1034 
1035  break;
1036  }
1037  }
1038  }
1039  }
1040 
1041  ASSERT(num1 == num2);
1042 
1043  for (int32_t i = 0; i < m_num_nodes; i++)
1044  {
1045  node_i = m_nodes + i;
1046 
1047  // test whether all edges in seach trees are non-saturated
1048  if (node_i->parent == NULL) {}
1049  else if (node_i->parent == ORPHAN_EDGE) {}
1050  else if (node_i->parent == TERMINAL_EDGE)
1051  {
1052  if (node_i->type_tree == SOURCE)
1053  ASSERT(node_i->tree_cap > 0)
1054  else
1055  ASSERT(node_i->tree_cap < 0)
1056  }
1057  else
1058  {
1059  if (node_i->type_tree == SOURCE)
1060  ASSERT(node_i->parent->reverse->residual_capacity > 0)
1061  else
1062  ASSERT(node_i->parent->residual_capacity > 0)
1063  }
1064 
1065  // test whether passive nodes in search trees have neighbors in
1066  // a different tree through non-saturated edges
1067  if (node_i->parent && !node_i->next)
1068  {
1069  if (node_i->type_tree == SOURCE)
1070  {
1071  ASSERT(node_i->tree_cap >= 0);
1072 
1073  for (edge = node_i->first; edge; edge = edge->next)
1074  {
1075  if (edge->residual_capacity > 0)
1076  {
1077  ASSERT(edge->head->parent && edge->head->type_tree == SOURCE);
1078  }
1079  }
1080  }
1081  else
1082  {
1083  ASSERT(node_i->tree_cap <= 0);
1084 
1085  for (edge = node_i->first; edge; edge = edge->next)
1086  {
1087  if (edge->reverse->residual_capacity > 0)
1088  {
1089  ASSERT(edge->head->parent && (edge->head->type_tree == SINK));
1090  }
1091  }
1092  }
1093  }
1094  // test marking invariants
1095  if (node_i->parent && node_i->parent != ORPHAN_EDGE && node_i->parent != TERMINAL_EDGE)
1096  {
1097  ASSERT(node_i->timestamp <= node_i->parent->head->timestamp);
1098 
1099  if (node_i->timestamp == node_i->parent->head->timestamp)
1100  {
1101  ASSERT(node_i->dist_terminal > node_i->parent->head->dist_terminal);
1102  }
1103  }
1104  }
1105 }
float64_t tree_cap
Definition: GraphCut.h:82
virtual const char * get_name() const
Definition: GraphCut.h:133
float64_t evaluate_energy(const SGVector< int32_t > state) const
float64_t residual_capacity
Definition: GraphCut.h:55
GCEdge * parent
Definition: GraphCut.h:68
#define INFINITE_D
Definition: GraphCut.h:27
#define ORPHAN_EDGE
Definition: GraphCut.h:25
const SGVector< int32_t > get_variables() const
Definition: Factor.cpp:107
Graph cuts node.
Definition: GraphCut.h:61
int32_t get_num_factors() const
void print_assignment()
Definition: GraphCut.cpp:984
ETerminalType
Definition: GraphCut.h:32
float64_t m_map_energy
Definition: GraphCut.h:307
CDynamicObjectArray * get_factors() const
#define SG_ERROR(...)
Definition: SGIO.h:128
#define REQUIRE(x,...)
Definition: SGIO.h:181
void add_edge(int32_t i, int32_t j, float64_t capacity, float64_t reverse_capacity)
Definition: GraphCut.cpp:451
GCEdge * next
Definition: GraphCut.h:51
Graph guts node pointer.
Definition: GraphCut.h:88
Class CMAPInferImpl abstract class of MAP inference implementation.
Definition: MAPInference.h:98
ETerminalType type_tree
Definition: GraphCut.h:78
void build_st_graph(int32_t num_nodes, int32_t num_edges)
Definition: GraphCut.cpp:122
SGVector< float64_t > get_energies() const
Definition: Factor.cpp:169
GCEdge * reverse
Definition: GraphCut.h:53
GCNode * ptr
Definition: GraphCut.h:91
GCNodePtr * next
Definition: GraphCut.h:93
float64_t compute_maxflow()
Definition: GraphCut.cpp:540
CFactorGraph * m_fg
Definition: MAPInference.h:128
#define SG_SPRINT(...)
Definition: SGIO.h:165
#define ASSERT(x)
Definition: SGIO.h:176
void add_tweights(int32_t i, float64_t cap_source, float64_t cap_sink)
Definition: GraphCut.cpp:431
SGIO * sg_io
Definition: init.cpp:44
GCNode * next
Definition: GraphCut.h:72
int32_t size() const
Definition: SGVector.h:156
int32_t get_num_vars() const
double float64_t
Definition: common.h:60
int32_t dist_terminal
Definition: GraphCut.h:76
int32_t id
Definition: GraphCut.h:64
Dynamic array class for CSGObject pointers that creates an array that can be used like a list or an a...
GCEdge * first
Definition: GraphCut.h:66
EMessageType get_loglevel() const
Definition: SGIO.cpp:189
#define TERMINAL_EDGE
Definition: GraphCut.h:24
const int32_t get_num_vars() const
Definition: Factor.cpp:112
#define SG_UNREF(x)
Definition: SGObject.h:53
#define SG_DEBUG(...)
Definition: SGIO.h:106
Graph cuts edge.
Definition: GraphCut.h:44
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
Class CFactorGraph a factor graph is a structured input in general.
Definition: FactorGraph.h:27
virtual ~CGraphCut()
Definition: GraphCut.cpp:42
CSGObject * get_element(int32_t index) const
SGVector< int32_t > get_cardinalities() const
int32_t id
Definition: GraphCut.h:47
#define SG_UNSTABLE(func,...)
Definition: SGIO.h:131
int32_t timestamp
Definition: GraphCut.h:74
Class CFactor A factor is defined on a clique in the factor graph. Each factor can have its own data...
Definition: Factor.h:89
GCNode * head
Definition: GraphCut.h:49
const SGVector< int32_t > get_cardinalities() const
Definition: Factor.cpp:122
virtual float64_t inference(SGVector< int32_t > assignment)
Definition: GraphCut.cpp:193
ETerminalType get_assignment(int32_t i, ETerminalType default_terminal=SOURCE)
Definition: GraphCut.cpp:944

SHOGUN Machine Learning Toolbox - Documentation