Limbo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
LPColoringOld.h
Go to the documentation of this file.
1 
17 
19 #ifndef LIMBO_ALGORITHMS_COLORING_LP
20 #define LIMBO_ALGORITHMS_COLORING_LP
21 
22 #include <iostream>
23 #include <vector>
24 #include <set>
25 #include <cassert>
26 #include <cmath>
27 #include <stdlib.h>
28 #include <cstdio>
29 #include <sstream>
30 #include <boost/cstdint.hpp>
31 #include <boost/unordered_map.hpp>
32 #include <boost/graph/graph_concepts.hpp>
33 #include <boost/graph/subgraph.hpp>
34 #include <boost/property_map/property_map.hpp>
35 #include <boost/graph/bipartite.hpp>
36 #include "gurobi_c++.h"
37 
38 using std::vector;
39 using std::set;
40 using std::cin;
41 using std::cout;
42 using std::endl;
43 using std::ofstream;
44 using std::pair;
45 using std::ostringstream;
46 using std::string;
47 using boost::unordered_map;
48 using boost::uint32_t;
49 using boost::int32_t;
50 
52 namespace limbo
53 {
55 namespace algorithms
56 {
58 namespace coloring
59 {
60 
61 template <typename GraphType>
62 class LPColoring
63 {
64  public:
65  typedef GraphType graph_type;
66  //typedef boost::subgraph_type<graph_type> subgraph_type;
67  typedef typename boost::graph_traits<graph_type>::vertex_descriptor graph_vertex_type;
68  typedef typename boost::graph_traits<graph_type>::edge_descriptor graph_edge_type;
72  typedef typename boost::property_map<graph_type, boost::edge_weight_t>::const_type edge_weight_map_type;
74  typedef typename boost::property_map<graph_type, boost::vertex_color_t>::const_type vertex_color_map_type;
75 
76  enum RoundingScheme
77  {
78  DIRECT_ILP, // directly solve ILP to get optimal solution
79  FIXED_ILP, // this rounding scheme is not in use now, because feasible solution is not guaranteed
80  ITERATIVE_ILP,
81  GREEDY,
82  POST_ILP
83  };
84  enum ColorNum
85  {
86  THREE,
87  FOUR
88  };
89  // hasher class for graph_edge_type
90  struct edge_hash_type : std::unary_function<graph_edge_type, std::size_t>
91  {
92  std::size_t operator()(graph_edge_type const& e) const
93  {
94  std::size_t seed = 0;
95  boost::hash_combine(seed, e.m_source);
96  boost::hash_combine(seed, e.m_target);
97  return seed;
98  }
99  };
100 
103  LPColoring(graph_type const& g);
105  ~LPColoring() {};
106 
108  void operator()() {this->graphColoring();}
109 
110  bool conflictCost() const {return m_conflict_cost;}
111  void conflictCost(bool flag) {m_conflict_cost = flag;}
112  double stitchWeight() const {return m_stitch_weight;}
113  void stitchWeight(double w) {m_stitch_weight = w;}
114  RoundingScheme roundingScheme() const {return m_rounding_scheme;}
115  void roundingScheme(RoundingScheme rs) {m_rounding_scheme = rs;}
116  ColorNum colorNum() const {return m_color_num;}
117  void colorNum(ColorNum cn) {m_color_num = cn;}
118 
120  void oddCycles(graph_vertex_type& v);
121 
123  void graphColoring();
128  void ILPColoring(vector<GRBVar>& coloringBits, vector<GRBVar>& vEdgeBit);
129 
131  void rounding_bindingAnalysis(GRBModel& opt_model, vector<GRBVar>& coloringBits, vector<GRBVar>& vEdgeBit);
132 
138  void rounding_ILP(GRBModel& opt_model, vector<GRBVar>& coloringBits, vector<GRBVar>& vEdgeBit);
140  void rounding_Greedy_v0(vector<GRBVar>& coloringBits);
141  void rounding_Greedy(vector<GRBVar>& coloringBits);
142 
144  uint32_t vertexColor(graph_vertex_type& node) const;
145  uint32_t conflictNum() const;
146  uint32_t stitchNum() const;
148  uint32_t lpConflictNum() const;
149  uint32_t lpStitchNum() const;
151  void store_lp_coloring(vector<GRBVar>& coloringBits);
152 
153  pair<uint32_t, uint32_t> nonIntegerNum(const vector<GRBVar>& coloringBits, const vector<GRBVar>& vEdgeBit) const;
155  void printBoolVec(vector<bool>& vec) const;
156  void printBoolArray(bool* vec, uint32_t vec_size) const;
158  void write_graph_dot(graph_vertex_type& v) const;
159  void write_graph_dot() const;
160  void write_graph_color() const;
161 
164  const uint32_t COLORING_BITS;
165  protected:
167  void setMap();
169  bool isInteger(double value) const {return value == floor(value);}
171  bool sameColor(graph_vertex_type v1, graph_vertex_type v2) const;
173  bool lpSameColor(graph_vertex_type v1, graph_vertex_type v2) const;
175  bool hasConflict(graph_edge_type e) const;
177  bool hasLPConflict(graph_edge_type e) const;
179  bool hasStitch(graph_edge_type e) const;
181  bool hasLPStitch(graph_edge_type e) const;
182 
185  graph_type const& m_graph;
187  unordered_map<graph_vertex_type, uint32_t> m_vertex_map;
188  unordered_map<uint32_t, graph_vertex_type> m_inverse_vertex_map;
190  unordered_map<graph_edge_type, uint32_t, edge_hash_type> m_edge_map;
192  edge_weight_map_type m_edge_weight_map;
194  vertex_color_map_type m_vertex_color_map;
196  unordered_map<graph_vertex_type, uint32_t> m_coloring;
198  vector<double> m_lp_coloring;
200  vector< vector<graph_vertex_type> > m_odd_cycles;
203  bool m_conflict_cost;
205  double m_stitch_weight;
207  RoundingScheme m_rounding_scheme;
209  ColorNum m_color_num;
211  uint32_t m_lp_iter_cnt;
213  uint32_t m_ilp_iter_cnt;
214 
215  //linear programming constraints property
216  unordered_map<string, bool> m_stitch_constrs;
217  uint32_t m_constrs_num;
218 };
219 
221 template <typename GraphType>
222 LPColoring<GraphType>::LPColoring(graph_type const& g) : COLORING_BITS(2), m_graph(g)
223 {
224  m_conflict_cost = true;
225  m_stitch_weight = 0.1;
226  m_rounding_scheme = FIXED_ILP;
227  m_color_num = THREE;
228  m_lp_iter_cnt = 0;
229  m_ilp_iter_cnt = 0;
230  m_constrs_num = 0;
231  this->setMap();
232 }
233 
234 //set the vertex map
235 template<typename GraphType>
236 void LPColoring<GraphType>::setMap()
237 {
238  // build the vertex-index map
239  pair<typename graph_type::vertex_iterator, typename graph_type::vertex_iterator> vertex_range = vertices(m_graph);
240  m_vertex_map.clear();
241  m_inverse_vertex_map.clear();
242  uint32_t vertex_num = 0;
243  for(BOOST_AUTO(itr, vertex_range.first); itr != vertex_range.second; ++itr)
244  {
245  m_vertex_map[*itr] = vertex_num;
246  m_inverse_vertex_map[vertex_num] = *itr;
247  ++vertex_num;
248  }
249  // build edge-index map
250  pair<typename graph_type::edge_iterator, typename graph_type::edge_iterator> edge_range = edges(m_graph);
251  m_edge_map.clear();
252  uint32_t edge_num = 0;
253  for (BOOST_AUTO(itr, edge_range.first); itr != edge_range.second; ++itr)
254  {
255  m_edge_map[*itr] = edge_num;
256  ++edge_num;
257  }
258  // build edge weight map and vertex color map
259  m_edge_weight_map = get(boost::edge_weight, m_graph);
260  m_vertex_color_map = get(boost::vertex_color, m_graph);
261 }
262 
263 //DFS to search for the odd cycles, stored in m_odd_cycles
264 template<typename GraphType>
265 void LPColoring<GraphType>::oddCycles(graph_vertex_type& v)
266 {
267  //if(m_vertex_map.empty() || m_inverse_vertex_map.empty()) this->setMap();
268 #ifdef DEBUG_LPCOLORING
269  assert(m_vertex_map.find(v) != m_vertex_map.end());
270 #endif
271  //odd_cycle results
272  this->m_odd_cycles.clear();
273 
274  //the array denoting the distancevisiting of the graph
275  uint32_t vertex_num = num_edges(m_graph);
276  int32_t nodeDist[vertex_num];
277  bool nodeVisited[vertex_num];
278  for(uint32_t k = 0; k < vertex_num; k++)
279  {
280  nodeDist[k] = -1;
281  nodeVisited[k] = false;
282  }
283 
284  //inCycle flag
285  bool inCycle[vertex_num];
286  for(uint32_t k = 0; k < vertex_num; k++) inCycle[k] = false;
287 
288  //dfs_stack for DFS
289  vector<graph_vertex_type> dfs_stack;
290  dfs_stack.reserve(vertex_num);
291  uint32_t v_index = m_vertex_map[v];
292  nodeVisited[v_index] = true;
293  nodeDist[v_index] = 0;
294  dfs_stack.push_back(v);
295  while(false == dfs_stack.empty())
296  {
297  //determine whether the top element needs to be popped
298  bool popFlag = true;
299  //access the top element on the dfs_stack
300  graph_vertex_type curr_v = dfs_stack.back();
301  uint32_t curr_index = m_vertex_map[curr_v];
302  //access the neighbors
303  typename boost::graph_traits<graph_type>::adjacency_iterator vi_begin, vi_end;
304  boost::tie(vi_begin, vi_end) = adjacent_vertices(curr_v, m_graph);
305  for(BOOST_AUTO(vi, vi_begin); vi != vi_end; ++vi)
306  {
307  BOOST_AUTO(found_edge, boost::edge(curr_v, *vi, m_graph));
308 #ifdef DEBUG_LPCOLORING
309  assert(found_edge.second);
310 #endif
311  // skip stitch edges
312  if (m_edge_weight_map[found_edge.first] < 0) continue;
313 
314  uint32_t next_index = m_vertex_map[*vi];
315  if(nodeDist[next_index] < 0)
316  {
317  nodeDist[next_index] = 1 - nodeDist[curr_index];
318  nodeVisited[next_index] = true;
319  //push to the dfs_stack
320  dfs_stack.push_back(*vi);
321  popFlag = false;
322  break;
323  }
324  } //end for
325 
326  if(true == popFlag)
327  {
328  //detect the odd cycle
329  for(BOOST_AUTO(vi, vi_begin); vi != vi_end; ++vi)
330  {
331  uint32_t next_index = m_vertex_map[*vi];
332  if(true == nodeVisited[next_index] && (nodeDist[next_index] == nodeDist[curr_index]))
333  {
334  //suppose v/v_index is not in the current cycle
335  inCycle[v_index] = true;
336  //detect the cycle between curr_v and *vi
337  vector<graph_vertex_type> cycle;
338  int cnt = dfs_stack.size();
339  for(int k = cnt - 1; k >= 0; k--)
340  {
341  cycle.push_back(dfs_stack[k]);
342  //flag the nodes in cycle
343  inCycle[m_vertex_map[dfs_stack[k]]] = true;
344  if(dfs_stack[k] == (*vi)) break;
345  }
346  //store the cycle, when v/v_index is in cycle
347  if(cycle.size() > 0 && inCycle[v_index])
348  {
349  this->m_odd_cycles.push_back(cycle);
350  }
351  else if(cycle.size() == 0)
352  {
353  cout << "ERROR: the cycle detected contains no nodes" << endl;
354  }
355  }//end if
356  }//end for vi
357 
358  //pop the top element
359  dfs_stack.pop_back();
360  nodeVisited[curr_index] = false;
361  }//end if popFlag
362 
363  }//end while
364 #ifdef DEBUG_LPCOLORING
365  if(m_odd_cycles.size() > 0) cout << m_odd_cycles.size() << " cycles detected." << endl;
366  else cout << "No cycles detected." << endl;
367 #endif
368 }
369 
370 //relaxed linear programming based coloring for the conflict graph (m_graph)
371 template<typename GraphType>
372 void LPColoring<GraphType>::graphColoring()
373 {
374  //build the vertex-index map
375  //pair<typename graph_type::vertex_iterator, typename graph_type::vertex_iterator> vertex_range = vertices(m_graph);
376 #if 0
377  //unordered_map<graph_vertex_type, uint32_t> m_vertex_map;
378  //unordered_map<uint32_t, graph_vertex_type> m_inverse_vertex_map;
379  m_vertex_map.clear(), m_inverse_vertex_map.clear();
380  uint32_t vertex_num = 0;
381  for(BOOST_AUTO(itr, vertex_range.first); itr != vertex_range.second; ++itr)
382  {
383  m_vertex_map[*itr] = vertex_num;
384  m_inverse_vertex_map[vertex_num] = *itr;
385  vertex_num++;
386  }
387 #endif
388  cout << "---------------------------------------Iterative LP based coloring scheme------------------------------------\n";
389  // reset iteration counters
390  m_lp_iter_cnt = m_ilp_iter_cnt = 0;
391 
392  uint32_t vertex_num = boost::num_vertices(m_graph);
393  uint32_t edge_num = boost::num_edges(m_graph);
394  uint32_t coloring_bits_num = COLORING_BITS*vertex_num;
395  //uint32_t variable_num = coloring_bits_num+edge_num;
396 
397  //set up the LP environment
398  GRBEnv env = GRBEnv();
399  //mute the log from the LP solver
400  //env.set(GRB_IntParam_OutputFlag, 0);
401  GRBModel opt_model = GRBModel(env);
402  //set up the LP variables
403  vector<GRBVar> coloringBits;
404  vector<GRBVar> vEdgeBit;
405  // vertex and edge variables
406  coloringBits.reserve(coloring_bits_num);
407  for (uint32_t i = 0; i != coloring_bits_num; ++i)
408  {
409  ostringstream oss;
410  oss << "v" << i;
411  if (roundingScheme() == DIRECT_ILP)
412  coloringBits.push_back(opt_model.addVar(0.0, 1.0, 0.0, GRB_INTEGER, oss.str()));
413  else
414  coloringBits.push_back(opt_model.addVar(0.0, 1.0, 0.0, GRB_CONTINUOUS, oss.str()));
415  }
416  vEdgeBit.reserve(edge_num);
417  for (uint32_t i = 0; i != edge_num; ++i)
418  {
419  // some variables here may not be used
420  ostringstream oss;
421  oss << "e" << i;
422  vEdgeBit.push_back(opt_model.addVar(0.0, 1.0, 0.0, GRB_CONTINUOUS, oss.str()));
423  }
424  // conflict edge variables
425  // they are used in objective function to minimize conflict cost
426  pair<typename graph_type::edge_iterator, typename graph_type::edge_iterator> edge_range = edges(m_graph);
427  //Integrate new variables
428  opt_model.update();
429 
430  //set up the objective
431  GRBLinExpr obj_init (0);
432  GRBLinExpr obj (0); // maybe useful
433  if (this->conflictCost())
434  {
435  for (BOOST_AUTO(itr, edge_range.first); itr != edge_range.second; ++itr)
436  {
437  BOOST_AUTO(w, m_edge_weight_map[*itr]);
438  uint32_t edge_idx = m_edge_map[*itr];
439  if (w >= 0) // conflict
440  obj_init += vEdgeBit[edge_idx];
441  }
442  }
443  for (BOOST_AUTO(itr, edge_range.first); itr != edge_range.second; ++itr)
444  {
445  BOOST_AUTO(w, m_edge_weight_map[*itr]);
446  uint32_t edge_idx = m_edge_map[*itr];
447  if (w < 0) // stitch
448  obj_init += m_stitch_weight*vEdgeBit[edge_idx];
449  }
450  obj = obj_init;
451  opt_model.setObjective(obj, GRB_MINIMIZE);
452 
453  //set up the LP constraints
454  //typename graph_type::edges_size_type edge_num = num_edges(m_graph);
455  //typename graph_type::vertices_size_type vertex_num = num_vertices(m_graph);
456  for(BOOST_AUTO(itr, edge_range.first); itr != edge_range.second; ++itr)
457  {
458  BOOST_AUTO(from, source(*itr, m_graph));
459  uint32_t f_ind = m_vertex_map[from];
460  BOOST_AUTO(to, target(*itr, m_graph));
461  uint32_t t_ind = m_vertex_map[to];
462  //coloring conflict constraints
463  uint32_t vertex_idx1 = f_ind<<1;
464  uint32_t vertex_idx2 = t_ind<<1;
465 
466  BOOST_AUTO(w, m_edge_weight_map[*itr]);
467  GRBConstr tmpConstr;
468  char buf[100];
469  string tmpConstr_name;
470  if (w >= 0) // constraints for conflict edges
471  {
472  if (!conflictCost())
473  {
474  sprintf(buf, "%u", m_constrs_num++);
475  tmpConstr_name = "R" + string(buf);
476  tmpConstr = opt_model.addConstr(
477  coloringBits[vertex_idx1] + coloringBits[vertex_idx1+1]
478  + coloringBits[vertex_idx2] + coloringBits[vertex_idx2+1] >= 1
479  , tmpConstr_name);
480 #ifdef DEBUG_LPCOLORING
481  assert(m_stitch_constrs.find(tmpConstr_name) == m_stitch_constrs.end());
482 #endif
483  m_stitch_constrs[tmpConstr_name] = false;
484 
485  sprintf(buf, "%u", m_constrs_num++);
486  tmpConstr_name = "R" + string(buf);
487  tmpConstr = opt_model.addConstr(
488  1 - coloringBits[vertex_idx1] + coloringBits[vertex_idx1+1]
489  + 1 - coloringBits[vertex_idx2] + coloringBits[vertex_idx2+1] >= 1
490  , tmpConstr_name);
491 #ifdef DEBUG_LPCOLORING
492  assert(m_stitch_constrs.find(tmpConstr_name) == m_stitch_constrs.end());
493 #endif
494  m_stitch_constrs[tmpConstr_name] = false;
495 
496  sprintf(buf, "%u", m_constrs_num++);
497  tmpConstr_name = "R" + string(buf);
498  tmpConstr = opt_model.addConstr(
499  coloringBits[vertex_idx1] + 1 - coloringBits[vertex_idx1+1]
500  + coloringBits[vertex_idx2] + 1 - coloringBits[vertex_idx2+1] >= 1
501  , tmpConstr_name);
502 #ifdef DEBUG_LPCOLORING
503  assert(m_stitch_constrs.find(tmpConstr_name) == m_stitch_constrs.end());
504 #endif
505  m_stitch_constrs[tmpConstr_name] = false;
506 
507  sprintf(buf, "%u", m_constrs_num++);
508  tmpConstr_name = "R" + string(buf);
509  tmpConstr = opt_model.addConstr(
510  1 - coloringBits[vertex_idx1] + 1 - coloringBits[vertex_idx1+1]
511  + 1 - coloringBits[vertex_idx2] + 1 - coloringBits[vertex_idx2+1] >= 1
512  , tmpConstr_name);
513 #ifdef DEBUG_LPCOLORING
514  assert(m_stitch_constrs.find(tmpConstr_name) == m_stitch_constrs.end());
515 #endif
516  m_stitch_constrs[tmpConstr_name] = false;
517 
518  }
519  else
520  {
521  uint32_t edge_idx = m_edge_map[*itr];
522 
523  sprintf(buf, "%u", m_constrs_num++);
524  tmpConstr_name = "R" + string(buf);
525  tmpConstr = opt_model.addConstr(
526  coloringBits[vertex_idx1] + coloringBits[vertex_idx1+1]
527  + coloringBits[vertex_idx2] + coloringBits[vertex_idx2+1]
528  + vEdgeBit[edge_idx] >= 1
529  , tmpConstr_name);
530 #ifdef DEBUG_LPCOLORING
531  assert(m_stitch_constrs.find(tmpConstr_name) == m_stitch_constrs.end());
532 #endif
533  m_stitch_constrs[tmpConstr_name] = false;
534 
535  sprintf(buf, "%u", m_constrs_num++);
536  tmpConstr_name = "R" + string(buf);
537  tmpConstr = opt_model.addConstr(
538  1 - coloringBits[vertex_idx1] + coloringBits[vertex_idx1+1]
539  + 1 - coloringBits[vertex_idx2] + coloringBits[vertex_idx2+1]
540  + vEdgeBit[edge_idx] >= 1
541  , tmpConstr_name);
542 #ifdef DEBUG_LPCOLORING
543  assert(m_stitch_constrs.find(tmpConstr_name) == m_stitch_constrs.end());
544 #endif
545  m_stitch_constrs[tmpConstr_name] = false;
546 
547  sprintf(buf, "%u", m_constrs_num++);
548  tmpConstr_name = "R" + string(buf);
549  tmpConstr = opt_model.addConstr(
550  coloringBits[vertex_idx1] + 1 - coloringBits[vertex_idx1+1]
551  + coloringBits[vertex_idx2] + 1 - coloringBits[vertex_idx2+1]
552  + vEdgeBit[edge_idx] >= 1
553  , tmpConstr_name);
554 #ifdef DEBUG_LPCOLORING
555  assert(m_stitch_constrs.find(tmpConstr_name) == m_stitch_constrs.end());
556 #endif
557  m_stitch_constrs[tmpConstr_name] = false;
558 
559  sprintf(buf, "%u", m_constrs_num++);
560  tmpConstr_name = "R" + string(buf);
561  tmpConstr = opt_model.addConstr(
562  1 - coloringBits[vertex_idx1] + 1 - coloringBits[vertex_idx1+1]
563  + 1 - coloringBits[vertex_idx2] + 1 - coloringBits[vertex_idx2+1]
564  + vEdgeBit[edge_idx] >= 1
565  , tmpConstr_name);
566 #ifdef DEBUG_LPCOLORING
567  assert(m_stitch_constrs.find(tmpConstr_name) == m_stitch_constrs.end());
568 #endif
569  m_stitch_constrs[tmpConstr_name] = false;
570 
571  }
572  }
573  else // constraints for stitch edges
574  {
575  uint32_t edge_idx = m_edge_map[*itr];
576 
577  sprintf(buf, "%u", m_constrs_num++);
578  tmpConstr_name = "R" + string(buf);
579  tmpConstr = opt_model.addConstr(
580  coloringBits[vertex_idx1] - coloringBits[vertex_idx2] - vEdgeBit[edge_idx] <= 0
581  , tmpConstr_name);
582 #ifdef DEBUG_LPCOLORING
583  assert(m_stitch_constrs.find(tmpConstr_name) == m_stitch_constrs.end());
584 #endif
585  m_stitch_constrs[tmpConstr_name] = true;
586 
587  sprintf(buf, "%u", m_constrs_num++);
588  tmpConstr_name = "R" + string(buf);
589  tmpConstr = opt_model.addConstr(
590  coloringBits[vertex_idx2] - coloringBits[vertex_idx1] - vEdgeBit[edge_idx] <= 0
591  , tmpConstr_name);
592 #ifdef DEBUG_LPCOLORING
593  assert(m_stitch_constrs.find(tmpConstr_name) == m_stitch_constrs.end());
594 #endif
595  m_stitch_constrs[tmpConstr_name] = true;
596 
597  sprintf(buf, "%u", m_constrs_num++);
598  tmpConstr_name = "R" + string(buf);
599  tmpConstr = opt_model.addConstr(
600  coloringBits[vertex_idx1+1] - coloringBits[vertex_idx2+1] - vEdgeBit[edge_idx] <= 0
601  , tmpConstr_name);
602 #ifdef DEBUG_LPCOLORING
603  assert(m_stitch_constrs.find(tmpConstr_name) == m_stitch_constrs.end());
604 #endif
605  m_stitch_constrs[tmpConstr_name] = true;
606 
607  sprintf(buf, "%u", m_constrs_num++);
608  tmpConstr_name = "R" + string(buf);
609  tmpConstr = opt_model.addConstr(
610  coloringBits[vertex_idx2+1] - coloringBits[vertex_idx1+1] - vEdgeBit[edge_idx] <= 0
611  , tmpConstr_name);
612 #ifdef DEBUG_LPCOLORING
613  assert(m_stitch_constrs.find(tmpConstr_name) == m_stitch_constrs.end());
614 #endif
615  m_stitch_constrs[tmpConstr_name] = true;
616 
617  }
618  }
619 
620  if (colorNum() == THREE)
621  {
622  GRBConstr tmpConstr;
623  string tmpConstr_name;
624  char buf[100];
625  for(uint32_t k = 0; k < coloring_bits_num; k += COLORING_BITS)
626  {
627  sprintf(buf, "%u", m_constrs_num++);
628  tmpConstr_name = "R" + string(buf);
629  tmpConstr = opt_model.addConstr(coloringBits[k] + coloringBits[k+1] <= 1, tmpConstr_name);
630 #ifdef DEBUG_LPCOLORING
631  assert(m_stitch_constrs.find(tmpConstr_name) == m_stitch_constrs.end());
632 #endif
633  m_stitch_constrs[tmpConstr_name] = false;
634  }
635  }
636  //optimize model
637  opt_model.update();
638  cout << "================== LP iteration #" << m_lp_iter_cnt++ << " ==================\n";
639  opt_model.optimize();
640 #ifdef DEBUG_LPCOLORING
641  opt_model.write("graph.lp");
642  opt_model.write("graph.sol");
643 #endif
644  int optim_status = opt_model.get(GRB_IntAttr_Status);
645  if(optim_status == GRB_INFEASIBLE)
646  {
647  cout << "ERROR: The model is infeasible... EXIT" << endl;
648  exit(1);
649  }
650  cout << endl;
651 
652  //iteratively scheme
653  while(true)
654  {
655  uint32_t non_integer_num, half_integer_num;
656  BOOST_AUTO(pair, this->nonIntegerNum(coloringBits, vEdgeBit));
657  non_integer_num = pair.first;
658  half_integer_num = pair.second;
659  if(non_integer_num == 0) break;
660  //store the lp coloring results
661  this->store_lp_coloring(coloringBits);
662  /*
663  m_lp_coloring.clear();
664  m_lp_coloring.reserve(coloring_bits_num);
665  for(uint32_t k = 0; k < coloring_bits_num; k++)
666  {
667  double value = coloringBits[k].get(GRB_DoubleAttr_X);
668  m_lp_coloring.push_back(value);
669  }
670  */
671 #ifdef DEBUG_LPCOLORING
672  // compare unrounded conflict number is more fair
673  this->lpConflictNum();
674  this->write_graph_dot();
675 #endif
676 
677  vector<bool> non_integer_flag(coloring_bits_num, false);
678  //set the new objective
679  //push the non-half_integer to 0/1
680  // only check for vertices
681 
682  for(uint32_t k = 0; k < coloring_bits_num; k += COLORING_BITS)
683  {
684  double value1 = coloringBits[k].get(GRB_DoubleAttr_X);
685  double value2 = coloringBits[k+1].get(GRB_DoubleAttr_X);
686  if (!isInteger(value1) || !isInteger(value2))
687  {
688  if (value1 > value2)
689  {
690  if (!conflictCost())
691  obj += coloringBits[k+1]-coloringBits[k];
692  }
693  else if (value1 < value2)
694  {
695  if (!conflictCost())
696  obj += coloringBits[k]-coloringBits[k+1];
697  }
698  if (value1 == 0.5) non_integer_flag[k] = true;
699  if (value2 == 0.5) non_integer_flag[k+1] = true;
700  }
701 #if 0
702  if(value < 0.5)
703  {
704  if(!conflictCost())
705  obj = obj + 2*coloringBits[k];
706  non_integer_flag[k] = false;
707  }
708  else if (value > 0.5)
709  {
710  if(!conflictCost())
711  obj = obj - 2*coloringBits[k];
712  non_integer_flag[k] = false;
713  }
714  else
715  {
716  non_integer_flag[k] = true;
717  }
718 #endif
719  }//end for
720 
721  //minimize the conflict number
722  if (!conflictCost())
723  {
724  for(BOOST_AUTO(itr, edge_range.first); itr != edge_range.second; ++itr)
725  {
726  BOOST_AUTO(from, source(*itr, m_graph));
727  uint32_t f_ind = m_vertex_map[from];
728  BOOST_AUTO(to, target(*itr, m_graph));
729  uint32_t t_ind = m_vertex_map[to];
730  uint32_t vertex_idx1 = f_ind<<1;
731  uint32_t vertex_idx2 = t_ind<<1;
732  if (m_lp_coloring[vertex_idx1] > m_lp_coloring[vertex_idx2])
733  {
734  obj += coloringBits[vertex_idx2] - coloringBits[vertex_idx1];
735  }
736  else if (m_lp_coloring[vertex_idx1] < m_lp_coloring[vertex_idx2])
737  {
738  obj += coloringBits[vertex_idx1] - coloringBits[vertex_idx2];
739  }
740 
741  vertex_idx1 += 1;
742  vertex_idx2 += 1;
743  if (m_lp_coloring[vertex_idx1] > m_lp_coloring[vertex_idx2])
744  {
745  obj += coloringBits[vertex_idx2] - coloringBits[vertex_idx1];
746  }
747  else if (m_lp_coloring[vertex_idx1] < m_lp_coloring[vertex_idx2])
748  {
749  obj += coloringBits[vertex_idx1] - coloringBits[vertex_idx2];
750  }
751  }//end for
752  }
753 
754  opt_model.setObjective(obj);
755 
756  //add the new constraints
757  //odd cycle trick from Prof. Baldick
758  for(uint32_t k = 0; k < coloring_bits_num; k++)
759  {
760  //if the current bit is already handled
761  if(false == non_integer_flag[k]) continue;
762 
763  uint32_t vertex_index = k/COLORING_BITS;
764 #ifdef DEBUG_LPCOLORING
765  assert(m_inverse_vertex_map.find(vertex_index) != m_inverse_vertex_map.end());
766 #endif
767  graph_vertex_type curr_v = m_inverse_vertex_map[vertex_index];
768  //detect the odd cycles related, stored in the m_odd_cycles;
769  this->oddCycles(curr_v);
770 #ifdef DEBUG_LPCOLORING
771  this->write_graph_dot(curr_v);
772 #endif
773 
774  uint32_t odd_cycle_cnt = m_odd_cycles.size();
775  for(uint32_t m = 0; m < odd_cycle_cnt; m++)
776  {
777  uint32_t cycle_len = m_odd_cycles[m].size();
778  GRBLinExpr constraint1 = 0;
779  GRBLinExpr constraint2 = 0;
780 #ifdef DEBUG_LPCOLORING
781  assert(cycle_len%2 == 1);
782 #endif
783  for(uint32_t n = 0; n < cycle_len; ++n)
784  {
785 #ifdef DEBUG_LPCOLORING
786  assert(m_vertex_map.find(m_odd_cycles[m][n]) != m_vertex_map.end());
787 #endif
788  uint32_t vi_index = m_vertex_map[m_odd_cycles[m][n]];
789  constraint1 += coloringBits[vi_index<<1];
790  constraint2 += coloringBits[(vi_index<<1)+1];
791 
793  for(uint32_t x = 0; x < COLORING_BITS; x++)
794  {
795  non_integer_flag[vi_index*COLORING_BITS + x] = false;
796  }//end for x
797  }//end for
798  opt_model.addConstr(constraint1 >= 1);
799  opt_model.addConstr(constraint2 >= 1);
800  }//end for m
801  //for(uint32_t m = 0; m < COLORING_BITS; m++)
802  //{
803  // non_integer_flag[vertex_index*COLORING_BITS + m] = false;
804  //}
805  }//end for k
806 
807  //optimize the new model
808  opt_model.update();
809  cout << "================== LP iteration #" << m_lp_iter_cnt++ << " ==================\n";
810  opt_model.optimize();
811 #ifdef DEBUG_LPCOLORING
812  opt_model.write("graph.lp");
813  opt_model.write("graph.sol");
814 #endif
815  optim_status = opt_model.get(GRB_IntAttr_Status);
816  if(optim_status == GRB_INFEASIBLE)
817  {
818  cout << "ERROR: the model is infeasible... EXIT" << endl;
819  exit(1);
820  }
821  cout << endl;
822 
823  //no more non-integers are removed
824  uint32_t non_integer_num_updated, half_integer_num_updated;
825  pair = this->nonIntegerNum(coloringBits, vEdgeBit);
826  non_integer_num_updated = pair.first;
827  half_integer_num_updated = pair.second;
828 
829  if (/*non_integer_num_updated == 0*/ half_integer_num_updated == 0 ||
830  (non_integer_num_updated >= non_integer_num && half_integer_num_updated >= half_integer_num)) break;
831 
832  cout << "\n\n-----------------------------------Next iteration------------------------------" << endl;
833  }//end while
834 
835  //store the lp coloring results
836  this->store_lp_coloring(coloringBits);
837  /*
838  m_lp_coloring.clear();
839  m_lp_coloring.reserve(coloring_bits_num);
840  for(uint32_t k = 0; k < coloring_bits_num; k++)
841  {
842  double value = coloringBits[k].get(GRB_DoubleAttr_X);
843  m_lp_coloring.push_back(value);
844  }
845  */
846 
847 #ifdef DEBUG_LPCOLORING
848  // compare unrounded conflict number is more fair
849  this->lpConflictNum();
850  this->write_graph_dot();
851 #endif
852 
853  //rounding the coloring results
854  this->rounding_bindingAnalysis(opt_model, coloringBits, vEdgeBit);
855  if (roundingScheme() == DIRECT_ILP)
856  this->rounding_Greedy(coloringBits);
857  else if (roundingScheme() == FIXED_ILP)
858  this->ILPColoring(coloringBits, vEdgeBit);
859  else if (roundingScheme() == ITERATIVE_ILP)
860  this->rounding_ILP(opt_model, coloringBits, vEdgeBit);
861  else if (roundingScheme() == GREEDY)
862  this->rounding_Greedy(coloringBits);
863 
864  this->conflictNum();
865  this->stitchNum();
866  this->write_graph_color();
867 #ifdef DEBUG_LPCOLORING
868  //this->write_graph_dot();
869 #endif
870  switch (roundingScheme())
871  {
872  case DIRECT_ILP:
873  cout << "DIRECT_ILP mode "; break;
874  case FIXED_ILP:
875  cout << "FIXED_ILP mode "; break;
876  case ITERATIVE_ILP:
877  cout << "ITERATIVE_ILP mode "; break;
878  case GREEDY:
879  cout << "GREEDY mode "; break;
880  default:
881  cout << "Unknown mode "; assert(0);
882  }
883  cout << "problem solved in " << m_lp_iter_cnt << " LP iterations, " << m_ilp_iter_cnt << " ILP iterations" << endl;
884 
885 }//end coloring
886 
887 // ILP based coloring
888 // all integer bits are fixed
889 // non-integer bits are set to integer in the ILP problem
890 template<typename GraphType>
891 void LPColoring<GraphType>::ILPColoring(vector<GRBVar>& coloringBits, vector<GRBVar>& vEdgeBit)
892 {
893  uint32_t coloring_bits_num = coloringBits.size();
894 
895  //set up the LP environment
896  GRBEnv env = GRBEnv();
897  //mute the log from the LP solver
898  //env.set(GRB_IntParam_OutputFlag, 0);
899  GRBModel opt_model_updated = GRBModel(env);
900  //set up the LP variables
901  vector<GRBVar> coloringBits_updated;
902  vector<GRBVar> vEdgeBit_updated;
903  // vertex and edge variables
904  coloringBits_updated.reserve(coloring_bits_num);
905  //set the non-integer solution of the
906  for(uint32_t k = 0; k < coloring_bits_num; k++)
907  {
908  ostringstream oss;
909  oss << "v" << k;
910  double value = coloringBits[k].get(GRB_DoubleAttr_X);
911  if(value == 0.0)
912  {
913  //integer solution 0.0
914  coloringBits_updated.push_back(opt_model_updated.addVar(0.0, 0.0, 0.0, GRB_CONTINUOUS, oss.str()));
915  }
916  else if (value == 1.0)
917  {
918  //integer solution 1.0
919  coloringBits_updated.push_back(opt_model_updated.addVar(1.0, 1.0, 0.0, GRB_CONTINUOUS, oss.str()));
920  }
921  else
922  {
923  //non-integer solution
924  coloringBits_updated.push_back(opt_model_updated.addVar(0.0, 1.0, 0.0, GRB_INTEGER, oss.str()));
925  }
926  }//end for
927  uint32_t edge_num = vEdgeBit.size();
928  vEdgeBit_updated.reserve(edge_num);
929  for (uint32_t i = 0; i != edge_num; ++i)
930  {
931  // some variables here may not be used
932  ostringstream oss;
933  oss << "e" << i;
934  //vEdgeBit_updated.push_back(opt_model_updated.addVar(0.0, 1.0, 0.0, GRB_INTEGER));
935  vEdgeBit_updated.push_back(opt_model_updated.addVar(0.0, 1.0, 0.0, GRB_CONTINUOUS, oss.str()));
936  }
937 
938  opt_model_updated.update();
939  //iterate through the edges
940  pair<typename graph_type::edge_iterator, typename graph_type::edge_iterator> edge_range = edges(m_graph);
941  //set new objective
942  GRBLinExpr obj_init (0);
943  for (BOOST_AUTO(itr, edge_range.first); itr != edge_range.second; ++itr)
944  {
945  BOOST_AUTO(w, m_edge_weight_map[*itr]);
946  uint32_t edge_idx = m_edge_map[*itr];
947  if (w >= 0) // conflict
948  obj_init += vEdgeBit_updated[edge_idx];
949  else // stitch
950  obj_init += m_stitch_weight*vEdgeBit_updated[edge_idx];
951  }
952  opt_model_updated.setObjective(obj_init, GRB_MINIMIZE);
953 
954  for(BOOST_AUTO(itr, edge_range.first); itr != edge_range.second; ++itr)
955  {
956  BOOST_AUTO(from, source(*itr, m_graph));
957  uint32_t f_ind = m_vertex_map[from];
958  BOOST_AUTO(to, target(*itr, m_graph));
959  uint32_t t_ind = m_vertex_map[to];
960  //coloring conflict constraints
961  uint32_t vertex_idx1 = f_ind<<1;
962  uint32_t vertex_idx2 = t_ind<<1;
963 
964  BOOST_AUTO(w, m_edge_weight_map[*itr]);
965  if (w >= 0) // constraints for conflict edges
966  {
967  uint32_t edge_idx = m_edge_map[*itr];
968  opt_model_updated.addConstr(
969  coloringBits_updated[vertex_idx1] + coloringBits_updated[vertex_idx1+1]
970  + coloringBits_updated[vertex_idx2] + coloringBits_updated[vertex_idx2+1]
971  + vEdgeBit_updated[edge_idx] >= 1
972  );
973  opt_model_updated.addConstr(
974  1 - coloringBits_updated[vertex_idx1] + coloringBits_updated[vertex_idx1+1]
975  + 1 - coloringBits_updated[vertex_idx2] + coloringBits_updated[vertex_idx2+1]
976  + vEdgeBit_updated[edge_idx] >= 1
977  );
978  opt_model_updated.addConstr(
979  coloringBits_updated[vertex_idx1] + 1 - coloringBits_updated[vertex_idx1+1]
980  + coloringBits_updated[vertex_idx2] + 1 - coloringBits_updated[vertex_idx2+1]
981  + vEdgeBit_updated[edge_idx] >= 1
982  );
983  opt_model_updated.addConstr(
984  1 - coloringBits_updated[vertex_idx1] + 1 - coloringBits_updated[vertex_idx1+1]
985  + 1 - coloringBits_updated[vertex_idx2] + 1 - coloringBits_updated[vertex_idx2+1]
986  + vEdgeBit_updated[edge_idx] >= 1
987  );
988  }
989  else // constraints for stitch edges
990  {
991  uint32_t edge_idx = m_edge_map[*itr];
992  opt_model_updated.addConstr(
993  coloringBits_updated[vertex_idx1] - coloringBits_updated[vertex_idx2] - vEdgeBit_updated[edge_idx] <= 0
994  );
995  opt_model_updated.addConstr(
996  coloringBits_updated[vertex_idx2] - coloringBits_updated[vertex_idx1] - vEdgeBit_updated[edge_idx] <= 0
997  );
998  opt_model_updated.addConstr(
999  coloringBits_updated[vertex_idx1+1] - coloringBits_updated[vertex_idx2+1] - vEdgeBit_updated[edge_idx] <= 0
1000  );
1001  opt_model_updated.addConstr(
1002  coloringBits_updated[vertex_idx2+1] - coloringBits_updated[vertex_idx1+1] - vEdgeBit_updated[edge_idx] <= 0
1003  );
1004  }
1005  }//end for itr
1006  if (colorNum() == THREE)
1007  {
1008  for(uint32_t k = 0; k < coloring_bits_num; k += COLORING_BITS)
1009  {
1010  opt_model_updated.addConstr(coloringBits_updated[k] + coloringBits_updated[k+1] <= 1);
1011  }
1012  }
1013 
1014  opt_model_updated.update();
1015 
1016 #ifdef DEBUG_LPCOLORING
1017  opt_model_updated.write("graph_ilp.lp");
1018 #endif
1019 
1020  //optimize model
1021  cout << "================== ILP iteration #" << m_ilp_iter_cnt++ << " ==================\n";
1022  opt_model_updated.optimize();
1023  int optim_status = opt_model_updated.get(GRB_IntAttr_Status);
1024  if(optim_status == GRB_INFEASIBLE)
1025  {
1026  cout << "ERROR: The model is infeasible... EXIT" << endl;
1027  exit(1);
1028  }
1029  cout << endl;
1030  cout << "ILP solved to seek a coloring solution." << endl;
1031 
1032  //store the lp coloring results
1033  this->store_lp_coloring(coloringBits_updated);
1034  //store the coloring results
1035  for(uint32_t k = 0; k < coloring_bits_num; k = k + COLORING_BITS) {
1036  uint32_t vertex_index = k/COLORING_BITS;
1037  graph_vertex_type vertex_key = this->m_inverse_vertex_map[vertex_index];
1038  uint32_t color = 0;
1039  uint32_t base = 1;
1040  for(uint32_t m = 0; m < COLORING_BITS; m++) {
1041  double value = coloringBits_updated[k+m].get(GRB_DoubleAttr_X);
1042 #if DEBUG_LPCOLORING
1043  assert(value == 0.0 || value == 1.0);
1044 #endif
1045  if(value == 1.0) color = color + base;
1046  base = base<<1;
1047  }
1048  //set the coloring
1049  this->m_coloring[vertex_key] = color;
1050  }
1051 }
1052 
1054 template<typename GraphType>
1055 void LPColoring<GraphType>::rounding_bindingAnalysis(GRBModel& opt_model, vector<GRBVar>& coloringBits, vector<GRBVar>& vEdgeBit) {
1056  cout << "\n\n---------------------------------------Optimal rounding from binding analysis------------------------------------\n";
1057  /*
1058 //{{{
1059  GRBConstr* constrs_ptr = opt_model.getConstrs();
1060  int constrs_num = opt_model.get(GRB_IntAttr_NumConstrs);
1061  for(int k = 0; k < constrs_num; k++) {
1062  double slack = constrs_ptr[k].get(GRB_DoubleAttr_Slack);
1063  cout << "The slack for the " << k << "th constraint - " << constrs_ptr[k].get(GRB_StringAttr_ConstrName) << ": " << slack << endl;
1064  }//end for k
1065  */
1066  uint32_t non_integer_num, half_integer_num;
1067  uint32_t non_integer_num_updated, half_integer_num_updated;
1068 /*
1069  //single digit rounding scheme
1070  while(true) {
1071  BOOST_AUTO(pair, this->nonIntegerNum(coloringBits, vEdgeBit));
1072  non_integer_num = pair.first;
1073  half_integer_num = pair.second;
1074  if(non_integer_num == 0) return;
1075  //store the lp coloring results
1076  this->store_lp_coloring(coloringBits);
1077 
1078  //rounding scheme
1079  bool rounding_flag = true;
1080  uint32_t variable_num = coloringBits.size();
1081  for(uint32_t k = 0; k < variable_num; k++) {
1082  double value = coloringBits[k].get(GRB_DoubleAttr_X);
1083  if(isInteger(value)) continue;
1084 #ifdef DEBUG_LPCOLORING
1085  cout << endl << endl << k << "th non-integer variable" << endl;
1086  //string check_point;
1087  //cin >> check_point;
1088 #endif
1089  rounding_flag = true;
1090  char pre_sense = '=';
1091  bool pre_sense_flag = false;
1092  GRBColumn column = opt_model.getCol(coloringBits[k]);
1093  int column_size = column.size();
1094  for(int m = 0; m < column_size; m++) {
1095  GRBConstr constr = column.getConstr(m);
1096  double slack = constr.get(GRB_DoubleAttr_Slack);
1097  string name = constr.get(GRB_StringAttr_ConstrName);
1098  //ignore the non-binding constr
1099  if(slack != 0.0) continue;
1100  //ignore the stitch constraint
1101  if(true == m_stitch_constrs[name]) continue;
1102  double coeff = opt_model.getCoeff(constr, coloringBits[k]);
1103  char sense = constr.get(GRB_CharAttr_Sense);
1104 #ifdef DEBUG_LPCOLORING
1105  cout << "The slack for " << constr.get(GRB_StringAttr_ConstrName) << ": " << slack << endl;
1106  cout << "The constraint sense for " << constr.get(GRB_StringAttr_ConstrName) << ": " << sense << endl;
1107  cout << "The coefficient for " << constr.get(GRB_StringAttr_ConstrName) << ": " << coeff << endl << endl;
1108 #endif
1109 #ifdef DEBUG_LPCOLORING
1110  assert(coeff == 1.0 || coeff == -1.0);
1111  assert(sense == '>' || sense == '<' || sense == '=');
1112 #endif
1113  //flip the sense based on coeff
1114  if(coeff == -1.0) {
1115  if(sense == '>') sense = '<';
1116  else if(sense == '<') sense = '>';
1117  }//end if
1118 
1119  //set the rounding_flag
1120  if(false == pre_sense_flag) {
1121  pre_sense = sense;
1122  pre_sense_flag = true;
1123  } else {
1124  if(sense != pre_sense) {
1125  rounding_flag = false;
1126  //break;
1127  }
1128  }//end if
1129  }//end for m
1130 
1131  //rounding the first feasible variable
1132  if(true == rounding_flag) {
1133 #ifdef DEBUG_LPCOLORING
1134  assert(pre_sense == '<' || pre_sense == '>');
1135 #endif
1136  if(pre_sense == '<') coloringBits[k].set(GRB_DoubleAttr_UB, 0.0);
1137  else coloringBits[k].set(GRB_DoubleAttr_LB, 1.0);
1138  break;
1139  }
1140  }//end for k
1141 
1142  if(true == rounding_flag) {
1143  //update the model and optimize
1144  opt_model.update();
1145  opt_model.optimize();
1146 #ifdef DEBUG_LPCOLORING
1147  opt_model.write("graph.lp");
1148  opt_model.write("graph.sol");
1149 #endif
1150  int optim_status = opt_model.get(GRB_IntAttr_Status);
1151  if(optim_status == GRB_INFEASIBLE)
1152  {
1153  cout << "ERROR: the ILP model is infeasible... Quit ILP based rounding" << endl;
1154  exit(1);
1155  }
1156  }//end if
1157 
1158  //no more non-integers are removed
1159  pair = this->nonIntegerNum(coloringBits, vEdgeBit);
1160  non_integer_num_updated = pair.first;
1161  half_integer_num_updated = pair.second;
1162  if (half_integer_num_updated == 0 ||
1163  (non_integer_num_updated >= non_integer_num && half_integer_num_updated >= half_integer_num)) break;
1164 
1165  }//end while
1166 //}}}
1167 */
1168 
1169  //a pair of digits rounding scheme, essentially this is a hard coding implementation
1170  while(true) {
1171  BOOST_AUTO(pair, this->nonIntegerNum(coloringBits, vEdgeBit));
1172  non_integer_num = pair.first;
1173  half_integer_num = pair.second;
1174  if(non_integer_num == 0) return;
1175  //store the lp coloring results
1176  this->store_lp_coloring(coloringBits);
1177 
1178  //rounding scheme
1179  bool rounding_flag = true;
1180  double bit1 = 0.0, bit2 = 0.0;
1181  uint32_t variable_num = coloringBits.size();
1182  for(uint32_t k = 0; k < variable_num; k = k + COLORING_BITS) {
1183  double value1 = coloringBits[k].get(GRB_DoubleAttr_X);
1184  double value2 = coloringBits[k+1].get(GRB_DoubleAttr_X);
1185  if(isInteger(value1) && isInteger(value2)) continue;
1186 #ifdef DEBUG_LPCOLORING
1187  //cout << endl << endl << k << "th and " << (k + 1) << "th non-integer variable" << endl;
1188  //string check_point;
1189  //cin >> check_point;
1190 #endif
1191  GRBColumn column1 = opt_model.getCol(coloringBits[k]);
1192  GRBColumn column2 = opt_model.getCol(coloringBits[k+1]);
1193  int column1_size = column1.size();
1194  int column2_size = column2.size();
1195  //constraints flag, shared constraints betweeen column1 & 2 are flagged true, otherwise false
1196  //coefficient, -1 or 1
1197  //sense, < or >
1198  unordered_map<string, bool> constrs_flag1;
1199  unordered_map<string, double> constrs_coeff1;
1200  unordered_map<string, char> constrs_sense1;
1201  unordered_map<string, bool> constrs_flag2;
1202  unordered_map<string, double> constrs_coeff2;
1203  unordered_map<string, char> constrs_sense2;
1204  //process the first set of binding constraints
1205 #ifdef DEBUG_LPCOLORING
1206  cout << endl << endl << k << "th non-integer variable" << endl;
1207 #endif
1208  for(int m = 0; m < column1_size; m++) {
1209  GRBConstr constr = column1.getConstr(m);
1210  double slack = constr.get(GRB_DoubleAttr_Slack);
1211  //ignore the non-binding constr
1212  if(slack != 0.0) continue;
1213  double coeff = opt_model.getCoeff(constr, coloringBits[k]);
1214  char sense = constr.get(GRB_CharAttr_Sense);
1215 #ifdef DEBUG_LPCOLORING
1216  cout << "The slack for " << constr.get(GRB_StringAttr_ConstrName) << ": " << slack << endl;
1217  cout << "The constraint sense for " << constr.get(GRB_StringAttr_ConstrName) << ": " << sense << endl;
1218  cout << "The coefficient for " << constr.get(GRB_StringAttr_ConstrName) << ": " << coeff << endl << endl;
1219 #endif
1220 #ifdef DEBUG_LPCOLORING
1221  assert(coeff == 1.0 || coeff == -1.0);
1222  assert(sense == '>' || sense == '<');
1223 #endif
1224  string name = constr.get(GRB_StringAttr_ConstrName);
1225  constrs_flag1[name] = false;
1226  constrs_coeff1[name] = coeff;
1227  constrs_sense1[name] = sense;
1228  }//end for m
1229 
1230 #ifdef DEBUG_LPCOLORING
1231  cout << endl << endl << (k+1) << "th non-integer variable" << endl;
1232 #endif
1233  //process the second set of binding constraints
1234  for(int m = 0; m < column2_size; m++) {
1235  GRBConstr constr = column2.getConstr(m);
1236  double slack = constr.get(GRB_DoubleAttr_Slack);
1237  //ignore the non-binding constr
1238  if(slack != 0.0) continue;
1239  double coeff = opt_model.getCoeff(constr, coloringBits[k+1]);
1240  char sense = constr.get(GRB_CharAttr_Sense);
1241 #ifdef DEBUG_LPCOLORING
1242  cout << "The slack for " << constr.get(GRB_StringAttr_ConstrName) << ": " << slack << endl;
1243  cout << "The constraint sense for " << constr.get(GRB_StringAttr_ConstrName) << ": " << sense << endl;
1244  cout << "The coefficient for " << constr.get(GRB_StringAttr_ConstrName) << ": " << coeff << endl << endl;
1245 #endif
1246 #ifdef DEBUG_LPCOLORING
1247  assert(coeff == 1.0 || coeff == -1.0);
1248  assert(sense == '>' || sense == '<');
1249 #endif
1250  string name = constr.get(GRB_StringAttr_ConstrName);
1251  if(constrs_flag1.find(name) == constrs_flag1.end())
1252  {
1253  constrs_flag2[name] = false;
1254  constrs_coeff2[name] = coeff;
1255  constrs_sense2[name] = sense;
1256  }
1257  else
1258  {
1259  constrs_flag1[name] = true;
1260  constrs_flag2[name] = true;
1261  constrs_coeff2[name] = coeff;
1262  constrs_sense2[name] = sense;
1263  }//end if
1264  }//end for m
1265 
1266  //flip the sense of the un-shared constraints if necessary
1267  unordered_map<string, double> shared_constrs_coeff1;
1268  unordered_map<string, char> shared_constrs_sense1;
1269  BOOST_AUTO(end, constrs_flag1.end());
1270  for(BOOST_AUTO(itr, constrs_flag1.begin()); itr != end; ++itr) {
1271  //ignore shared constraint
1272  if(true == itr->second)
1273  {
1274  shared_constrs_coeff1[itr->first] = constrs_coeff1[itr->first];
1275  shared_constrs_sense1[itr->first] = constrs_sense1[itr->first];
1276  continue;
1277  }
1278  if(-1.0 == constrs_coeff1[itr->first])
1279  {
1280  //flip sense
1281  constrs_coeff1[itr->first] = 1.0;
1282  if(constrs_sense1[itr->first] == '<') constrs_sense1[itr->first] = '>';
1283  else constrs_sense1[itr->first] = '<';
1284  }//end if
1285  }
1286 
1287  unordered_map<string, double> shared_constrs_coeff2;
1288  unordered_map<string, double> shared_constrs_sense2;
1289  end = constrs_flag2.end();
1290  for(BOOST_AUTO(itr, constrs_flag2.begin()); itr != end; ++itr) {
1291  //ignore shared constraint
1292  if(true == itr->second)
1293  {
1294  shared_constrs_coeff2[itr->first] = constrs_coeff2[itr->first];
1295  shared_constrs_sense2[itr->first] = constrs_sense2[itr->first];
1296  continue;
1297  }
1298  if(-1.0 == constrs_coeff2[itr->first])
1299  {
1300  //flip sense
1301  constrs_coeff2[itr->first] = 1.0;
1302  if(constrs_sense2[itr->first] == '<') constrs_sense2[itr->first] = '>';
1303  else constrs_sense2[itr->first] = '<';
1304  }//end if
1305  }
1306 #ifdef DEBUG_LPCOLORING
1307  assert(shared_constrs_coeff1.size() == shared_constrs_coeff2.size());
1308 #endif
1309 
1310  //check whether value1 can be rounded
1311  bool rounding_flag1 = true;
1312  char pre_sense1 = '=';
1313  bool pre_sense_flag = false;
1314  end = constrs_flag1.end();
1315  for(BOOST_AUTO(itr, constrs_flag1.begin()); itr != end; ++itr) {
1316  //ignore shared constraint
1317  if(true == itr->second) continue;
1318 #ifdef DEBUG_LPCOLORING
1319  assert(1.0 == constrs_coeff1[itr->first]);
1320 #endif
1321  if(false == pre_sense_flag)
1322  {
1323  pre_sense1 = constrs_sense1[itr->first];
1324  pre_sense_flag = true;
1325  }
1326  else
1327  {
1328  if(pre_sense1 != constrs_sense1[itr->first]) {
1329  rounding_flag1 = false;
1330  break;
1331  }
1332  }//end if
1333  }//end for itr
1334 
1335  //check whether value2 can be rounded
1336  bool rounding_flag2 = true;
1337  char pre_sense2 = '=';
1338  pre_sense_flag = false;
1339  end = constrs_flag2.end();
1340  for(BOOST_AUTO(itr, constrs_flag2.begin()); itr != end; ++itr) {
1341  //ignore the shared constraint
1342  if(true == itr->second) continue;
1343 #ifdef DEBUG_LPCOLORING
1344  assert(1.0 == constrs_coeff2[itr->first]);
1345 #endif
1346  if(false == pre_sense_flag)
1347  {
1348  pre_sense2 = constrs_sense2[itr->first];
1349  pre_sense_flag = true;
1350  }
1351  else
1352  {
1353  if(pre_sense2 != constrs_sense2[itr->first]) {
1354  rounding_flag2 = false;
1355  break;
1356  }
1357  }//end if
1358  }//end for itr
1359 
1360  //check value1 and value2 can be rounded on un-shared constraints
1361  rounding_flag = rounding_flag1 & rounding_flag2;
1362  if(false == rounding_flag) break;
1363 
1364  //check whether value1&2 can be rounded simultaneously for shared constraints
1365  //flag whether a feasible flag is found or not
1366  bool bit_pair_flag = false;
1367  bit1 = 0.0, bit2 = 0.0;
1368  for(bit1 = 0.0; bit1 <= 1.0; bit1 = bit1 + 1.0) {
1369  if(pre_sense1 == '<' && bit1 == 1.0) continue;
1370  if(pre_sense1 == '>' && bit1 == 0.0) continue;
1371  for(bit2 = 0.0; bit2 <= 1.0; bit2 = bit2 + 1.0) {
1372  if(pre_sense2 == '<' && bit2 == 1.0) continue;
1373  if(pre_sense2 == '>' && bit2 == 0.0) continue;
1374  cout << "bit1: " << bit1 << "; bit2: " << bit2 << endl;
1375  if (colorNum() == THREE && (bit1 == 1.0) && (bit2 == 1.0)) continue;
1376  //check each shared constraint
1377  BOOST_AUTO(itr_end, shared_constrs_coeff1.end());
1378  for(BOOST_AUTO(itr, shared_constrs_coeff1.begin()); itr != itr_end; ++itr) {
1379  double delta_value = itr->second * (bit1 - value1) + shared_constrs_coeff2[itr->first] * (bit2 - value2);
1380  char sense = shared_constrs_sense1[itr->first];
1381  if(delta_value < 0 && sense == '>') {
1382  rounding_flag = false;
1383  break;
1384  }
1385  if(delta_value > 0 && sense == '<') {
1386  rounding_flag = false;
1387  break;
1388  }
1389  }//end for itr
1390  bit_pair_flag = true;
1391  if(true == rounding_flag) break;
1392  }//end for bit2
1393  if(true == rounding_flag) break;
1394  }//end for bit1
1395  //rounding the first feasible variable
1396  if(true == rounding_flag && true == bit_pair_flag) {
1397 #ifdef DEBUG_LPCOLORING
1398 #endif
1399  coloringBits[k].set(GRB_DoubleAttr_UB, bit1);
1400  coloringBits[k].set(GRB_DoubleAttr_LB, bit1);
1401  coloringBits[k+1].set(GRB_DoubleAttr_UB, bit2);
1402  coloringBits[k+1].set(GRB_DoubleAttr_LB, bit2);
1403  break;
1404  }
1405  }//end for k
1406 
1407  if(true == rounding_flag) {
1408  //update the model and optimize
1409  opt_model.update();
1410  opt_model.optimize();
1411 #ifdef DEBUG_LPCOLORING
1412  opt_model.write("graph.lp");
1413  opt_model.write("graph.sol");
1414  this->store_lp_coloring(coloringBits);
1415  this->write_graph_dot();
1416 #endif
1417  int optim_status = opt_model.get(GRB_IntAttr_Status);
1418  if(optim_status == GRB_INFEASIBLE)
1419  {
1420  cout << "ERROR: the model after binding analysis is infeasible... Quit ILP based rounding" << endl;
1421  exit(1);
1422  }
1423  }//end if
1424 
1425  //no more non-integers are removed
1426  pair = this->nonIntegerNum(coloringBits, vEdgeBit);
1427  non_integer_num_updated = pair.first;
1428  half_integer_num_updated = pair.second;
1429  if (/*non_integer_num_updated == 0*/ half_integer_num_updated == 0 ||
1430  (non_integer_num_updated >= non_integer_num && half_integer_num_updated >= half_integer_num)) break;
1431 
1432  }//end while
1433 
1434 
1435  this->store_lp_coloring(coloringBits);
1436  return;
1437 }
1438 
1440 template<typename GraphType>
1441 void LPColoring<GraphType>::rounding_ILP(GRBModel& opt_model, vector<GRBVar>& coloringBits, vector<GRBVar>& vEdgeBit)
1442 {
1443  cout << "\n\n---------------------------------------ILP rounding scheme------------------------------------\n";
1444  uint32_t non_integer_num, half_integer_num;
1445  uint32_t non_integer_num_updated, half_integer_num_updated;
1446  while(true)
1447  {
1448  BOOST_AUTO(pair, this->nonIntegerNum(coloringBits, vEdgeBit));
1449  non_integer_num = pair.first;
1450  half_integer_num = pair.second;
1451  if(non_integer_num == 0) break;
1452  //store the lp coloring results
1453  this->store_lp_coloring(coloringBits);
1454 #ifdef DEBUG_LPCOLORING
1455  // compare unrounded conflict number is more fair
1456  this->lpConflictNum();
1457 #endif
1458  //update the ILP model
1459  uint32_t coloring_bits_num = coloringBits.size();
1460  for(uint32_t k = 0; k < coloring_bits_num; ++k)
1461  {
1462  double value = coloringBits[k].get(GRB_DoubleAttr_X);
1463  if(value == 1.0 || value == 0.0)
1464  {
1465  coloringBits[k].set(GRB_CharAttr_VType, 'C');
1466  }
1467  else
1468  {
1469  coloringBits[k].set(GRB_CharAttr_VType, 'I');
1470  }
1471  }//end for k
1472  opt_model.update();
1473 #ifdef DEBUG_LPCOLORING
1474  opt_model.write("graph_ilp.lp");
1475 #endif
1476  cout << "================== ILP iteration #" << m_ilp_iter_cnt++ << " ==================\n";
1477  opt_model.optimize();
1478  int optim_status = opt_model.get(GRB_IntAttr_Status);
1479  if(optim_status == GRB_INFEASIBLE)
1480  {
1481  cout << "ERROR: the ILP model is infeasible... Quit ILP based rounding" << endl;
1482 #ifdef DEBUG_LPCOLORING
1483  opt_model.computeIIS();
1484  opt_model.write("graph_ilp.ilp");
1485 #endif
1486  exit(1);
1487  }
1488 
1489  //no more non-integers are removed
1490  pair = this->nonIntegerNum(coloringBits, vEdgeBit);
1491  non_integer_num_updated = pair.first;
1492  half_integer_num_updated = pair.second;
1493  if (/*non_integer_num_updated == 0*/ half_integer_num_updated == 0 ||
1494  (non_integer_num_updated >= non_integer_num && half_integer_num_updated >= half_integer_num)) break;
1495  }//end while
1496 
1497  //final greedy rounding
1498  this->rounding_Greedy(coloringBits);
1499 }
1500 
1501 //greedy rounding scheme, need better scheme later on
1502 template<typename GraphType>
1503 void LPColoring<GraphType>::rounding_Greedy_v0(vector<GRBVar>& coloringBits)
1504 {
1505  cout << "\n\n---------------------------------------Greedy rounding scheme------------------------------------\n";
1506 
1507  //greedy rounding scheme
1508  uint32_t vec_size = coloringBits.size();
1509  //the rounding results
1510  vector<bool> coloringBinary(vec_size, false);
1511  //the flag denoting whether current bit is rounded or not
1512  vector<bool> roundingFlag(vec_size, false);
1513  //rounding by range
1514  for(uint32_t k = 0; k < vec_size; k++)
1515  {
1516  double value = coloringBits[k].get(GRB_DoubleAttr_X);
1517  if (0.0 == value)
1518  {
1519  coloringBinary[k] = false;
1520  roundingFlag[k] = true;
1521  }
1522  else if (1.0 == value)
1523  {
1524  coloringBinary[k] = true;
1525  roundingFlag[k] = true;
1526  }
1527  else
1528  {
1529  coloringBinary[k] = false;
1530  roundingFlag[k] = false;
1531  }//end if
1532  }//end for
1533 
1534  //rounding of the half integer
1535  //greedy rounding schme should minimize the conflict and stitch number
1536  for(uint32_t k = 0; k < vec_size; k++)
1537  {
1538  if(true == roundingFlag[k]) continue;
1539 #ifdef DEBUG_LPCOLORING
1540  //this->printBoolArray(roundingFlag, vec_size);
1541  //this->printBoolArray(coloringBinary, vec_size);
1542  //cout << endl << endl;
1543 #endif
1544  //greedy rounding scheme
1545  uint32_t vertex_index = k/COLORING_BITS;
1546  vector<bool> color_bits;
1547  color_bits.reserve(COLORING_BITS);
1548  vector<bool> best_bits;
1549  best_bits.reserve(COLORING_BITS);
1550  //initialize the color_bits
1551  for(uint32_t m = 0; m < COLORING_BITS; ++m)
1552  {
1553  color_bits.push_back(coloringBinary[vertex_index*COLORING_BITS + m]);
1554  }//end for
1555 
1556  //get the neighbors
1557  typename boost::graph_traits<graph_type>::adjacency_iterator vi_begin, vi_end;
1558  BOOST_AUTO(vertex_key, m_inverse_vertex_map[vertex_index]);
1559  boost::tie(vi_begin, vi_end) = adjacent_vertices(vertex_key, m_graph);
1560  //calculate the current
1561  uint32_t same_color_bound = std::numeric_limits<uint32_t>::max();
1562  uint32_t same_color_count = 0;
1563  uint32_t color = 0;
1564  uint32_t base = 1;
1565  while(true)
1566  {
1567  same_color_bound = std::numeric_limits<uint32_t>::max();
1568  same_color_count = 0;
1569  color = 0;
1570  base = 1;
1571  for(uint32_t m = 0; m < COLORING_BITS; ++m)
1572  {
1573  if(color_bits[m]) color = color + base;
1574  base = base<<1;
1575  }//end for
1576  //check the same color neighbors
1577  for(BOOST_AUTO(vi, vi_begin); vi != vi_end; ++vi)
1578  {
1579  if(this->m_coloring.find(*vi) == m_vertex_map.end()) continue;
1580  if(this->m_coloring[*vi] == color) same_color_count++;
1581  }//end for
1582  //assign better color
1583  if(same_color_count < same_color_bound)
1584  {
1585  same_color_bound = same_color_count;
1586  best_bits = color_bits;
1587  }
1588 
1589  //explore the next color_bits
1590  bool nextFlag = false;
1591  for(uint32_t m = 0; m < COLORING_BITS; ++m)
1592  {
1593  if(color_bits[m] == false && roundingFlag[vertex_index*COLORING_BITS + m] == false)
1594  {
1595  //flip the color bit that has not be rounded
1596  color_bits[m] = true;
1597  nextFlag = true;
1598  //for Triple patterning mode
1599  if((colorNum() == THREE) && (m == COLORING_BITS - 1)) nextFlag = false;
1600  break;
1601  }
1602  }//end for m
1603  //all the color_bits are explored
1604  if(nextFlag == false) break;
1605  }//end while
1606 
1607  //the vertex is colored
1608  color = 0;
1609  base = 1;
1610  for(uint32_t m = 0; m < COLORING_BITS; m++)
1611  {
1612  coloringBinary[vertex_index*COLORING_BITS + m] = best_bits[m];
1613  roundingFlag[vertex_index*COLORING_BITS + m] = true;
1614  if(best_bits[m]) color = color + base;
1615  base = base<<1;
1616  }
1617 #ifdef DEBUG_LPCOLORING
1618  assert(color < 4);
1619 #endif
1620  this->m_coloring[vertex_key] = color;
1621  }//end for k
1622 
1623  //assign the coloring results
1624  for(uint32_t k = 0; k < vec_size; k = k + COLORING_BITS)
1625  {
1626  BOOST_AUTO(vertex_key, this->m_inverse_vertex_map[(k/COLORING_BITS)]);
1627  uint32_t color = 0;
1628  uint32_t base = 1;
1629  for(uint32_t m = 0; m < COLORING_BITS; ++m)
1630  {
1631  if(coloringBinary[k + m]) {
1632  this->m_lp_coloring[k + m] = 1.0;
1633  color = color + base;
1634  } else {
1635  this->m_lp_coloring[k + m] = 0.0;
1636  }
1637  base = base<<1;
1638  }//end for k
1639  if(this->m_coloring.find(vertex_key) == this->m_coloring.end())
1640  this->m_coloring[vertex_key] = color;
1641 #ifdef DEBUG_LPCOLORING
1642  else
1643  {
1644  //cout << "stored color-" << this->m_coloring[vertex_key] << " vs assigned color-" << color << endl;
1645  assert(this->m_coloring[vertex_key] == color);
1646  }
1647 #endif
1648  }//end for
1649 }
1650 
1651 template<typename GraphType>
1652 void LPColoring<GraphType>::rounding_Greedy(vector<GRBVar>& coloringBits)
1653 {
1654  cout << "\n\n---------------------------------------Greedy rounding scheme------------------------------------\n";
1655 
1656  //greedy rounding scheme
1657  uint32_t vec_size = coloringBits.size();
1658  this->store_lp_coloring(coloringBits);
1659 
1660  //greedy rounding schme should minimize the conflict and stitch number
1661  for(uint32_t k = 0; k < vec_size; k = k + COLORING_BITS)
1662  {
1663  double value1 = coloringBits[k].get(GRB_DoubleAttr_X);
1664  double value2 = coloringBits[k+1].get(GRB_DoubleAttr_X);
1665  if(isInteger(value1) && isInteger(value2)) continue;
1666  //greedy rounding scheme
1667  uint32_t vertex_index = k/COLORING_BITS;
1668 
1669  //get the neighbors
1670  typename boost::graph_traits<graph_type>::adjacency_iterator vi_begin, vi_end;
1671  BOOST_AUTO(vertex_key, m_inverse_vertex_map[vertex_index]);
1672  boost::tie(vi_begin, vi_end) = adjacent_vertices(vertex_key, m_graph);
1673  //calculate the current
1674  uint32_t same_color_bound = std::numeric_limits<uint32_t>::max();
1675  uint32_t same_color_count = 0;
1676  double best_bit1 = 0.0, best_bit2 = 0.0;
1677  //greedy selection
1678  for(double bit1 = 0.0; bit1 <= 1.0; bit1 = bit1 + 1.0)
1679  {
1680  if(isInteger(value1) && (bit1 != value1)) continue;
1681  for(double bit2 = 0.0; bit2 <= 1.0; bit2 = bit2 + 1.0)
1682  {
1683  if(isInteger(value2) && (bit2 != value2)) continue;
1684  if (colorNum() == THREE && (bit1 == 1.0) && (bit2 == 1.0)) continue;
1685  same_color_count = 0;
1686  //cout << endl << "current_index: " << vertex_index << ", bits: " << bit1 << bit2 << endl;
1687  //check the same color neighbors
1688  for(BOOST_AUTO(vi, vi_begin); vi != vi_end; ++vi)
1689  {
1690 #ifdef DEBUG_LPCOLORING
1691  assert(this->m_vertex_map.find(*vi) != m_vertex_map.end());
1692 #endif
1693  uint32_t neighbor_index = this->m_vertex_map[*vi];
1694  //same color neighbor
1695  if((m_lp_coloring[2*neighbor_index] == bit1) && (m_lp_coloring[2*neighbor_index+1] == bit2))
1696  {
1697  same_color_count++;
1698  //cout << "conflict neighbor_index: " << neighbor_index << endl;
1699  }
1700  }//end for
1701  //assign better color
1702  if(same_color_count < same_color_bound)
1703  {
1704  same_color_bound = same_color_count;
1705  best_bit1 = bit1;
1706  best_bit2 = bit2;
1707  }
1708  }//end for bit2
1709  }//end for bit1
1710 
1711  //the vertex is colored
1712  this->m_lp_coloring[k] = best_bit1;
1713  this->m_lp_coloring[k+1] = best_bit2;
1714  }//end for k
1715 
1716  //assign the coloring results
1717  for(uint32_t k = 0; k < vec_size; k = k + COLORING_BITS)
1718  {
1719  BOOST_AUTO(vertex_key, this->m_inverse_vertex_map[(k/COLORING_BITS)]);
1720  uint32_t color = 0;
1721  uint32_t base = 1;
1722  for(uint32_t m = 0; m < COLORING_BITS; ++m)
1723  {
1724 #ifdef DEBUG_LPCOLORING
1725  assert(isInteger(this->m_lp_coloring[k + m]));
1726 #endif
1727  if(this->m_lp_coloring[k + m] == 1.0) color = color + base;
1728  base = base<<1;
1729  }//end for k
1730  if(this->m_coloring.find(vertex_key) == this->m_coloring.end())
1731  this->m_coloring[vertex_key] = color;
1732 #ifdef DEBUG_LPCOLORING
1733  else
1734  {
1735  //cout << "stored color-" << this->m_coloring[vertex_key] << " vs assigned color-" << color << endl;
1736  assert(this->m_coloring[vertex_key] == color);
1737  }
1738 #endif
1739  }//end for
1740 }
1741 
1742 
1743 
1744 //get the vertex color
1745 template<typename GraphType>
1746 uint32_t LPColoring<GraphType>::vertexColor(graph_vertex_type& node) const
1747 {
1748  //the graph is not colored
1749  if(this->m_coloring.empty()) return 0; //this->graphColoring(m_graph);
1750  return this->m_coloring.at(node);
1751 }
1752 
1753 template <typename GraphType>
1754 bool LPColoring<GraphType>::sameColor(typename LPColoring<GraphType>::graph_vertex_type v1,
1755  typename LPColoring<GraphType>::graph_vertex_type v2) const
1756 {
1757  assert(!m_coloring.empty());
1758 #ifdef DEBUG_LPCOLORING
1759  assert(this->m_coloring.find(v1) != this->m_coloring.end());
1760  assert(this->m_coloring.find(v2) != this->m_coloring.end());
1761 #endif
1762  return m_coloring.at(v1) == m_coloring.at(v2);
1763 }
1764 
1765 template <typename GraphType>
1766 bool LPColoring<GraphType>::hasConflict(typename LPColoring<GraphType>::graph_edge_type e) const
1767 {
1768  BOOST_AUTO(w, m_edge_weight_map[e]);
1769  // skip stitch edges
1770  if (w < 0) return false;
1771 
1772  BOOST_AUTO(from, source(e, m_graph));
1773  BOOST_AUTO(to, target(e, m_graph));
1774  return sameColor(from, to);
1775 }
1776 
1777 template <typename GraphType>
1778 bool LPColoring<GraphType>::hasStitch(typename LPColoring<GraphType>::graph_edge_type e) const
1779 {
1780  BOOST_AUTO(w, m_edge_weight_map[e]);
1781  // skip conflict edges
1782  if (w >= 0) return false;
1783 
1784  BOOST_AUTO(from, source(e, m_graph));
1785  BOOST_AUTO(to, target(e, m_graph));
1786  return !sameColor(from, to);
1787 }
1788 
1789 //report the conflict number
1790 template<typename GraphType>
1791 uint32_t LPColoring<GraphType>::conflictNum() const
1792 {
1793  //the graph is not colored
1794  if (this->m_coloring.empty()) return 0; //this->graphColoring();
1795  //check the coloring results
1796  uint32_t conflict_edge_num = 0;
1797  uint32_t conflict_num = 0;
1798  pair<typename graph_type::edge_iterator, typename graph_type::edge_iterator> edge_range = edges(m_graph);
1799  for(BOOST_AUTO(itr, edge_range.first); itr != edge_range.second; ++itr)
1800  {
1801  if (hasConflict(*itr))
1802  ++conflict_num;
1803  ++conflict_edge_num;
1804  }//end for itr
1805  cout << "Conflict number: " << conflict_num << " out of " << conflict_edge_num << " conflict edges" << endl;
1806  return conflict_num;
1807 }
1808 
1809 //report the stitch number
1810 template<typename GraphType>
1811 uint32_t LPColoring<GraphType>::stitchNum() const
1812 {
1813  //the graph is not colored
1814  if(this->m_coloring.empty()) return 0; //this->graphColoring(m_graph);
1815 
1816  //check the coloring results
1817  uint32_t stitch_edge_num = 0;
1818  uint32_t stitch_num = 0;
1819  pair<typename graph_type::edge_iterator, typename graph_type::edge_iterator> edge_range = edges(m_graph);
1820  for(BOOST_AUTO(itr, edge_range.first); itr != edge_range.second; ++itr)
1821  {
1822  if (hasStitch(*itr))
1823  ++stitch_num;
1824  ++stitch_edge_num;
1825  }//end for itr
1826  cout << "Stitch number: " << stitch_num << " out of " << stitch_edge_num << " stitch edges" << endl;
1827  return stitch_num;
1828 }
1829 
1830 template <typename GraphType>
1831 bool LPColoring<GraphType>::lpSameColor(typename LPColoring<GraphType>::graph_vertex_type v1,
1832  typename LPColoring<GraphType>::graph_vertex_type v2) const
1833 {
1834  assert(!m_lp_coloring.empty());
1835 
1836  uint32_t f_ind = m_vertex_map.at(v1);
1837  uint32_t t_ind = m_vertex_map.at(v2);
1838  uint32_t vertex_idx1 = f_ind<<1;
1839  uint32_t vertex_idx2 = t_ind<<1;
1840 
1841  return m_lp_coloring.at(vertex_idx1) == m_lp_coloring.at(vertex_idx2)
1842  && m_lp_coloring.at(vertex_idx1+1) == m_lp_coloring.at(vertex_idx2+1);
1843 }
1844 
1845 template <typename GraphType>
1846 bool LPColoring<GraphType>::hasLPConflict(typename LPColoring<GraphType>::graph_edge_type e) const
1847 {
1848  BOOST_AUTO(w, m_edge_weight_map[e]);
1849  // skip stitch edges
1850  if (w < 0) return false;
1851 
1852  BOOST_AUTO(from, source(e, m_graph));
1853  BOOST_AUTO(to, target(e, m_graph));
1854  return lpSameColor(from, to);
1855 }
1856 
1857 template <typename GraphType>
1858 bool LPColoring<GraphType>::hasLPStitch(typename LPColoring<GraphType>::graph_edge_type e) const
1859 {
1860  BOOST_AUTO(w, m_edge_weight_map[e]);
1861  // skip conflict edges
1862  if (w >= 0) return false;
1863 
1864  BOOST_AUTO(from, source(e, m_graph));
1865  BOOST_AUTO(to, target(e, m_graph));
1866  return !lpSameColor(from, to);
1867 }
1868 
1869 //report the unrounded LP conflict number
1870 template<typename GraphType>
1871 uint32_t LPColoring<GraphType>::lpConflictNum() const
1872 {
1873  if (m_lp_coloring.empty()) return 0;
1874  //check the coloring results
1875  uint32_t conflict_edge_num = 0;
1876  uint32_t conflict_num = 0;
1877  pair<typename graph_type::edge_iterator, typename graph_type::edge_iterator> edge_range = edges(m_graph);
1878  for(BOOST_AUTO(itr, edge_range.first); itr != edge_range.second; ++itr)
1879  {
1880  if (hasLPConflict(*itr))
1881  ++conflict_num;
1882  ++conflict_edge_num;
1883  }//end for itr
1884  cout << "LP Conflict number: " << conflict_num << " out of " << conflict_edge_num << " conflict edges" << endl;
1885  return conflict_num;
1886 }
1887 
1888 //report the unrounded LP stitch number
1889 template<typename GraphType>
1890 uint32_t LPColoring<GraphType>::lpStitchNum() const
1891 {
1892  if (m_lp_coloring.empty()) return 0;
1893  //check the coloring results
1894  uint32_t stitch_edge_num = 0;
1895  uint32_t stitch_num = 0;
1896  pair<typename graph_type::edge_iterator, typename graph_type::edge_iterator> edge_range = edges(m_graph);
1897  for(BOOST_AUTO(itr, edge_range.first); itr != edge_range.second; ++itr)
1898  {
1899  if (hasLPStitch(*itr))
1900  {
1901  ++stitch_num;
1902  }
1903  ++stitch_edge_num;
1904  }//end for itr
1905  cout << "LP Stitch number: " << stitch_num << " out of " << stitch_edge_num << " stitch edges" << endl;
1906  return stitch_num;
1907 }
1908 
1909 //store the lp coloring results
1910 template<typename GraphType>
1911 void LPColoring<GraphType>::store_lp_coloring(vector<GRBVar>& coloringBits) {
1912  uint32_t coloring_bits_num = coloringBits.size();
1913  m_lp_coloring.clear();
1914  m_lp_coloring.reserve(coloring_bits_num);
1915  for(uint32_t k = 0; k < coloring_bits_num; k++)
1916  {
1917  double value = coloringBits[k].get(GRB_DoubleAttr_X);
1918  m_lp_coloring.push_back(value);
1919  }
1920  return;
1921 }
1922 
1923 //for debug use
1924 template<typename GraphType>
1925 pair<uint32_t, uint32_t> LPColoring<GraphType>::nonIntegerNum(const vector<GRBVar>& coloringBits, const vector<GRBVar>& vEdgeBit) const
1926 {
1927  uint32_t non_integer_num = 0;
1928  uint32_t half_num = 0;
1929  uint32_t vec_size = coloringBits.size();
1930  for(uint32_t k = 0; k < vec_size; k = k + COLORING_BITS)
1931  {
1932  for(uint32_t m = 0; m < COLORING_BITS; m++)
1933  {
1934  double value = coloringBits[k + m].get(GRB_DoubleAttr_X);
1935  if(value != 0.0 && value != 1.0)
1936  {
1937  non_integer_num++;
1938  //break;
1939  }
1940  if(value == 0.5) half_num++;
1941  cout << k+m << "-" << value << " ";
1942  }
1943  }//end for k
1944  cout << endl;
1945  cout << "NonInteger count: " << non_integer_num << ", half integer count: " << half_num << ", out of " << vec_size << " vertex variable" << endl;
1946 
1947  uint32_t non_integer_num_edge = 0;
1948  uint32_t half_num_edge = 0;
1949  uint32_t conflict_num = 0;
1950  vec_size = vEdgeBit.size();
1951  for(uint32_t k = 0; k < vec_size; ++k)
1952  {
1953  double value = vEdgeBit[k].get(GRB_DoubleAttr_X);
1954  if(value != 0.0 && value != 1.0)
1955  {
1956  non_integer_num_edge++;
1957  //break;
1958  }
1959  if(value == 0.5) half_num_edge++;
1960  if(value == 1) conflict_num++;
1961  cout << k << "-" << value << " ";
1962  }//end for k
1963  cout << endl;
1964  cout << "NonInteger count: " << non_integer_num_edge << ", half integer count: " << half_num_edge << ", out of " << vec_size << " edge variable" << endl;
1965  //cout << "conflict variable number: " << conflict_num << endl;
1966  return pair<uint32_t, uint32_t>(non_integer_num+non_integer_num_edge, half_num+half_num_edge);
1967 }
1968 
1969 template<typename GraphType>
1970 void LPColoring<GraphType>::printBoolVec(vector<bool>& vec) const
1971 {
1972  uint32_t vec_size = vec.size();
1973  for(uint32_t k = 0; k < vec_size; k++)
1974  {
1975  //cout << k << "-" << vec[k] << "; ";
1976  cout << vec[k];
1977  }//end for
1978  cout << endl;
1979 }
1980 
1981 template<typename GraphType>
1982 void LPColoring<GraphType>::printBoolArray(bool* vec, uint32_t vec_size) const
1983 {
1984  for(uint32_t k = 0; k < vec_size; k++)
1985  {
1986  //cout << k << "-" << vec[k] << "; ";
1987  cout << vec[k];
1988  }//end for
1989  cout << endl;
1990 }
1991 
1992 //print graphviz
1993 template<typename GraphType>
1994 void LPColoring<GraphType>::write_graph_dot(graph_vertex_type& v) const
1995 {
1996 #ifdef DEBUG_LPCOLORING
1997  // assert(m_graph_ptr);
1998 #endif
1999  //if(m_vertex_map.empty() || m_inverse_vertex_map.empty()) this->setMap();
2000 
2001  ofstream dot_file("graph.dot");
2002  dot_file << "graph D { \n"
2003  << " randir = LR\n"
2004  << " size=\"4, 3\"\n"
2005  << " ratio=\"fill\"\n"
2006  << " edge[style=\"bold\",fontsize=120]\n"
2007  << " node[shape=\"circle\",fontsize=120]\n";
2008 
2009  //output nodes
2010  uint32_t vertex_num = num_vertices(m_graph);
2011  //check cycles
2012  bool inCycle[vertex_num];
2013  for(uint32_t k = 0; k < vertex_num; k++) inCycle[k] = false;
2014  uint32_t cycle_cnt = m_odd_cycles.size();
2015  for(uint32_t k = 0; k < cycle_cnt; k++)
2016  {
2017  uint32_t cycle_len = m_odd_cycles[k].size();
2018  for(uint32_t m = 0; m < cycle_len; m++)
2019  {
2020  uint32_t index = m_vertex_map.at(m_odd_cycles[k][m]);
2021  inCycle[index] = true;
2022  }
2023  }//end for k
2024 
2025 #ifdef DEBUG_LPCOLORING
2026  assert(m_vertex_map.find(v) != m_vertex_map.end());
2027 #endif
2028  uint32_t start_index = m_vertex_map.at(v);
2029  for(uint32_t k = 0; k < vertex_num; k++)
2030  {
2031  if(k == start_index) dot_file << " " << k << "[shape=\"square\"";
2032  else dot_file << " " << k << "[shape=\"circle\"";
2033  //output coloring label
2034  dot_file << ",label=\"" << k << ":(" << m_lp_coloring[2*k] << "," << m_lp_coloring[2*k+1] << ")\"";
2035  if(inCycle[k]) dot_file << ",color=\"red\"]\n";
2036  else dot_file << "]\n";
2037  }//end for
2038 
2039  //output edges
2040  pair<typename graph_type::edge_iterator, typename graph_type::edge_iterator> edge_range = edges(m_graph);
2041  for(BOOST_AUTO(itr, edge_range.first); itr != edge_range.second; ++itr)
2042  {
2043  BOOST_AUTO(from, source(*itr, m_graph));
2044  uint32_t f_ind = m_vertex_map.at(from);
2045  BOOST_AUTO(to, target(*itr, m_graph));
2046  uint32_t t_ind = m_vertex_map.at(to);
2047  bool conflict_flag = false;
2048  if(m_lp_coloring[2*f_ind] == m_lp_coloring[2*t_ind] && m_lp_coloring[2*f_ind+1] == m_lp_coloring[2*t_ind+1]) {
2049  conflict_flag = true;
2050  }//end if
2051  if (m_edge_weight_map[*itr] >= 0) // conflict edge
2052  {
2053  if(conflict_flag)
2054  dot_file << " " << f_ind << "--" << t_ind << "[color=\"blue\",style=\"solid\"]\n";
2055  else
2056  dot_file << " " << f_ind << "--" << t_ind << "[color=\"black\",style=\"solid\"]\n";
2057  }
2058  else // stitch edge
2059  dot_file << " " << f_ind << "--" << t_ind << "[color=\"black\",style=\"dashed\"]\n";
2060  }
2061  dot_file << "}";
2062  dot_file.close();
2063  system("dot -Tpdf graph.dot -o graph.pdf");
2064 }
2065 
2066 //print graphviz
2067 template<typename GraphType>
2068 void LPColoring<GraphType>::write_graph_dot() const
2069 {
2070 #ifdef DEBUG_LPCOLORING
2071  // assert(m_graph_ptr);
2072 #endif
2073  //if(m_vertex_map.empty() || m_inverse_vertex_map.empty()) this->setMap();
2074 
2075  ofstream dot_file("graph.dot");
2076  dot_file << "graph D { \n"
2077  << " randir = LR\n"
2078  << " size=\"4, 3\"\n"
2079  << " ratio=\"fill\"\n"
2080  << " edge[style=\"bold\",fontsize=200]\n"
2081  << " node[shape=\"circle\",fontsize=200]\n";
2082 
2083  //output nodes
2084  uint32_t vertex_num = num_vertices(m_graph);
2085  for(uint32_t k = 0; k < vertex_num; k++)
2086  {
2087  dot_file << " " << k << "[shape=\"circle\"";
2088  //output coloring label
2089  //dot_file << ",label=\"" << k << ":(" << m_lp_coloring[2*k] << "," << m_lp_coloring[2*k+1] << ")\"";
2090  dot_file << ",label=\"(" << m_lp_coloring[2*k] << "," << m_lp_coloring[2*k+1] << ")\"";
2091  dot_file << "]\n";
2092  }//end for
2093 
2094  //output edges
2095  pair<typename graph_type::edge_iterator, typename graph_type::edge_iterator> edge_range = edges(m_graph);
2096  for(BOOST_AUTO(itr, edge_range.first); itr != edge_range.second; ++itr)
2097  {
2098  BOOST_AUTO(from, source(*itr, m_graph));
2099  uint32_t f_ind = m_vertex_map.at(from);
2100  BOOST_AUTO(to, target(*itr, m_graph));
2101  uint32_t t_ind = m_vertex_map.at(to);
2102  bool conflict_flag = false;
2103  if(m_lp_coloring[2*f_ind] == m_lp_coloring[2*t_ind] && m_lp_coloring[2*f_ind+1] == m_lp_coloring[2*t_ind+1])
2104  {
2105  conflict_flag = true;
2106  }//end if
2107  if (m_edge_weight_map[*itr] >= 0) // conflict edge
2108  {
2109  if(conflict_flag)
2110  {
2111  dot_file << " " << f_ind << "--" << t_ind << "[color=\"blue\",style=\"solid\",penwidth=3]\n";
2112  cout << "conflict (" << f_ind << ", " << t_ind << ")" << endl;
2113  }
2114  else
2115  dot_file << " " << f_ind << "--" << t_ind << "[color=\"black\",style=\"solid\",penwidth=3]\n";
2116  }
2117  else // stitch edge
2118  dot_file << " " << f_ind << "--" << t_ind << "[color=\"black\",style=\"dashed\",penwidth=3]\n";
2119  }
2120  dot_file << "}";
2121  dot_file.close();
2122  system("dot -Tpdf graph.dot -o graph.pdf");
2123 }
2124 
2125 
2126 
2127 template<typename GraphType>
2128 void LPColoring<GraphType>::write_graph_color() const
2129 {
2130 #ifdef DEBUG_LPCOLORING
2131  // assert(m_graph_ptr);
2132 #endif
2133  //if(m_vertex_map.empty() || m_inverse_vertex_map.empty()) this->setMap();
2134  ofstream dot_file("color_graph.dot");
2135  dot_file << "graph D { \n"
2136  << " randir = LR\n"
2137  << " size=\"4, 3\"\n"
2138  << " ratio=\"fill\"\n"
2139  << " edge[style=\"bold\",fontsize=120]\n"
2140  << " node[shape=\"circle\",fontsize=120]\n";
2141 
2142  //output nodes
2143  uint32_t vertex_num = num_vertices(m_graph);
2144  //check cycles
2145  bool inCycle[vertex_num];
2146  for(uint32_t k = 0; k < vertex_num; k++) inCycle[k] = false;
2147  uint32_t cycle_cnt = m_odd_cycles.size();
2148  for(uint32_t k = 0; k < cycle_cnt; k++)
2149  {
2150  uint32_t cycle_len = m_odd_cycles[k].size();
2151  for(uint32_t m = 0; m < cycle_len; m++)
2152  {
2153  uint32_t index = m_vertex_map.at(m_odd_cycles[k][m]);
2154  inCycle[index] = true;
2155  }
2156  }//end for k
2157 
2158  for(uint32_t k = 0; k < vertex_num; k++)
2159  {
2160  dot_file << " " << k << "[shape=\"circle\"";
2161  //output coloring label
2162  dot_file << ",label=\"" << m_coloring.at(m_inverse_vertex_map.at(k)) << "\"";
2163  if(inCycle[k]) dot_file << ",color=\"red\"]\n";
2164  else dot_file << "]\n";
2165  }//end for
2166 
2167  //output edges
2168  pair<typename graph_type::edge_iterator, typename graph_type::edge_iterator> edge_range = edges(m_graph);
2169  for(BOOST_AUTO(itr, edge_range.first); itr != edge_range.second; ++itr)
2170  {
2171  BOOST_AUTO(from, source(*itr, m_graph));
2172  uint32_t f_ind = m_vertex_map.at(from);
2173  BOOST_AUTO(to, target(*itr, m_graph));
2174  uint32_t t_ind = m_vertex_map.at(to);
2175 
2176  if (m_edge_weight_map[*itr] >= 0) // conflict edge
2177  {
2178  if(m_coloring.at(from) != m_coloring.at(to))
2179  dot_file << " " << f_ind << "--" << t_ind << "[color=\"black\",style=\"solid\",penwidth=3]\n";
2180  else
2181  dot_file << " " << f_ind << "--" << t_ind << "[color=\"blue\",style=\"solid\",penwidth=3]\n";
2182  }
2183  else // stitch edge
2184  {
2185  if(m_coloring.at(from) == m_coloring.at(to))
2186  dot_file << " " << f_ind << "--" << t_ind << "[color=\"black\",style=\"dashed\",penwidth=3]\n";
2187  else
2188  dot_file << " " << f_ind << "--" << t_ind << "[color=\"blue\",style=\"dashed\",penwidth=3]\n";
2189  }
2190  }
2191  dot_file << "}";
2192  dot_file.close();
2193  system("dot -Tpdf color_graph.dot -o color_graph.pdf");
2194 }
2195 
2196 } // namespace coloring
2197 } // namespace algorithms
2198 } // namespace limbo
2199 
2200 #endif
2201 
limbo::algorithms::coloring::LPColoring::m_constrs_num
uint32_t m_constrs_num
record number of constraints
Definition: LPColoring.h:223
limbo
namespace for Limbo
Definition: GraphUtility.h:22
limbo::algorithms::coloring::Coloring::m_graph
graph_type const & m_graph
initial graph
Definition: Coloring.h:229
limbo::algorithms::coloring::LPColoring::LPColoring
LPColoring(graph_type const &g)
constructor
Definition: LPColoring.h:229
limbo::algorithms::coloring::Coloring::EdgeHashType::operator()
std::size_t operator()(graph_edge_type const &e) const
Definition: Coloring.h:142
limbo::algorithms::coloring::Coloring::m_color_num
ColorNumType m_color_num
number of colors
Definition: Coloring.h:232
limbo::algorithms::coloring::Coloring::m_stitch_weight
double m_stitch_weight
stitch weight
Definition: Coloring.h:233
Generated on Wed Jul 3 2019 09:58:19 for Limbo by   doxygen 1.8.8