/home/coin/SVN-release/OS-2.4.2/Bcp/examples/MaxCut/CG/MC_shortest_path_cycle.cpp

Go to the documentation of this file.
00001 // Copyright (C) 2000, International Business Machines
00002 // Corporation and others.  All Rights Reserved.
00003 
00004 #include <queue>
00005 #include <vector>
00006 
00007 #include "CoinSort.hpp"
00008 
00009 #include "MC.hpp"
00010 #include "MC_cut.hpp"
00011 
00012 #include "BCP_matrix.hpp"
00013 
00014 //#############################################################################
00015 
00016 void
00017 MC_create_shortest_path_cut(const int n, const int sp_tree_root,
00018                             const MC_graph_edge* edges,
00019                             const MC_path_node* bp_nodes,
00020                             BCP_vec<BCP_cut*>& new_cuts,
00021                             BCP_vec<BCP_row*>& new_rows,
00022                             int * itmp, double * dtmp)
00023 {
00024   int pos = 0, neg = 0;
00025   int* posedges = itmp;
00026   int* negedges = itmp + n;
00027 
00028   int i = sp_tree_root + n;
00029   int parent;
00030   do {
00031     parent = bp_nodes[i].parent;
00032     if ((parent < n && i < n) || (parent >= n && i >= n)) {
00033       // side edge thus not in F
00034       negedges[neg++] = bp_nodes[i].edge_to_parent;
00035     } else {
00036       // crossing edge thus in F
00037       posedges[pos++] = bp_nodes[i].edge_to_parent;
00038     }
00039     i = parent;
00040   } while (parent != sp_tree_root);
00041 
00042   if ((pos % 2) != 1)
00043     abort();
00044   const int len = pos + neg;
00045   CoinDisjointCopyN(negedges, neg, posedges + pos);
00046   new_cuts.push_back(new MC_cycle_cut(posedges, pos, len));
00047   CoinFillN(dtmp, pos, 1.0);
00048   CoinFillN(dtmp + pos, neg, -1.0);
00049   CoinSort_2(posedges, posedges + len, dtmp);
00050   new_rows.push_back(new BCP_row(posedges, posedges + len,
00051                                  dtmp, dtmp + len, -BCP_DBL_MAX, pos - 1.0));
00052 }
00053 
00054 //#############################################################################
00055 
00056 int
00057 MC_find_components(const int nodenum, const MC_path_node* bp_nodes,
00058                    int * component)
00059 {
00060   int* node_stack = new int[nodenum];
00061   int stack_size = 0;
00062   int comp_num = 0;
00063   int i, j;
00064   for (i = 0; i < nodenum; ++i) {
00065     if (component[i] >= 0) // it's already encountered
00066       continue;
00067     node_stack[stack_size++] = i;
00068     ++comp_num;
00069     while (stack_size > 0) {
00070       const int node = node_stack[--stack_size];
00071       component[node] = i;
00072       const MC_path_adj_entry* neighbors = bp_nodes[node].adj_list;
00073       for (j = bp_nodes[node].degree - 1; j >= 0; --j) {
00074         const int other = neighbors[j].neighbor;
00075         if (component[other] == -1) {
00076           component[other] = i;
00077           node_stack[stack_size++] = other;
00078         }
00079       }
00080     }
00081   }
00082   delete[] node_stack;
00083   return comp_num;
00084 }
00085 
00086 //#############################################################################
00087 
00088 void
00089 MC_generate_shortest_path_cycles(const MC_problem& mc, const double* x,
00090                                  const bool generate_all_cuts,
00091                                  const double minviol,
00092                                  BCP_vec<BCP_cut*>& new_cuts,
00093                                  BCP_vec<BCP_row*>& new_rows)
00094 {
00095   const double minviol1 = 1.0 - minviol;
00096 
00097   const int n = mc.num_nodes;
00098   const int m = mc.num_edges;
00099   double* dtmp = new double[n];
00100   int* itmp = new int[4 * n + m];
00101   // the first 2n entries will be used in MC_generate_shortest_path_cycles as
00102   // temporaries
00103   int* component = itmp + (2 * n);
00104   int* edges_to_create = itmp + (4 * n);
00105   const MC_graph_edge* edges = mc.edges;
00106 
00107   MC_path_node* bp_nodes = new MC_path_node[2 * n];
00108   MC_path_adj_entry* bp_all_adj_list = new MC_path_adj_entry[8 * m];
00109 
00110   int i, j;
00111 
00112   for (i = 0; i < n; ++i)
00113     bp_nodes[i].degree = 0;
00114 
00115   // Build the bipartite graph
00116   for (i = 0; i < m; ++i) {
00117     const double xi = x[i];
00118     const int t = edges[i].tail;
00119     const int h = edges[i].head;
00120     edges_to_create[i] = 0;
00121     if (xi >= minviol) {
00122       // the length of the crossing edges is < 1 - minviol, so there can be a
00123       // violated cycle through them, so they exist
00124       edges_to_create[i] |= 1;
00125       ++bp_nodes[t].degree;
00126       ++bp_nodes[h].degree;
00127     }
00128     if (xi < minviol1) {
00129       // the length of the side edges is < 1 - minviol, so there can be a
00130       // violated cycle through them, so they exist
00131       edges_to_create[i] |= 2;
00132       ++bp_nodes[t].degree;
00133       ++bp_nodes[h].degree;
00134     }
00135   }
00136 
00137   int total_degree = 0;
00138   for (i = 0; i < n; ++i) {
00139     bp_nodes[i].adj_list = bp_all_adj_list + total_degree;
00140     total_degree += bp_nodes[i].degree;
00141     bp_nodes[i + n].adj_list = bp_all_adj_list + total_degree;
00142     total_degree += bp_nodes[i].degree;
00143   }
00144 
00145   for (i = 0; i < n; ++i)
00146     bp_nodes[i].degree = 0;
00147 
00148   for (i = 0; i < m; ++i) {
00149     const double xi = CoinMax(x[i], 0.0);
00150     const double xi1 = CoinMax(1 - xi, 0.0);
00151     const int t = edges[i].tail;
00152     const int h = edges[i].head;
00153     const int nt = n + t;
00154     const int nh = n + h;
00155     if ((edges_to_create[i] & 1) != 0) { // create the crossing edges
00156       const int degt = bp_nodes[t].degree;
00157       const int degh = bp_nodes[h].degree;
00158       MC_path_adj_entry& t_entry = bp_nodes[t].adj_list[degt];
00159       MC_path_adj_entry& h_entry = bp_nodes[h].adj_list[degh];
00160       MC_path_adj_entry& nt_entry = bp_nodes[nt].adj_list[degt];
00161       MC_path_adj_entry& nh_entry = bp_nodes[nh].adj_list[degh];
00162       
00163       t_entry.neighbor = nh;
00164       h_entry.neighbor = nt;
00165       nt_entry.neighbor = h;
00166       nh_entry.neighbor = t;
00167 
00168       t_entry.orig_edge = i;
00169       h_entry.orig_edge = i;
00170       nt_entry.orig_edge = i;
00171       nh_entry.orig_edge = i;
00172 
00173       t_entry.cost = xi1;
00174       h_entry.cost = xi1;
00175       nt_entry.cost = xi1;
00176       nh_entry.cost = xi1;
00177 
00178       ++bp_nodes[t].degree;
00179       ++bp_nodes[h].degree;
00180     }
00181     if ((edges_to_create[i] & 2) != 0) { // create the side edges
00182       const int degt = bp_nodes[t].degree;
00183       const int degh = bp_nodes[h].degree;
00184       MC_path_adj_entry& t_entry = bp_nodes[t].adj_list[degt];
00185       MC_path_adj_entry& h_entry = bp_nodes[h].adj_list[degh];
00186       MC_path_adj_entry& nt_entry = bp_nodes[nt].adj_list[degt];
00187       MC_path_adj_entry& nh_entry = bp_nodes[nh].adj_list[degh];
00188 
00189       t_entry.neighbor = h;
00190       h_entry.neighbor = t;
00191       nt_entry.neighbor = nh;
00192       nh_entry.neighbor = nt;
00193 
00194       t_entry.orig_edge = i;
00195       h_entry.orig_edge = i;
00196       nt_entry.orig_edge = i;
00197       nh_entry.orig_edge = i;
00198 
00199       t_entry.cost = xi;
00200       h_entry.cost = xi;
00201       nt_entry.cost = xi;
00202       nh_entry.cost = xi;
00203 
00204       ++bp_nodes[t].degree;
00205       ++bp_nodes[h].degree;
00206     }
00207   }
00208   for (i = 0; i < n; ++i)
00209     bp_nodes[n+i].degree = bp_nodes[i].degree;
00210 
00211   // Create a priority queue that will hold the shortest path tree
00212   std::priority_queue<MC_path_node*,
00213                       std::vector<MC_path_node*>,
00214                       MC_path_node_ptr_compare> pqueue;
00215 
00216   CoinFillN(component, 2 * n, -1);
00217 
00218   const int comp_num = MC_find_components(2*n, bp_nodes, component);
00219   printf("MC: sp based exact: there are %i components in the graph\n",
00220          comp_num);
00221 
00222   /* Check the nodes where we already have a sufficiently violated cycle
00223      through the node. If such a node and it's pair is in the same component
00224      then set the component id of one to be -1 so we skip the shortest path
00225      computation for them. */
00226   double viol;
00227   int skip = 0;
00228   for (i = new_cuts.size() - 1; i >= 0; --i) {
00229     const MC_cycle_cut* cut = dynamic_cast<const MC_cycle_cut*>(new_cuts[i]);
00230     if (! cut)
00231       continue;
00232     viol = - (cut->pos_edges - 1);
00233     for (j = 0; j < cut->pos_edges; ++j) {
00234       viol += x[cut->edges[j]];
00235     }
00236     for ( ; j < cut->cycle_len; ++j) {
00237       viol -= x[cut->edges[j]];
00238     }
00239     if (viol <= minviol)
00240       continue;
00241     for (j = 0; j < cut->cycle_len; ++j) {
00242       const int t = edges[cut->edges[j]].tail;
00243       const int h = edges[cut->edges[j]].head;
00244       if (component[t] == component[n + t]) {
00245         ++skip;
00246         component[t] = -2;
00247       }
00248       if (component[h] == component[n + h]) {
00249         ++skip;
00250         component[h] = -2;
00251       }
00252     }
00253   }
00254   printf("MC: sp based exact: skipping %i sp computations\n", skip);
00255   printf("                    because of already found cycles\n");
00256 
00257   int split_num = 0;
00258   for (i = 0; i < n; ++i) {
00259     if (component[i] != component[n + i]) {
00260       ++split_num;
00261       continue;
00262     }
00263     // do the shortest path
00264     double limit = minviol1;
00265 
00266     // set the original distance of everything to 1.0 (that's larger than the
00267     // longest edge, and we will anyway terminate every path that's longer
00268     // than limit). Note that a node can be in the priority queue more than
00269     // once, but it's going to be processed only once. All this is necessary
00270     // because the std priority queue installation cannot "relocate" a node
00271     // when it's value changes.
00272     for (j = 2*n - 1; j >= 0; --j) {
00273       bp_nodes[j].distance = 1.0;
00274       bp_nodes[j].processed = false;
00275     }
00276     bp_nodes[i].distance = 0.0;
00277 
00278     pqueue.push(bp_nodes + i);
00279 
00280     while (! pqueue.empty()) {
00281       MC_path_node* this_node_nonconst = pqueue.top();
00282       const MC_path_node* this_node = this_node_nonconst;
00283       pqueue.pop();
00284       if (this_node->processed)
00285         continue;
00286       const int this_node_index = this_node - bp_nodes;
00287       if (this_node_index == n+i) {
00288         // got the shortest path to the pair, need not go further
00289         if (! generate_all_cuts) {
00290           // wee needed only the best cut for each node, here it is
00291           MC_create_shortest_path_cut(n, this_node_index, edges, bp_nodes,
00292                                       new_cuts, new_rows, itmp, dtmp);
00293         }
00294         break;
00295       }
00296       this_node_nonconst->processed = true;
00297       const double this_dist = this_node->distance;
00298       const MC_path_adj_entry* this_alist = this_node->adj_list;
00299       for (j = this_node->degree - 1; j >= 0; --j) {
00300         const double neighborcost = this_dist + this_alist[j].cost;
00301         if (neighborcost >= limit)
00302           continue;
00303         const int neighbor = this_alist[j].neighbor;
00304         if (bp_nodes[neighbor].distance > neighborcost) {
00305           MC_path_node* neighbor_node = bp_nodes + neighbor;
00306           pqueue.push(neighbor_node);
00307           neighbor_node->distance = neighborcost;
00308           neighbor_node->parent = this_node_index;
00309           neighbor_node->edge_to_parent = this_alist[j].orig_edge;
00310           if (neighbor == n+i) {
00311             // found a new violated cycle...
00312             limit = neighborcost;
00313             if (generate_all_cuts) {
00314               MC_create_shortest_path_cut(n, i, edges, bp_nodes,
00315                                           new_cuts, new_rows, itmp, dtmp);
00316             }
00317           }
00318         }
00319       }
00320     }
00321   }
00322 
00323   printf("MC: sp based exact: computed %i sp out of %i\n", n - split_num, n);
00324 
00325   delete[] bp_all_adj_list;
00326   delete[] bp_nodes;
00327   delete[] itmp;
00328   delete[] dtmp;
00329 }

Generated on Wed Nov 30 03:03:46 2011 by  doxygen 1.4.7