00001
00002
00003
00004
00005 #include <cstdio>
00006 #include <fstream>
00007
00008 #include "CoinHelperFunctions.hpp"
00009
00010 #include "BCP_math.hpp"
00011 #include "BCP_string.hpp"
00012 #include "BCP_buffer.hpp"
00013 #include "BCP_vector.hpp"
00014
00015 #include "MC.hpp"
00016 #include "MC_solution.hpp"
00017
00018 #define MC_PRINT_ITER_COUNT 0
00019
00020 #define MC_TEST_COST_AFTER_EVERY_SWITCH 0
00021
00022
00023
00024 BCP_buffer&
00025 MC_solution::pack(BCP_buffer& buf) const {
00026 buf.pack(cost).pack(sig);
00027 return buf;
00028 }
00029
00030
00031
00032 BCP_buffer&
00033 MC_solution::unpack(BCP_buffer& buf) {
00034 buf.unpack(cost).unpack(sig);
00035 return buf;
00036 }
00037
00038
00039
00040 MC_solution&
00041 MC_solution::operator=(const MC_solution& sol) {
00042 cost = sol.cost;
00043 sig = sol.sig;
00044 return *this;
00045 }
00046
00047
00048
00049 double
00050 MC_solution::compute_cost(const int m, const MC_graph_edge* edges)
00051 {
00052 cost = 0.0;
00053 for (int i = 0; i < m; ++i) {
00054 const MC_graph_edge& edge = edges[i];
00055 if (sig[edge.tail] != sig[edge.head]) {
00056 cost += edge.cost;
00057 }
00058 }
00059 return cost;
00060 }
00061
00062
00063
00064 static inline double
00065 MC_switch_cost(const MC_adjacency_entry* adj_list, const int degree,
00066 const int* sig, const int this_sig)
00067 {
00068 double sum = 0;
00069 for (int j = degree - 1; j >= 0; --j) {
00070 #if 1
00071
00072 sum += sig[adj_list[j].neighbor] * this_sig * adj_list[j].cost;
00073 #else
00074 if (sig[adj_list[j].neighbor] == this_sig) {
00075 sum += adj_list[j].cost;
00076 } else {
00077 sum -= adj_list[j].cost;
00078 }
00079 #endif
00080 }
00081 return sum;
00082 }
00083
00084
00085
00086 double
00087 MC_solution::switch_improve(const MC_problem& mc, const int maxiter)
00088 {
00089 const int n = mc.num_nodes;
00090 bool flag = true;
00091 int i, iter;
00092 double newcost = cost;
00093 for (iter = 0; iter < maxiter && flag; ++iter) {
00094 flag = false;
00095 for (i = 0; i < n; ++i) {
00096 const double sum = MC_switch_cost(mc.nodes[i].adj_list,
00097 mc.nodes[i].degree,
00098 sig.begin(), sig[i]);
00099 if (sum < 0) {
00100 sig[i] *= -1;
00101 newcost += sum;
00102 flag = true;
00103 #if MC_TEST_COST_AFTER_EVERY_SWITCH
00104 compute_cost(mc.num_edges, mc.edges);
00105 if (cost != newcost)
00106 abort();
00107 #endif
00108 }
00109 }
00110 }
00111 #if MC_PRINT_ITER_COUNT
00112 printf("MC_solution: switch_improve: %i.\n", iter);
00113 #endif
00114
00115 compute_cost(mc.num_edges, mc.edges);
00116 if (cost != newcost)
00117 abort();
00118
00119 return cost;
00120 }
00121
00122
00123
00124
00125
00126
00127
00128 double
00129 MC_solution::edge_switch_improve(const MC_problem& mc, const int maxiter)
00130 {
00131 const int m = mc.num_edges;
00132 bool flag = true;
00133 int k, iter;
00134 double newcost = cost;
00135 const MC_graph_edge* edges = mc.edges;
00136
00137 for (iter = 0; iter < maxiter && flag; ++iter) {
00138 flag = false;
00139 for (k = 0; k < m; ++k) {
00140 const int i = edges[k].tail;
00141 const int j = edges[k].head;
00142 const double costk = 2 * sig[i] * sig[j] * edges[k].cost;
00143
00144
00145
00146 if (costk <= 0.0)
00147 continue;
00148 const double sum =
00149 MC_switch_cost(mc.nodes[i].adj_list, mc.nodes[i].degree,
00150 sig.begin(), sig[i]) +
00151 MC_switch_cost(mc.nodes[j].adj_list, mc.nodes[j].degree,
00152 sig.begin(), sig[j]) -
00153 costk;
00154 if (sum < 0) {
00155 sig[i] *= -1;
00156 sig[j] *= -1;
00157 newcost += sum;
00158 flag = true;
00159 #if MC_TEST_COST_AFTER_EVERY_SWITCH
00160 compute_cost(mc.num_edges, mc.edges);
00161 if (cost != newcost)
00162 abort();
00163 #endif
00164 }
00165 }
00166 }
00167 #if MC_PRINT_ITER_COUNT
00168 printf("MC_solution: edge_switch_improve: %i.\n", iter);
00169 #endif
00170
00171 compute_cost(mc.num_edges, mc.edges);
00172 if (cost != newcost)
00173 abort();
00174
00175 return cost;
00176 }
00177
00178
00179
00180
00181
00182
00183 double
00184 MC_solution::ising_with_external_edge_switch_improve(const MC_problem& mc,
00185 const int maxiter)
00186 {
00187 const int m = mc.num_edges;
00188 const int grid_edge_num = m - mc.num_nodes + 1;
00189 bool flag = true;
00190 int k, iter;
00191 double newcost = cost;
00192 const MC_graph_edge* edges = mc.edges;
00193
00194 const int plus_node = mc.num_nodes - 1;
00195 const MC_adjacency_entry* plus_node_adj_list = mc.nodes[plus_node].adj_list;
00196
00197 for (iter = 0; iter < maxiter && flag; ++iter) {
00198 flag = false;
00199
00200 for (k = 0; k < grid_edge_num; ++k) {
00201 const int i = edges[k].tail;
00202 const int j = edges[k].head;
00203 const double costk = 2 * sig[i] * sig[j] * edges[k].cost;
00204
00205
00206
00207 if (costk <= 0.0)
00208 continue;
00209 const double sum =
00210 MC_switch_cost(mc.nodes[i].adj_list, mc.nodes[i].degree,
00211 sig.begin(), sig[i]) +
00212 MC_switch_cost(mc.nodes[j].adj_list, mc.nodes[j].degree,
00213 sig.begin(), sig[j]) -
00214 costk;
00215 if (sum < 0) {
00216 sig[i] *= -1;
00217 sig[j] *= -1;
00218 newcost += sum;
00219 flag = true;
00220 #if MC_TEST_COST_AFTER_EVERY_SWITCH
00221 compute_cost(mc.num_edges, mc.edges);
00222 if (cost != newcost)
00223 abort();
00224 #endif
00225 }
00226 }
00227
00228
00229 double switch_cost_of_plus_node =
00230 MC_switch_cost(plus_node_adj_list, mc.nodes[plus_node].degree,
00231 sig.begin(), sig[plus_node]);
00232 for ( ; k < m; ++k) {
00233 const int i = edges[k].tail;
00234 const double costk = 2 * sig[i] * sig[plus_node] * edges[k].cost;
00235
00236
00237
00238 if (costk <= 0.0)
00239 continue;
00240 const double sum =
00241 MC_switch_cost(mc.nodes[i].adj_list, mc.nodes[i].degree,
00242 sig.begin(), sig[i]) +
00243 switch_cost_of_plus_node -
00244 costk;
00245 if (sum < 0) {
00246 sig[i] *= -1;
00247 sig[plus_node] *= -1;
00248 switch_cost_of_plus_node = - switch_cost_of_plus_node + costk;
00249 newcost += sum;
00250 flag = true;
00251 #if MC_TEST_COST_AFTER_EVERY_SWITCH
00252 compute_cost(mc.num_edges, mc.edges);
00253 if (cost != newcost)
00254 abort();
00255 #endif
00256 }
00257 }
00258 }
00259 #if MC_PRINT_ITER_COUNT
00260 printf("MC_solution: edge_switch_improve: %i.\n", iter);
00261 #endif
00262
00263 compute_cost(mc.num_edges, mc.edges);
00264 if (cost != newcost)
00265 abort();
00266
00267 return cost;
00268 }
00269
00270
00271
00272 double
00273 MC_solution::structure_switch_improve(const MC_problem& mc,
00274 const int struct_ind, const int maxiter)
00275 {
00276 if (mc.num_switch_structures[struct_ind] == 0)
00277 return cost;
00278
00279 const int s = mc.num_switch_structures[struct_ind];
00280 const MC_switch_structure* sw_structs = mc.switch_structures[struct_ind];
00281 bool flag = true;
00282 int i, k, iter;
00283 double sum;
00284 double newcost = cost;
00285 const MC_graph_edge* edges = mc.edges;
00286
00287 for (iter = 0; iter < maxiter && flag; ++iter) {
00288 flag = false;
00289 for (k = 0; k < s; ++k) {
00290 const MC_switch_structure& sw_struct = sw_structs[k];
00291 const int * nb = sw_struct.neighbors;
00292 sum = 0.0;
00293 for (i = sw_struct.num_neighbors - 1; i >= 0; --i) {
00294 const MC_graph_edge& e = edges[nb[i]];
00295 sum += sig[e.tail] * sig[e.head] * e.cost;
00296 }
00297 if (sum < 0) {
00298 const int * nodes = sw_struct.nodes;
00299 for (i = sw_struct.num_nodes - 1; i >= 0; --i) {
00300 sig[nodes[i]] *= -1;
00301 }
00302 newcost += sum;
00303 flag = true;
00304 #if MC_TEST_COST_AFTER_EVERY_SWITCH
00305 compute_cost(mc.num_edges, mc.edges);
00306 if (cost != newcost)
00307 abort();
00308 #endif
00309 }
00310 }
00311 }
00312 #if MC_PRINT_ITER_COUNT
00313 printf("MC_solution: sw_struct_improve: %i.\n", iter);
00314 #endif
00315
00316 compute_cost(mc.num_edges, mc.edges);
00317 if (cost != newcost)
00318 abort();
00319
00320 return cost;
00321 }
00322
00323
00324
00325 double
00326 MC_solution::lk_switch_improve(const MC_problem& mc, const int maxiter)
00327 {
00328 const int n = mc.num_nodes;
00329 int i, start, end, iter = 0;
00330 double sum;
00331 double best_sum;
00332
00333 start = 0;
00334 for (iter = 0; iter < 1000 * maxiter; ++iter) {
00335 best_sum = 0.0;
00336 sum = 0.0;
00337
00338
00339 for (i = start; i < n; ++i) {
00340 sum = MC_switch_cost(mc.nodes[i].adj_list, mc.nodes[i].degree,
00341 sig.begin(), sig[i]);
00342 if (sum < 0.0)
00343 break;
00344 }
00345 if (sum == 0.0) {
00346 for (i = 0; i < start; ++i) {
00347 sum = MC_switch_cost(mc.nodes[i].adj_list, mc.nodes[i].degree,
00348 sig.begin(), sig[i]);
00349 if (sum < 0.0)
00350 break;
00351 }
00352 }
00353
00354 if (sum >= 0.0)
00355 break;
00356
00357 sig[i] *= -1;
00358 start = i;
00359 best_sum = sum;
00360 end = i;
00361
00362
00363 for (i = start+1; i < n; ++i) {
00364 sum += MC_switch_cost(mc.nodes[i].adj_list, mc.nodes[i].degree,
00365 sig.begin(), sig[i]);
00366 if (sum >= 0)
00367 break;
00368 sig[i] *= -1;
00369 if (sum < best_sum) {
00370 best_sum = sum;
00371 end = i;
00372 }
00373 }
00374 if (sum < 0) {
00375 for (i = 0; i < start; ++i) {
00376 sum += MC_switch_cost(mc.nodes[i].adj_list, mc.nodes[i].degree,
00377 sig.begin(), sig[i]);
00378 if (sum >= 0)
00379 break;
00380 sig[i] *= -1;
00381 if (sum < best_sum) {
00382 best_sum = sum;
00383 end = i;
00384 }
00385 }
00386 }
00387
00388
00389
00390
00391
00392 if (end > i) {
00393 for (--i; i >= 0; --i)
00394 sig[i] *= -1;
00395 i = n;
00396 }
00397 for (--i; i > end; --i)
00398 sig[i] *= -1;
00399
00400 const double newcost = cost + best_sum;
00401 compute_cost(mc.num_edges, mc.edges);
00402 if (cost != newcost)
00403 abort();
00404
00405 start = (end + 1) %n;
00406 }
00407 #if MC_PRINT_ITER_COUNT
00408 printf("MC_solution: lk_switch_improve: %i.\n", iter);
00409 #endif
00410
00411 return cost;
00412 }
00413
00414
00415
00416 void
00417 MC_solution::display(const BCP_string& fname) const
00418 {
00419 const int num_nodes = sig.size();
00420 int objval = static_cast<int>(objective_value());
00421
00422 printf("MC_solution:\n");
00423 printf("%i\n", objval);
00424 for (int i = 0; i < num_nodes; ++i) {
00425 printf("%2i\n", sig[i]);
00426 }
00427 printf("\n");
00428
00429 if (fname.length() != 0) {
00430 std::ofstream ofile(fname.c_str());
00431 ofile << "MC_solution:\n";
00432 ofile << objval << "\n";
00433 for (int i = 0; i < num_nodes; ++i) {
00434 ofile << sig[i] << "\n";
00435 }
00436 }
00437 }
00438
00439
00440
00441 MC_solution::MC_solution(const BCP_vec<int>& sign,
00442 const MC_problem& mc,
00443 const int heurswitchround,
00444 const bool do_edge_switch_heur,
00445 const int struct_switch_heur) :
00446 sig(sign)
00447 {
00448 const int m = mc.num_edges;
00449 const MC_graph_edge* edges = mc.edges;
00450 const int num_struct_type = mc.num_structure_type;
00451
00452 const double percentage[3] = {0.05, 0.01, 0.002};
00453
00454 BCP_vec<int> bestsig;
00455 double bestcost = BCP_DBL_MAX;
00456
00457 int i;
00458 bool do_plain_switch = true;
00459 bool do_edge_switch = true;
00460 bool* do_structure_switch = new bool[num_struct_type];
00461 bool do_switch = true;
00462
00463 printf("MC: Heur Sol Improvement: ");
00464 for (int k = 0; k < 1; ++k) {
00465 compute_cost(m, edges);
00466
00467 BCP_vec<double> obj;
00468 BCP_vec<char> type;
00469
00470 obj.push_back(objective_value());
00471 type.push_back('o');
00472
00473 do_plain_switch = true;
00474 do_edge_switch = true;
00475 CoinFillN(do_structure_switch, num_struct_type, true);
00476 do_switch = true;
00477
00478 while (do_switch) {
00479 if (do_plain_switch) {
00480 switch_improve(mc, heurswitchround);
00481 if (obj.back() > objective_value()) {
00482 obj.push_back(objective_value());
00483 type.push_back('P');
00484
00485 do_edge_switch = true;
00486 CoinFillN(do_structure_switch, num_struct_type, true);
00487 }
00488 do_plain_switch = false;
00489 }
00490 if (do_edge_switch) {
00491 if (do_edge_switch_heur) {
00492 if (!mc.ising_triangles)
00493 edge_switch_improve(mc, heurswitchround);
00494 else
00495 ising_with_external_edge_switch_improve(mc, heurswitchround);
00496 if (obj.back() > objective_value()) {
00497 obj.push_back(objective_value());
00498 type.push_back('E');
00499
00500 do_plain_switch = true;
00501 CoinFillN(do_structure_switch, num_struct_type, true);
00502 }
00503 }
00504 do_edge_switch = false;
00505 }
00506 for (i = 0; i < num_struct_type; ++i) {
00507 if (do_structure_switch[i]) {
00508 if ((struct_switch_heur & (1 << i)) != 0) {
00509 structure_switch_improve(mc, i, heurswitchround);
00510 if (obj.back() > objective_value()) {
00511 obj.push_back(objective_value());
00512 type.push_back('0'+i);
00513
00514 do_plain_switch = true;
00515 do_edge_switch = true;
00516 CoinFillN(do_structure_switch, num_struct_type, true);
00517 }
00518 }
00519 }
00520 do_structure_switch[i] = false;
00521 }
00522
00523 do_switch = do_plain_switch || do_edge_switch;
00524 for (i = 0; i < num_struct_type; ++i)
00525 do_switch |= do_structure_switch[i];
00526 }
00527 printf(" %.0f", obj[0]);
00528 const int rounds = obj.size();
00529 for (i = 1; i < rounds; ++i) {
00530 printf(" / %c %.0f", type[i], obj[i]);
00531 }
00532 printf("\n");
00533
00534 if (cost < bestcost) {
00535 bestsig = sig;
00536 bestcost = cost;
00537 }
00538
00539 for (i = sig.size() - 1; i >= 0; --i) {
00540 if (CoinDrand48() < percentage[k])
00541 sig[i] = -sig[i];
00542 }
00543 }
00544
00545 sig = bestsig;
00546 cost = bestcost;
00547
00548 delete[] do_structure_switch;
00549 }