// -*- mode: C++; c-file-style: "cc-mode" -*- //************************************************************************* // DESCRIPTION: Verilator: Implementation of Christofides' algorithm to // approximate the solution to the traveling salesman problem. // // ISSUES: This isn't exactly Christofides algorithm; see the TODO // in perfectMatching(). True minimum-weight perfect matching // would produce a better result. How much better is TBD. // // Code available from: http://www.veripool.org/verilator // //************************************************************************* // // This program is free software; you can redistribute it and/or modify // it under the terms of either the GNU Lesser General Public License // Version 3 or the Perl Artistic License Version 2.0. // // Verilator is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU General Public License for more details. // //************************************************************************* #include "config_build.h" #include "verilatedos.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "V3Error.h" #include "V3Global.h" #include "V3File.h" #include "V3Graph.h" #include "V3TSP.h" #include VL_INCLUDE_UNORDERED_SET #include VL_INCLUDE_UNORDERED_MAP //###################################################################### // Support classes namespace V3TSP { static unsigned edgeIdNext = 0; static void selfTestStates(); static void selfTestString(); VL_DEBUG_FUNC; // Declare debug() } // namespace V3TSP // Vertex that tracks a per-vertex key template class TspVertexTmpl : public V3GraphVertex { private: T_Key m_key; public: TspVertexTmpl(V3Graph* graphp, const T_Key& k) : V3GraphVertex(graphp), m_key(k) {} virtual ~TspVertexTmpl() {} const T_Key& key() const { return m_key; } private: VL_UNCOPYABLE(TspVertexTmpl); }; // TspGraphTmpl represents a complete graph, templatized to work with // different T_Key types. template class TspGraphTmpl : public V3Graph { public: // TYPES typedef TspVertexTmpl Vertex; // MEMBERS typedef vl_unordered_map VMap; VMap m_vertices; // T_Key to Vertex lookup map // CONSTRUCTORS TspGraphTmpl() : V3Graph() {} virtual ~TspGraphTmpl() {} // METHODS void addVertex(const T_Key &key) { typename VMap::iterator itr = m_vertices.find(key); UASSERT(itr == m_vertices.end(), "Vertex already exists with same key"); Vertex *v = new Vertex(this, key); m_vertices[key] = v; } // For purposes of TSP, we are using non-directional graphs. // Map that onto the normally-directional V3Graph by creating // a matched pairs of opposite-directional edges to represent // each non-directional edge: void addEdge(const T_Key& from, const T_Key& to, int cost) { UASSERT(from != to, "Adding edge would form a loop"); Vertex* fp = findVertex(from); Vertex* tp = findVertex(to); // No need to dedup edges. // The only time we may create duplicate edges is when // combining the MST with the perfect-matched pairs, // and in that case, we want to permit duplicate edges. unsigned edgeId = ++V3TSP::edgeIdNext; // Record the 'id' which identifies a single bidir edge // in the user field of each V3GraphEdge: (new V3GraphEdge(this, fp, tp, cost))->user(edgeId); (new V3GraphEdge(this, tp, fp, cost))->user(edgeId); } bool empty() const { return m_vertices.empty(); } std::list keysToVertexList(const std::vector& odds) { std::list vertices; for(unsigned i = 0; i < odds.size(); ++i) { vertices.push_back(findVertex(odds.at(i))); } return vertices; } class EdgeCmp { // Provides a deterministic compare for outgoing V3GraphEdge's // to be used in Prim's algorithm below. Also used in the // perfectMatching() routine. public: // CONSTRUCTORS EdgeCmp() {} // METHODS bool operator() (const V3GraphEdge* ap, const V3GraphEdge* bp) { int aCost = ap->weight(); int bCost = bp->weight(); // Sort first on cost, lowest cost edges first: if (aCost < bCost) return true; if (bCost < aCost) return false; // Costs are equal. Compare edgeId's which should be unique. return ap->user() < bp->user(); } private: VL_UNCOPYABLE(EdgeCmp); }; static Vertex* castVertexp(V3GraphVertex* vxp) { return dynamic_cast(vxp); } // From *this, populate *mstp with the minimum spanning tree. // *mstp must be initially empty. void makeMinSpanningTree(TspGraphTmpl* mstp) { UASSERT(mstp->empty(), "Output graph must start empty"); // Use Prim's algorithm to efficiently construct the MST. vl_unordered_set visited_set; EdgeCmp cmp; typedef std::set PendingEdgeSet; // This is the set of pending edges from visited to unvisited // nodes. PendingEdgeSet pendingEdges(cmp); vluint32_t vert_ct = 0; for (V3GraphVertex* vxp = verticesBeginp(); vxp; vxp = vxp->verticesNextp()) { mstp->addVertex(castVertexp(vxp)->key()); vert_ct++; } // Choose an arbitrary start vertex and visit it; // all incident edges from this vertex go into a pending edge set. Vertex* start_vertexp = castVertexp(verticesBeginp()); visited_set.insert(start_vertexp); for (V3GraphEdge* edgep = start_vertexp->outBeginp(); edgep; edgep = edgep->outNextp()) { pendingEdges.insert(edgep); } // Repeatedly find the least costly edge in the pending set. // If it connects to an unvisited node, visit that node and update // the pending edge set. If it connects to an already visited node, // discard it and repeat again. unsigned edges_made = 0; while (!pendingEdges.empty()) { typename PendingEdgeSet::iterator firstIt = pendingEdges.begin(); V3GraphEdge* bestEdgep = *firstIt; pendingEdges.erase(firstIt); // bestEdgep->fromp() should be already seen Vertex* from_vertexp = castVertexp(bestEdgep->fromp()); UASSERT(visited_set.find(from_vertexp) != visited_set.end(), "Can't find vertex"); // If the neighbor is not yet visited, visit it and add its edges // to the pending set. Vertex* neighborp = castVertexp(bestEdgep->top()); if (visited_set.find(neighborp) == visited_set.end()) { int bestCost = bestEdgep->weight(); UINFO(6, "bestCost = "<key() <<" to "<key()<addEdge(from_vertexp->key(), neighborp->key(), bestCost); edges_made++; // Mark this vertex as visited visited_set.insert(neighborp); // Update the pending edges with new edges for (V3GraphEdge* edgep = neighborp->outBeginp(); edgep; edgep = edgep->outNextp()) { pendingEdges.insert(edgep); } } else { UINFO(6, "Discarding edge to already-visited neighbor " <key()<& oddKeys, TspGraphTmpl* outp) { UASSERT(outp->empty(), "Output graph must start empty"); std::list odds = keysToVertexList(oddKeys); vl_unordered_set unmatchedOdds; typedef typename std::list::iterator VertexListIt; for (VertexListIt it = odds.begin(); it != odds.end(); ++it) { outp->addVertex((*it)->key()); unmatchedOdds.insert(*it); } UASSERT(odds.size() % 2 == 0, "number of odd-order nodes should be even"); // TODO: The true Chrisofides algorithm calls for minimum-weight // perfect matching. Instead, we have a simple greedy algorithm // which might get close to the minimum, maybe, with luck? // // TODO: Revisit this. It's possible to compute the true minimum in // N*N*log(N) time using variants of the Blossom algorithm. // Implementing Blossom looks hard, maybe we can use an existing // open source implementation -- for example the "LEMON" library // which has a TSP solver. // ----- // Reuse the comparator from Prim's routine. The logic is the same // here. Note that the two V3GraphEdge's representing a single // bidir edge will collide in the pendingEdges set here, but this // is OK, we'll ignore the direction on the edge anyway. EdgeCmp cmp; typedef std::set PendingEdgeSet; PendingEdgeSet pendingEdges(cmp); for (VertexListIt it = odds.begin(); it != odds.end(); ++it) { for (V3GraphEdge* edgep = (*it)->outBeginp(); edgep; edgep = edgep->outNextp()) { pendingEdges.insert(edgep); } } // Iterate over all edges, in order from low to high cost. // For any edge whose ends are both odd-order vertices which // haven't been matched yet, match them. for (typename PendingEdgeSet::iterator it = pendingEdges.begin(); it != pendingEdges.end(); ++it) { Vertex* fromp = castVertexp((*it)->fromp()); Vertex* top = castVertexp((*it)->top()); if ((unmatchedOdds.find(fromp) != unmatchedOdds.end()) && (unmatchedOdds.find(top) != unmatchedOdds.end())) { outp->addEdge(fromp->key(), top->key(), (*it)->weight()); unmatchedOdds.erase(fromp); unmatchedOdds.erase(top); } } UASSERT(unmatchedOdds.empty(), "Algorithm should have processed all vertices"); } void combineGraph(const TspGraphTmpl& g) { vl_unordered_set edges_done; for (V3GraphVertex* vxp = g.verticesBeginp(); vxp; vxp = vxp->verticesNextp()) { Vertex* fromp = castVertexp(vxp); for (V3GraphEdge* edgep = fromp->outBeginp(); edgep; edgep = edgep->outNextp()) { Vertex* top = castVertexp(edgep->top()); if (edges_done.find(edgep->user()) == edges_done.end()) { addEdge(fromp->key(), top->key(), edgep->weight()); edges_done.insert(edgep->user()); } } } } void findEulerTourRecurse(vl_unordered_set* markedEdgesp, Vertex* startp, std::vector* sortedOutp) { Vertex* cur_vertexp = startp; // Go on a random tour. Fun! std::vector tour; do { UINFO(6, "Adding "<key()<<" to tour.\n"); tour.push_back(cur_vertexp); // Look for an arbitrary edge we've not yet marked for (V3GraphEdge* edgep = cur_vertexp->outBeginp(); edgep; edgep = edgep->outNextp()) { vluint32_t edgeId = edgep->user(); if (markedEdgesp->end() == markedEdgesp->find(edgeId)) { // This edge is not yet marked, so follow it. markedEdgesp->insert(edgeId); Vertex* neighborp = castVertexp(edgep->top()); UINFO(6, "following edge "<key() <<" to "<key()<key()<::iterator it = tour.begin(); it != tour.end(); ++it) { Vertex* vxp = *it; bool recursed; do { recursed = false; // Look for an arbitrary edge at vxp we've not yet marked for (V3GraphEdge* edgep = vxp->outBeginp(); edgep; edgep = edgep->outNextp()) { vluint32_t edgeId = edgep->user(); if (markedEdgesp->end() == markedEdgesp->find(edgeId)) { UINFO(6, "Recursing.\n"); findEulerTourRecurse(markedEdgesp, vxp, sortedOutp); recursed = true; goto recursed; } } recursed: ; } while (recursed); sortedOutp->push_back(vxp->key()); } UINFO(6, "Tour was: "); for (typename std::vector::iterator it = tour.begin(); it != tour.end(); ++it) { Vertex* vxp = *it; UINFONL(6, " "<key()); } UINFONL(6, "\n"); } void dumpGraph(std::ostream& os, const string& nameComment) const { // UINFO(0) as controlled by caller os<<"At "<verticesNextp()) { Vertex* tspvp = castVertexp(vxp); os<<" "<key()<outBeginp(); edgep; edgep = edgep->outNextp()) { Vertex* neighborp = castVertexp(edgep->top()); os<<" has edge "<user()<<" to "<key()< logp(V3File::new_ofstream(filename)); if (logp->fail()) v3fatal("Can't write "<* sortedOutp) { UASSERT(sortedOutp->empty(), "Output graph must start empty"); if (debug() >= 6) dumpDotFilePrefixed("findEulerTour"); vl_unordered_set markedEdges; // Pick a start node Vertex* start_vertexp = castVertexp(verticesBeginp()); findEulerTourRecurse(&markedEdges, start_vertexp, sortedOutp); } std::vector getOddDegreeKeys() const { std::vector result; for (V3GraphVertex* vxp = verticesBeginp(); vxp; vxp = vxp->verticesNextp()) { Vertex* tspvp = castVertexp(vxp); vluint32_t degree = 0; for (V3GraphEdge* edgep = vxp->outBeginp(); edgep; edgep = edgep->outNextp()) { degree++; } if (degree & 1) { result.push_back(tspvp->key()); } } return result; } private: Vertex* findVertex(const T_Key& key) const { typename VMap::const_iterator it = m_vertices.find(key); UASSERT(it != m_vertices.end(), "Vertex not found"); return it->second; } VL_UNCOPYABLE(TspGraphTmpl); }; //###################################################################### // Main algorithm void V3TSP::tspSort(const V3TSP::StateVec& states, V3TSP::StateVec* resultp) { UASSERT(resultp->empty(), "Output graph must start empty"); // Make this TSP implementation work for graphs of size 0 or 1 // which, unfortunately, is a special case as the following // code assumes >= 2 nodes. if (states.size() == 0) { return; } if (states.size() == 1) { resultp->push_back(*(states.begin())); return; } // Build the initial graph from the starting state set. typedef TspGraphTmpl Graph; Graph graph; for (V3TSP::StateVec::const_iterator it = states.begin(); it != states.end(); ++it) { graph.addVertex(*it); } for (V3TSP::StateVec::const_iterator it = states.begin(); it != states.end(); ++it) { for (V3TSP::StateVec::const_iterator jt = it; jt != states.end(); ++jt) { if (it == jt) continue; graph.addEdge(*it, *jt, (*it)->cost(*jt)); } } // Create the minimum spanning tree Graph minGraph; graph.makeMinSpanningTree(&minGraph); if (debug() >= 6) minGraph.dumpGraphFilePrefixed("minGraph"); std::vector oddDegree = minGraph.getOddDegreeKeys(); Graph matching; graph.perfectMatching(oddDegree, &matching); if (debug() >= 6) matching.dumpGraphFilePrefixed("matching"); // Adds edges to minGraph, the resulting graph will have even number of // edge counts at every vertex: minGraph.combineGraph(matching); V3TSP::StateVec prelim_result; minGraph.findEulerTour(&prelim_result); UASSERT(prelim_result.size() >= states.size(), "Algorithm size error"); // Discard duplicate nodes that the Euler tour might contain. { vl_unordered_set seen; for (V3TSP::StateVec::iterator it = prelim_result.begin(); it != prelim_result.end(); ++it) { const TspStateBase* elemp = *it; if (seen.find(elemp) == seen.end()) { seen.insert(elemp); resultp->push_back(elemp); } } } UASSERT(resultp->size() == states.size(), "Algorithm size error"); // Find the most expensive arc and rotate the list so that the most // expensive arc connects the last and first elements. (Since we're not // modeling something that actually cycles back, we don't need to pay // that cost at all.) { unsigned max_cost = 0; unsigned max_cost_idx = 0; for (unsigned i = 0; i < resultp->size(); ++i) { const TspStateBase* ap = (*resultp)[i]; const TspStateBase* bp = (i+1 == resultp->size()) ? (*resultp)[0] : (*resultp)[i+1]; unsigned cost = ap->cost(bp); if (cost > max_cost) { max_cost = cost; max_cost_idx = i; } } if (max_cost_idx == resultp->size() - 1) { // List is already rotated for minimum cost. stop. UASSERT(resultp->size() == resultp->size(), "sizes should match"); return; } V3TSP::StateVec new_result; unsigned i = max_cost_idx + 1; UASSERT(i < resultp->size(), "Algorithm size error"); while(i != max_cost_idx) { new_result.push_back((*resultp)[i]); i++; if (i >= resultp->size()) { i = 0; } } new_result.push_back((*resultp)[i]); UASSERT(resultp->size() == new_result.size(), "Algorithm size error"); *resultp = new_result; } } //###################################################################### // Self Tests class TspTestState : public V3TSP::TspStateBase { public: TspTestState(unsigned xpos, unsigned ypos) : m_xpos(xpos), m_ypos(ypos), m_serial(++m_serialNext) {} ~TspTestState() {} virtual int cost(const TspStateBase* otherp) const { return cost(dynamic_cast(otherp)); } static unsigned diff(unsigned a, unsigned b) { if (a>b) return a-b; return b-a; } virtual int cost(const TspTestState* otherp) const { // For test purposes, each TspTestState is merely a point // on the cartesian plane; cost is the linear distance // between two points. unsigned xabs, yabs; xabs = diff(otherp->m_xpos, m_xpos); yabs = diff(otherp->m_ypos, m_ypos); return (int)(0.5 + sqrt(xabs*xabs + yabs*yabs)); } unsigned xpos() const { return m_xpos; } unsigned ypos() const { return m_ypos; } bool operator< (const TspStateBase& other) const { return operator< (dynamic_cast(other)); } bool operator< (const TspTestState& other) const { return m_serial < other.m_serial; } private: unsigned m_xpos; unsigned m_ypos; unsigned m_serial; static unsigned m_serialNext; }; unsigned TspTestState::m_serialNext = 0; void V3TSP::selfTestStates() { // Linear test -- coords all along the x-axis { V3TSP::StateVec states; TspTestState s10(10,0); TspTestState s60(60,0); TspTestState s20(20,0); TspTestState s100(100,0); TspTestState s5(5,0); states.push_back(&s10); states.push_back(&s60); states.push_back(&s20); states.push_back(&s100); states.push_back(&s5); V3TSP::StateVec result; tspSort(states, &result); V3TSP::StateVec expect; expect.push_back(&s100); expect.push_back(&s60); expect.push_back(&s20); expect.push_back(&s10); expect.push_back(&s5); if (expect != result) { for (V3TSP::StateVec::iterator it = result.begin(); it != result.end(); ++it) { const TspTestState* statep = dynamic_cast(*it); cout<xpos()<<" "; } cout<(*it); cout<xpos()<<","<ypos()<<" "; } cout< Graph; Graph graph; graph.addVertex("0"); graph.addVertex("1"); graph.addVertex("2"); graph.addVertex("3"); graph.addEdge("0", "1", 3943); graph.addEdge("0", "2", 3456); graph.addEdge("0", "3", 4920); graph.addEdge("1", "2", 2730); graph.addEdge("1", "3", 8199); graph.addEdge("2", "3", 4130); Graph minGraph; graph.makeMinSpanningTree(&minGraph); if (debug() >= 6) minGraph.dumpGraphFilePrefixed("minGraph"); std::vector oddDegree = minGraph.getOddDegreeKeys(); Graph matching; graph.perfectMatching(oddDegree, &matching); if (debug() >= 6) matching.dumpGraphFilePrefixed("matching"); minGraph.combineGraph(matching); std::vector result; minGraph.findEulerTour(&result); std::vector expect; expect.push_back("0"); expect.push_back("2"); expect.push_back("1"); expect.push_back("2"); expect.push_back("3"); if (expect != result) { for (std::vector::iterator it = result.begin(); it != result.end(); ++it) { cout<<*it<<" "; } cout<