angel
mercurial changeset:
|
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