angel  mercurial changeset:
angel_tools.cpp
Go to the documentation of this file.
00001 // $Id: angel_tools.cpp,v 1.4 2008/02/28 14:57:33 gottschling Exp $
00002 /*
00003 #############################################################
00004 # This file is part of angel released under the BSD license #
00005 # The full COPYRIGHT notice can be found in the top         #
00006 # level directory of the angel distribution                 #
00007 #############################################################
00008 */
00009 
00010 #include "angel/include/angel_tools.hpp"
00011 
00012 #include <queue>
00013 
00014 #include "angel/include/angel_exceptions.hpp"
00015 #include "angel/include/angel_io.hpp"
00016 #include "angel/include/eliminations.hpp"
00017 
00018 namespace angel {
00019 
00020 using namespace std;
00021 using namespace boost;
00022 
00023 bool reachable (const c_graph_t::vertex_t src,
00024                 const c_graph_t::vertex_t tgt,
00025                 c_graph_t& angelLCG) {
00026   property_map<pure_c_graph_t, VertexVisited>::type visited = get(VertexVisited(), angelLCG);
00027 
00028   c_graph_t::oei_t oei, oe_end;
00029   for (tie (oei, oe_end) = out_edges (src, angelLCG); oei != oe_end; ++oei) {
00030     if (target(*oei, angelLCG) == tgt) return true;
00031     if (!visited[target(*oei, angelLCG)]) { // call recursively on unvisited successors
00032       visited[target(*oei, angelLCG)] = true;
00033       if (reachable (target(*oei, angelLCG), tgt, angelLCG)) return true;
00034     }
00035   } // end all outedges
00036 
00037   return false;
00038 } // end
00039 
00040 void vertex_upset (const c_graph_t::vertex_t v,
00041                    const c_graph_t& angelLCG,
00042                    vertex_set_t& upset) {
00043   upset.clear();
00044   upset.insert(v);
00045 
00046   c_graph_t::oei_t oei, oe_end;
00047   for (tie (oei, oe_end) = out_edges (v, angelLCG); oei != oe_end; ++oei) {
00048     vertex_set_t successor_upset;
00049     vertex_upset (target(*oei, angelLCG), angelLCG, successor_upset);
00050     upset.insert(successor_upset.begin(), successor_upset.end());
00051   }
00052 }
00053 
00054 void vertex_downset (const c_graph_t::vertex_t v,
00055                      const c_graph_t& angelLCG,
00056                      vertex_set_t& downset) {
00057   downset.clear();
00058   downset.insert(v);
00059 
00060   c_graph_t::iei_t iei, ie_end;
00061   for (tie (iei, ie_end) = in_edges (v, angelLCG); iei != ie_end; ++iei) {
00062     vertex_set_t predecessor_downset;
00063     vertex_downset (source (*iei, angelLCG), angelLCG, predecessor_downset);
00064     downset.insert(predecessor_downset.begin(), predecessor_downset.end());
00065   }
00066 }
00067 
00068 c_graph_t::edge_t getEdge(unsigned int i,
00069                           unsigned int j,
00070                           const c_graph_t& angelLCG) {
00071   c_graph_t::edge_t e;
00072   bool found_e;
00073   tie (e, found_e) = edge(i, j, angelLCG);
00074   THROW_EXCEPT_MACRO(!found_e, consistency_exception, "getEdge: edge could not be found from src,tgt");
00075   return e;
00076 } // end getEdge()
00077 
00078 bool lex_less_face (line_graph_t::face_t e1, line_graph_t::face_t e2,
00079                     const line_graph_t& lg) {
00080 
00081   line_graph_t::edge_t s1= source (e1, lg), s2= source (e2, lg),
00082                        t1= target (e1, lg), t2= target (e2, lg);
00083 
00084   // vertices i, j and k of e1 and e2
00085   property_map<pure_line_graph_t, vertex_name_t>::const_type evn= get(vertex_name, lg);
00086   c_graph_t::vertex_t vi1= evn[s1].first, vj1= evn[s1].second, vk1= evn[t1].second,
00087                        vi2= evn[s2].first, vj2= evn[s2].second, vk2= evn[t2].second;
00088 
00089 #ifndef NDEBUG
00090   if (evn[s1].second != evn[t1].first || evn[s2].second != evn[t2].first) {
00091     cout << "face " << s1 << "," << t1 << " or " << s2 << "," << t2 << " mismatch\n";
00092     write_graph ("graph", lg);
00093     write_vertex_property (std::cout, "edge names ", evn, lg); }
00094 #endif
00095 
00096   return vj1 < vj2 || (vj1 == vj2 && vi1 < vi2) || (vj1 == vj2 && vi1 == vi2 && vk1 <= vk2);
00097 }
00098 
00099 // =====================================================
00100 // vertex properties
00101 // =====================================================
00102 
00103 void in_out_path_lengths (const c_graph_t& cg,
00104                           std::vector<int>& vni, std::vector<int>& vli, 
00105                           std::vector<int>& vno, std::vector<int>& vlo) {
00106 
00107   vni.resize (num_vertices (cg));
00108   vli.resize (num_vertices (cg));
00109   vno.resize (num_vertices (cg));
00110   vlo.resize (num_vertices (cg));
00111 
00112   graph_traits<c_graph_t>::vertex_iterator     vi, v_end;
00113   int c;
00114 
00115   // initialize independent vertices and set vi to first intermediate vertex
00116   tie(vi, v_end)= vertices(cg);
00117   for (c= 0; c < cg.x(); c++, ++vi) {
00118     vni[c]= 1; vli[c]= 0; }
00119 
00120   // set other vertices from predecessors
00121   for (; vi != v_end; c++, ++vi) {
00122     vni[c]= 0; vli[c]= 0;
00123     graph_traits<c_graph_t>::in_edge_iterator iei, ie_end;
00124     for (tie(iei, ie_end)= in_edges(*vi, cg); iei != ie_end; ++iei) {
00125       c_graph_t::vertex_descriptor p= source (*iei, cg);
00126       vni[c]+= vni[p]; vli[c]+= vni[p] + vli[p]; }
00127   }
00128 
00129   // initialize vertices without successors
00130   tie(vi, v_end)= vertices(cg); --v_end; // to the last not behind
00131   c= num_vertices (cg) - 1;
00132   for (; out_degree (*v_end, cg) > 0; c--, --v_end) {
00133     vno[c]= 1; vlo[c]= 0; }
00134 
00135   // set other vertices from successors
00136   for (; c >= 0; c--, --v_end) {
00137     vno[c]= 0; vlo[c]= 0;
00138     graph_traits<c_graph_t>::out_edge_iterator oei, oe_end;
00139     for (tie(oei, oe_end)= out_edges(*v_end, cg); oei != oe_end; ++oei) {
00140       c_graph_t::vertex_descriptor s= target (*oei, cg);
00141       vno[c]+= vno[s]; vlo[c]+= vno[s] + vlo[s]; }
00142   }
00143 }
00144 
00145 void number_dependent_vertices (const c_graph_t& cg, std::vector<int>& v) {
00146   typedef c_graph_t::vertex_t                  vertex_t;
00147 
00148   v.resize (num_vertices (cg));
00149   fill (v.begin(), v.end(), 0);
00150 
00151   for (size_t c= 0; c < cg.dependents.size(); c++) {
00152     vertex_t dep= cg.dependents[c];
00153     // which vertices are relevant for dep ?
00154     queue<vertex_t> q; q.push (dep);
00155     vector<int>     dv (v.size()); dv[dep]= 1;
00156 
00157     while (!q.empty()) {
00158       vertex_t            v= q.front();
00159       c_graph_t::iei_t    iei, ie_end;
00160       for (tie(iei, ie_end)= in_edges(v, cg); iei != ie_end; ++iei) {
00161         vertex_t vin= source (*iei, cg);
00162         if (!dv[vin]) {
00163           dv[vin]= 1; q.push (vin); } }
00164       q.pop(); }
00165 
00166     for (size_t i= 0; i < v.size(); i++)
00167       v[i]+= dv[i];
00168   }
00169 }
00170 
00171 void number_independent_vertices (const c_graph_t& cg, std::vector<int>& v) {
00172 
00173   typedef c_graph_t::vertex_descriptor                  vertex_t;
00174   typedef graph_traits<c_graph_t>::vertex_iterator      vi_t;
00175   //  typedef graph_traits<c_graph_t>::in_edge_iterator     ie_t;
00176   typedef graph_traits<c_graph_t>::adjacency_iterator   ai_t;
00177 
00178   v.resize (num_vertices (cg));
00179   std::fill (v.begin(), v.end(), 0);
00180 
00181   vi_t                          ivi, iv_end;
00182   int                           c= 0;
00183   for (tie(ivi, iv_end)= vertices(cg), c= 0; c < cg.x() && ivi != iv_end; ++c, ++ivi) {
00184     // which vertices are reachable from *ivi ?
00185     std::queue<vertex_t> q; q.push (*ivi);
00186     std::vector<int>     dv (v.size()); dv[*ivi]= 1;
00187 
00188     while (!q.empty()) {
00189       vertex_t v= q.front();
00190       ai_t     ai, a_end;
00191       for (tie(ai, a_end)= adjacent_vertices(v, cg); ai != a_end; ++ai)
00192         if (!dv[*ai]) {
00193           dv[*ai]= 1; q.push (*ai); }
00194       q.pop(); 
00195     }
00196 
00197     for (size_t i= 0; i < v.size(); i++)
00198       v[i]+= dv[i];
00199   }
00200 }
00201 
00202 // =====================================================
00203 // graph transformations
00204 // =====================================================
00205 
00206 
00207 // independent vertices shall be first ones after permutation
00208 void permutate_vertices (const c_graph_t& gin, const vector<c_graph_t::vertex_t>& p,
00209                          c_graph_t& gout) {
00210   
00211   typedef c_graph_t::vertex_t vertex_t;
00212   // permutate vector of dependent vertices
00213   vector<vertex_t> deps;
00214   for (size_t c= 0; c < gin.dependents.size(); c++)
00215     deps.push_back (p[gin.dependents[c]]);
00216 
00217   c_graph_t gtmp (gin.v(), gin.x(), deps);
00218   int en= 0; // counter for edge_number
00219   
00220   c_graph_t::const_ew_t ewc= get(edge_weight, gin);  
00221   c_graph_t::ew_t       ew= get(edge_weight, gout);  
00222   c_graph_t::ei_t       ei, e_end;
00223   for (tie (ei, e_end)= edges (gin); ei != e_end; ++ei) {
00224     bool ins; c_graph_t::edge_t edge;
00225     tie (edge, ins)= add_edge (p[source(*ei,gin)], p[target(*ei,gin)], en++, gtmp);
00226     ew[edge]= ewc[*ei]; }
00227 
00228   // it may be usefull to sort the out_edges of each vertex 
00229 
00230   gtmp.next_edge_number= renumber_edges (gtmp);
00231   gout.swap (gtmp);
00232 }
00233 
00234 void independent_vertices_to_front (const c_graph_t& gin, 
00235                                     const vector<c_graph_t::vertex_t>& indeps,
00236                                     c_graph_t& gout) {
00237   typedef c_graph_t::vertex_t vertex_t;
00238   THROW_EXCEPT_MACRO (gin.x() != int (indeps.size()), consistency_exception, 
00239                    "Wrong number of independents");
00240   vector<vertex_t> counter (gin.v());
00241 
00242   int cindep= 0, cnotin= 0;
00243   for (int c= 0; c < gin.v(); c++)
00244     if (std::find (indeps.begin(), indeps.end(), vertex (c, gin)) != indeps.end())
00245       counter[c]= cindep++;
00246     else
00247       counter[c]= cnotin++;
00248 
00249   for (int c= 0; c < gin.v(); c++)
00250     if (std::find (indeps.begin(), indeps.end(), vertex (c, gin)) == indeps.end())
00251       counter[c]+= cindep;
00252 
00253   permutate_vertices (gin, counter, gout);
00254 }
00255 
00256 int renumber_edges (c_graph_t& cg) {
00257 
00258   graph_traits<c_graph_t>::edge_iterator      ei, e_end;
00259   property_map<c_graph_t, edge_index_t>::type  eid= get (edge_index, cg);
00260 
00261   int       c;
00262   for (tie (ei, e_end)= edges (cg), c= 0; ei != e_end; ++ei, c++)
00263     eid[*ei]= c;
00264   return c;
00265 }
00266 
00267 void take_over_successors (c_graph_t::vertex_t v1, c_graph_t::vertex_t v2, 
00268                            int offset, const c_graph_t& g1, 
00269                            int& edge_number, c_graph_t& g2) {
00270 
00271   c_graph_t::oei_t oei, oe_end;
00272   for (tie (oei, oe_end)= out_edges (v1, g1); oei != oe_end; ++oei) {
00273     int t= int (target (*oei, g1)), it= offset + t;
00274     THROW_DEBUG_EXCEPT_MACRO (it <= 0 || it >= g2.v(), consistency_exception, 
00275                            "Target out of range");
00276     c_graph_t::vertex_t vt= vertex (it, g2); // equivalent vertex in g2
00277     add_edge (v2, vt, edge_number++, g2); 
00278   }
00279 }
00280 
00281 
00282 
00283 // -----------------------------------------------------
00284 // some preprocessing removals
00285 // -----------------------------------------------------
00286 
00287 int remove_hoisting_vertices (c_graph_t& cg) {
00288 
00289   int hv= 0;
00290   c_graph_t::vi_t   vi, v_end;
00291   for (tie (vi, v_end)= vertices (cg); vi != v_end; ++vi)
00292     if (cg.vertex_type (*vi) == intermediate 
00293         && in_degree (*vi, cg) == 1 
00294         && out_degree (*vi, cg) == 1)
00295         hv+= vertex_elimination (*vi, cg);
00296   return hv;
00297 }
00298 
00299 void remove_parallel_edges (c_graph_t& cg) {
00300 
00301   typedef std::pair<c_graph_t::edge_t,bool> eb_t;
00302 
00303   c_graph_t gtmp (cg.v(), cg.x(), cg.dependents);
00304   int edge_number= 0;
00305 
00306   property_map<c_graph_t, edge_weight_t>::type   ew1= get(edge_weight, cg),
00307                                                   ew2= get(edge_weight, gtmp);
00308 
00309   c_graph_t::vi_t      vi, v_end;
00310   for (boost::tie (vi, v_end)= vertices (cg); vi != v_end; vi++) {
00311     std::vector<c_graph_t::vertex_t> targets;
00312     c_graph_t::oei_t      oei, oe_end;
00313     for (boost::tie (oei, oe_end)= out_edges (*vi, cg); oei != oe_end; oei++) {
00314       c_graph_t::vertex_t t= target (*oei, cg);
00315       eb_t edge_double= edge (*vi, t, gtmp);         // already inserted ?
00316       if (edge_double.second)                       // if so
00317         ew2[edge_double.first]+= ew1[*oei];
00318       else {
00319         eb_t new_edge= add_edge (*vi, t, edge_number++, gtmp);
00320         ew2[new_edge.first]=  ew1[*oei];
00321         targets.push_back (t); } 
00322     } }
00323   
00324   cg.swap (gtmp);
00325 }
00326 
00327 void remove_trivial_edges (c_graph_t& cg) {
00328 
00329   c_graph_t::ew_t ew= get(edge_weight, cg);
00330   for (bool found= true; found; ) {
00331     // write_edge_property (std::cout, "edge weights ", ew, cg);
00332     found= false;
00333     c_graph_t::ei_t ei, e_end;
00334     for (tie (ei, e_end)= edges (cg); ei != e_end; ++ei) {
00335       c_graph_t::edge_t   e= *ei;
00336       if (ew[e] == 1) {
00337         int ee; // number of eliminated edges
00338         if (cg.vertex_type (source (e, cg)) == independent) 
00339           ee= front_edge_elimination (e, cg);
00340         else ee= back_edge_elimination (e, cg); 
00341         if (ee > 0) {found= true; break;} // after elimination iterators may be invalid
00342       } } }
00343 }
00344 
00345 // =====================================================
00346 // Functions for partial accumulation simulated annealing (scarcity exploitation)
00347 // =====================================================
00348 
00349 double gen_prob() {
00350   double a = rand();
00351   double b = rand();
00352   return (a > b) ? b/a
00353                  : a/b;
00354 } // end gen_prob()
00355 
00356 unsigned int chooseTarget_sa(std::vector<double>& deltaE) {
00357   #define ECONST 2.71828
00358 
00359   //write_vector("deltaE (before adjustment): ", deltaE);
00360 
00361   // prefer improvements, if there are any
00362   double best_improvement = 100;
00363   for(unsigned int i = 0; i < deltaE.size(); i++)
00364     if(best_improvement > deltaE[i])
00365       best_improvement = deltaE[i];
00366   if(best_improvement < 0) {
00367     //cout << "best_improvement of " << best_improvement << " was recognized" << endl;
00368     for (unsigned int i = 0; i < deltaE.size(); i++)
00369       if (deltaE[i] > -1)
00370         deltaE[i] = 0;
00371       else
00372         deltaE[i] = pow(ECONST, (-deltaE[i]));
00373   }
00374   else {
00375     for (unsigned int i = 0; i < deltaE.size(); i++)
00376       deltaE[i] = 1.0/pow(ECONST,deltaE[i] + 1);
00377   }
00378 
00379   //write_vector("deltaE (after adjustment): ", deltaE);
00380 
00381   // normalize the probabilities
00382   double deltasum = 0;
00383   for(unsigned int i = 0; i < deltaE.size(); i++)
00384     deltasum += deltaE[i];
00385   //cout << "deltasum = " << deltasum << endl;
00386   for(unsigned int i = 0; i < deltaE.size(); i++)
00387     deltaE[i] = deltaE[i]/deltasum;
00388   //write_vector("deltaE (normalized): ", deltaE);
00389 
00390   // choose a vector index randomly
00391   double Pr = gen_prob();
00392   double current_ptr = deltaE[0];
00393   unsigned int choice_index = 0;
00394   while(current_ptr < Pr) {
00395     choice_index++;
00396     current_ptr += deltaE[choice_index];
00397   }
00398   //cout << "got index " << choice_index << " with prob " << Pr <<  endl;
00399   return choice_index;
00400 } // end chooseTarget_sa()
00401 
00402 int chooseEdgeElimRandomly(std::vector<double>& deltaE) {
00403   if (gen_prob() > 0.95)
00404     return -1;
00405 
00406   #define ECONST 2.71828
00407   //write_vector("deltaE (before adjustment): ", deltaE);
00408   for (unsigned int i = 0; i < deltaE.size(); i++)
00409     deltaE[i] = 1.0/pow(ECONST,deltaE[i] + 1);
00410   //write_vector("deltaE (after adjustment): ", deltaE);
00411 
00412   // normalize the probabilities
00413   double deltasum = 0;
00414   for(unsigned int i = 0; i < deltaE.size(); i++)
00415     deltasum += deltaE[i];
00416   //cout << "deltasum = " << deltasum << endl;
00417   for(unsigned int i = 0; i < deltaE.size(); i++)
00418     deltaE[i] = deltaE[i]/deltasum;
00419   //write_vector("deltaE (normalized): ", deltaE);
00420 
00421   // choose a vector index randomly
00422   double Pr = gen_prob();
00423   double current_ptr = deltaE[0];
00424   unsigned int choice_index = 0;
00425   while(current_ptr < Pr) {
00426     choice_index++;
00427     current_ptr += deltaE[choice_index];
00428   }
00429   return choice_index;
00430 } // end chooseEdgeElimRandomly()
00431 
00432 int chooseEdgeElimRandomlyGreedy(std::vector<double>& deltaE) {
00433   #define ECONST 2.71828
00434   //write_vector("deltaE (before adjustment): ", deltaE);
00435 
00436   // select against those that arent best
00437   double best = deltaE[0];
00438   size_t best_index = 0;
00439   for (size_t i = 0; i < deltaE.size(); i++) {
00440     if (deltaE[i] < best) {
00441       best = deltaE[i];
00442       best_index = i;
00443     }
00444   }
00445   for (size_t i = 0; i < deltaE.size(); i++)
00446     if (deltaE[i] > best)
00447       deltaE[i] += 30000;
00448 
00449   for (unsigned int i = 0; i < deltaE.size(); i++)
00450     deltaE[i] = 1.0/pow(ECONST,deltaE[i] + 1);
00451   //write_vector("deltaE (after adjustment): ", deltaE);
00452 
00453   // normalize the probabilities
00454   double deltasum = 0;
00455   for(unsigned int i = 0; i < deltaE.size(); i++)
00456     deltasum += deltaE[i];
00457   //cout << "deltasum = " << deltasum << endl;
00458   for(unsigned int i = 0; i < deltaE.size(); i++)
00459     deltaE[i] = deltaE[i]/deltasum;
00460   //write_vector("deltaE (normalized): ", deltaE);
00461 
00462   // choose a vector index randomly
00463   double Pr = gen_prob();
00464   double current_ptr = deltaE[0];
00465   unsigned int choice_index = 0;
00466   while(current_ptr < Pr) {
00467     choice_index++;
00468     current_ptr += deltaE[choice_index];
00469   }
00470 
00471   return choice_index;
00472 } // end chooseEdgeElimRandomly()
00473 
00474 } // namespace angel
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines