00001
00002
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
00034 negedges[neg++] = bp_nodes[i].edge_to_parent;
00035 } else {
00036
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)
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
00102
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
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
00123
00124 edges_to_create[i] |= 1;
00125 ++bp_nodes[t].degree;
00126 ++bp_nodes[h].degree;
00127 }
00128 if (xi < minviol1) {
00129
00130
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) {
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) {
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
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
00223
00224
00225
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
00264 double limit = minviol1;
00265
00266
00267
00268
00269
00270
00271
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
00289 if (! generate_all_cuts) {
00290
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
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 }