/home/coin/SVN-release/OS-2.4.1/Bcp/examples/MaxCut/Member/MC_solution.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 <cfloat>
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     // Because of pipelined processors, this is probably much faster...
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 // This routine must be invoked after switch_improve(), since to speed things
00125 // up we assume that no single node switch can help, i.e., MC_switch_cost() is
00126 // non-negative for all nodes.
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       // we assume that no single node switch can help at this point ==>
00144       // The MC_switch_cost values will be non-negative ==> costk > 0 is a
00145       // necessary condition to find improvement.
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 // This routine must be invoked after switch_improve(), since to speed things
00180 // up we assume that no single node switch can help, i.e., MC_switch_cost() is
00181 // non-negative for all nodes.
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     // First try to switch edges not going to the extra node
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       // we assume that no single node switch can help at this point ==>
00205       // The MC_switch_cost values will be non-negative ==> costk > 0 is a
00206       // necessary condition to find improvement.
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     // Now try to switch edges going to the extra node
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       // we assume that no single node switch can help at this point ==>
00236       // The MC_switch_cost values will be non-negative ==> costk > 0 is a
00237       // necessary condition to find improvement.
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     // Find the first node that does improve
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) // No individual node can be switched
00355       break;
00356 
00357     sig[i] *= -1;
00358     start = i;
00359     best_sum = sum;
00360     end = i;
00361     
00362     // Do the Lin-Kerninghan thingy  
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     // Now best_sum must be < 0 and we got to switch the [start,end] part
00389     // only, i.e., we got to reverse the switching of (end,i).
00390     // Note that i cannot be start since when everything is switched then sum
00391     // is 0.0.
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           // improved, got to do the others, too
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             // improved, got to do the others, too
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               // improved, got to do the others, too
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 }

Generated on Thu Nov 10 03:05:37 2011 by  doxygen 1.4.7