Grafalgo
Library of useful data structures and algorithms
|
00001 #include "edmonds.h" 00002 00003 using namespace grafalgo; 00004 00005 // Find a maximum size matching in the graph graf and 00006 // return it as a list of edges. 00007 edmonds::edmonds(Graph& graf1, Dlist& match1, int &size) 00008 : graf(&graf1), match(&match1) { 00009 vertex u, v; edge e; 00010 blossoms = new Partition(graf->n()); // set per blossom 00011 augpath = new RlistSet(graf->m()); // reversible list 00012 origin = new vertex[graf->n()+1]; // original vertex for each blossom 00013 bridge = new BridgePair[graf->n()+1];// edge that formed a blossom 00014 state = new stype[graf->n()+1]; // state used in path search 00015 pEdge = new edge[graf->n()+1]; // edge to parent in tree 00016 mEdge = new edge[graf->n()+1]; // incident matching edge (if any) 00017 mark = new bool[graf->n()+1]; // mark bits used by nca 00018 00019 mSize = stepCount = blossomCount = pathInitTime = pathFindTime = 0; 00020 int t1, t2, t3; 00021 t1 = Util::getTime(); 00022 // Create initial maximal (not maximum) matching 00023 for (u = 1; u <= graf->n(); u++) { 00024 mEdge[u] = 0; mark[u] = false; 00025 } 00026 match->clear(); size = 0; 00027 for (e = graf->first(); e != 0; e = graf->next(e)) { 00028 u = graf->left(e); v = graf->right(e); 00029 if (mEdge[u] == 0 && mEdge[v] == 0) { 00030 mEdge[u] = mEdge[v] = e; 00031 match->addLast(e); mSize++; 00032 } 00033 } 00034 iSize = mSize; 00035 00036 t2 = Util::getTime(); 00037 imatchTime = (t2-t1); 00038 while((e = findpath()) != 0) { augment(e); mSize++; } 00039 size = mSize; 00040 t3 = Util::getTime(); 00041 rmatchTime = (t3-t2); 00042 00043 delete blossoms; delete augpath; delete [] origin; 00044 delete [] bridge; delete [] pEdge; delete [] mEdge; delete[] mark; 00045 } 00046 00047 void edmonds::augment(edge e) { 00048 // Modify the matching by augmenting along the path defined by 00049 // the list in the augpath data structure whose last element is e. 00050 00051 vertex u, v; edge e1; 00052 while (true) { 00053 e1 = augpath->first(e); 00054 if (match->member(e1)) match->remove(e1); 00055 else { 00056 match->addLast(e1); 00057 mEdge[graf->left(e1)] = mEdge[graf->right(e1)] = e1; 00058 } 00059 if (e == augpath->first(e)) break; 00060 e = augpath->pop(e); 00061 } 00062 } 00063 00064 // If u and v are in the same tree, return their nearest common 00065 // ancestor in the current "condensed graph". To avoid excessive 00066 // search time, search upwards from both vertices in parallel, using 00067 // mark bits to identify the nca. Before returning, clear the mark 00068 // bits by traversing the paths a second time. The mark bits are 00069 // initialized in the constructor. 00070 vertex edmonds::nca(vertex u, vertex v) { 00071 vertex x,px,y,py,result; 00072 00073 // first pass to find the nca 00074 x = u; px = (pEdge[x] != 0 ? graf->mate(x,pEdge[x]) : 0); 00075 y = v; py = (pEdge[y] != 0 ? graf->mate(y,pEdge[y]) : 0); 00076 while (true) { 00077 if (x == y) { result = x; break; } 00078 if (px == 0 && py == 0) { result = 0; break; } 00079 if (px != 0) { 00080 if (mark[x]) { result = x; break; } 00081 mark[x] = true; 00082 x = origin[blossoms->find(px)]; 00083 px = (pEdge[x] != 0 ? graf->mate(x,pEdge[x]) : 0); 00084 } 00085 if (py != 0) { 00086 if (mark[y]) { result = y; break; } 00087 mark[y] = true; 00088 y = origin[blossoms->find(py)]; 00089 py = (pEdge[y] != 0 ? graf->mate(y,pEdge[y]) : 0); 00090 } 00091 } 00092 // second pass to clear mark bits 00093 x = u, y = v; 00094 while (mark[x] || mark[y]) { 00095 mark[x] = mark[y] = false; 00096 px = (pEdge[x] != 0 ? graf->mate(x,pEdge[x]) : 0); 00097 py = (pEdge[y] != 0 ? graf->mate(y,pEdge[y]) : 0); 00098 x = (px == 0 ? x : origin[blossoms->find(px)]); 00099 y = (py == 0 ? y : origin[blossoms->find(py)]); 00100 } 00101 return result; 00102 } 00103 00104 edge edmonds::path(vertex a, vertex b) { 00105 // Find path joining a and b defined by parent pointers and bridges. 00106 // A is assumed to be a descendant of b, and the path to b from a 00107 // is assumed to pass through the matching edge incident to a. 00108 // Return path in the augpath data structure. 00109 vertex pa, p2a, da; edge e, e1, e2; 00110 if (a == b) return 0; 00111 if (state[a] == even) { 00112 e1 = pEdge[a]; 00113 pa = graf->mate(a,e1); 00114 if (pa == b) return e1; 00115 e2 = pEdge[pa]; 00116 p2a = graf->mate(pa,e2); 00117 e = augpath->join(e1,e2); 00118 e = augpath->join(e,path(p2a,b)); 00119 return e; 00120 } else { 00121 e = bridge[a].e; da = bridge[a].v; 00122 e = augpath->join(augpath->reverse(path(da,a)),e); 00123 e = augpath->join(e,path(graf->mate(da,e),b)); 00124 return e; 00125 } 00126 } 00127 00128 // Search for an augmenting path in the graph graf. 00129 // On success, the path data structure will include a list 00130 // that forms the augmenting path. In this case, the last 00131 // edge in the list is returned as the function value. 00132 // On failure, returns 0. 00133 edge edmonds::findpath() { 00134 vertex u,v,vp,w,wp,x,y; edge e, f; 00135 00136 int t1, t2, t3; 00137 t1 = Util::getTime(); 00138 00139 blossoms->clear(); 00140 for (u = 1; u <= graf->n(); u++) { 00141 state[u] = even; pEdge[u] = 0; origin[u] = u; 00142 } 00143 00144 for (e = match->first(); e != 0; e = match->next(e)) { 00145 u = graf->left(e); v = graf->right(e); 00146 state[u] = state[v] = unreached; 00147 } 00148 List q(graf->m()); // list of edges to be processed in main loop 00149 for (e = 1; e <= graf->m(); e++) { 00150 if (state[graf->left(e)] == even || 00151 state[graf->right(e)] == even) 00152 q.addLast(e); 00153 } 00154 00155 t2 = Util::getTime(); 00156 pathInitTime += (t2-t1); 00157 while (!q.empty()) { 00158 stepCount++; 00159 e = q.first(); q.removeFirst(); 00160 v = graf->left(e); vp = origin[blossoms->find(v)]; 00161 if (state[vp] != even) { 00162 v = graf->right(e); vp = origin[blossoms->find(v)]; 00163 } 00164 w = graf->mate(v,e); wp = origin[blossoms->find(w)]; 00165 if (vp == wp) continue; // skip internal edges in a blossom 00166 if (state[wp] == unreached) { 00167 // w is not contained in a blossom and is matched 00168 // so extend tree and add newly eligible edges to q 00169 x = graf->mate(w,mEdge[w]); 00170 state[w] = odd; pEdge[w] = e; 00171 state[x] = even; pEdge[x] = mEdge[w]; 00172 for (f = graf->firstAt(x); f != 0; 00173 f = graf->nextAt(x,f)) { 00174 if ((f != mEdge[x]) && !q.member(f)) 00175 q.addLast(f); 00176 } 00177 continue; 00178 } 00179 u = nca(vp,wp); 00180 if (state[wp] == even && u == 0) { 00181 // vp, wp are different trees - construct path & return 00182 x = vp; 00183 while (pEdge[x] != 0) { 00184 x = origin[blossoms->find( 00185 graf->mate(x,pEdge[x]))]; 00186 } 00187 y = wp; 00188 while (pEdge[y] != 0) { 00189 y = origin[blossoms->find( 00190 graf->mate(y,pEdge[y]))]; 00191 } 00192 e = augpath->join(augpath->reverse(path(v,x)),e); 00193 e = augpath->join(e,path(w,y)); 00194 t3 = Util::getTime(); 00195 pathFindTime += (t3-t2); 00196 return e; 00197 } else if (state[wp] == even) { 00198 // vp and wp are in same tree - collapse blossom 00199 blossomCount++; 00200 x = vp; 00201 int fu = blossoms->find(u); 00202 int fx = blossoms->find(x); 00203 while (fx != fu) { 00204 fu = blossoms->link(fx,fu); 00205 if (state[x] == odd) { 00206 bridge[x].e = e; bridge[x].v = v; 00207 for (f = graf->firstAt(x); f != 0; 00208 f = graf->nextAt(x,f)) 00209 if (!q.member(f)) q.addLast(f); 00210 } 00211 x = origin[blossoms->find( 00212 graf->mate(x,pEdge[x]))]; 00213 fx = blossoms->find(x); 00214 } 00215 y = wp; 00216 int fy = blossoms->find(y); 00217 while (fy != fu) { 00218 fu = blossoms->link(fy,fu); 00219 if (state[y] == odd) { 00220 bridge[y].e = e; bridge[y].v = w; 00221 for (f = graf->firstAt(y); f != 0; 00222 f = graf->nextAt(y,f)) 00223 if (!q.member(f)) q.addLast(f); 00224 } 00225 y = origin[blossoms->find( 00226 graf->mate(y,pEdge[y]))]; 00227 fy = blossoms->find(y); 00228 } 00229 origin[fu] = u; 00230 } 00231 } 00232 t3 = Util::getTime(); 00233 pathFindTime += (t3-t2); 00234 return 0; 00235 } 00236 00243 string& edmonds::statString(bool verbose, string& s) { 00244 string s1; 00245 if (verbose) { 00246 s = "iSize=" + Util::num2string(iSize,s1) + " "; 00247 s += "mSize=" + Util::num2string(mSize,s1) + " "; 00248 s += "stepCount=" + Util::num2string(stepCount,s1) + " "; 00249 s += "blossomCount=" + Util::num2string(blossomCount,s1) + " "; 00250 s += "imatchTime=" + Util::num2string(imatchTime,s1) + " "; 00251 s += "rmatchTime=" + Util::num2string(rmatchTime,s1) + " "; 00252 s += "pathInitTime=" + Util::num2string(pathInitTime,s1) + " "; 00253 s += "pathFindTime=" + Util::num2string(pathFindTime,s1); 00254 } else { 00255 s = Util::num2string(iSize,s1) + " "; 00256 s += Util::num2string(mSize,s1) + " "; 00257 s += Util::num2string(stepCount,s1) + " "; 00258 s += Util::num2string(blossomCount,s1) + " "; 00259 s += Util::num2string(imatchTime,s1) + " "; 00260 s += Util::num2string(rmatchTime,s1) + " "; 00261 s += Util::num2string(pathInitTime,s1) + " "; 00262 s += Util::num2string(pathFindTime,s1); 00263 } 00264 return s; 00265 }