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

Generated on Wed Mar 11 10:33:11 2009 for angel by  doxygen 1.5.3