本文原创
最短路问题乃图论经典问题之一,著名的算法有Dijkstra、Bellman-Ford、Floyd-Warshall等,而求第K短路径的最为经典的算法由JIN Y. YEN于1971年在其论文《Finding The K Shortest Loopless Paths In A Network》中提出,对于权值严格非负图,其算法的时间复杂度上限为K*N^3/2,空间复杂度为N^2 + KN。
在Ernesto Q.V. Martins和Marta M.B. Pascoal合著的论文《A new implementation of Yen’s ranking loopless paths algorithm》中,他们对Yen氏算法进行了优化改进,论文中的算法伪代码如下:
本文以邻接表的存储形式对上述算法进行了实现,实现中调用了Dijkstra算法,实现中还构建了反向图,相关源代码摘录如下,首先是结点和边以及路径解的存储形式:
class Arc { public: Arc() : ArcID(0), SourceID(0), DestinationID(0), cost(INF), nextArc(NULL) {} int ArcID; // int SourceID; // int DestinationID; // int cost; Arc *nextArc; }; class rArc // For Yen's kShortestPath { public: rArc() : ArcID(0), SourceID(0), DestinationID(0), cost(INF), nextArc(NULL) {} int ArcID; // int SourceID; // int DestinationID; // int cost; rArc *nextArc; }; class Node { public: Node() : ID(0), firstArc(NULL) {} int ID; vector<Arc*> previous; // For Yen's kShortestPath Arc *firstArc; }; class rNode // For Yen's kShortestPath { public: rNode() : ID(0), firstArc(NULL) {} int ID; vector<rArc*> previous; // For Yen's kShortestPath rArc *firstArc; }; class Path { public: Path(int d) : costSum(0), deviation(d), deviationFromWhichPath(0) {} Path& operator=(const Path& rhs) { nodePath.assign(rhs.nodePath.begin(), rhs.nodePath.end()); arcPath.assign(rhs.arcPath.begin(), rhs.arcPath.end()); costSum = rhs.costSum; deviation = rhs.deviation; deviationFromWhichPath = rhs.deviationFromWhichPath; return *this; } vector<int> nodePath; vector<int> arcPath; int costSum; int deviation; int deviationFromWhichPath; }; bool pathCmp(const Path* lhs, const Path* rhs) { if (lhs->costSum == rhs->costSum) { if (lhs->nodePath.size() == rhs->nodePath.size()) return lhs->deviationFromWhichPath < rhs->deviationFromWhichPath; return lhs->nodePath.size() < rhs->nodePath.size(); } return lhs->costSum < rhs->costSum; }
针对正向图和反向图的Dijkstra算法,及其路径的记录函数:
void dijkstra(Node *NodeList, int NodeNum, int SourceID) { int i, u, v; // set<int> S; priority_queue<pair<int, int>, vector<pair<int, int> >, cmp > priorityQ; for (i = 0; i < NodeNum; i++) { if (SourceID != i) { dist[i] = INF; pi[i] = NULL; } } dist[SourceID] = 0; pi[SourceID] = 0; arcRecord[SourceID] = 0; priorityQ.push(make_pair(SourceID, 0)); while (priorityQ.empty() != true) { u = priorityQ.top().first; priorityQ.pop(); // S.insert(u); Arc *parc; parc = NodeList[u].firstArc; while (parc != NULL) { v = parc->DestinationID; if (dist[v] > dist[u] + parc->cost) { dist[v] = dist[u] + parc->cost; pi[v] = u; arcRecord[v] = parc->ArcID; priorityQ.push(make_pair(v, dist[v])); } parc = parc->nextArc; } } } void rDijkstra(rNode *rNodeList, int *Pi, int NodeNum, int SourceID) { int i, u, v; // set<int> S; priority_queue<pair<int, int>, vector<pair<int, int> >, cmp > priorityQ; for (i = 0; i < NodeNum; i++) if (SourceID != i) { Pi[i] = INF; pi[i] = NULL; } Pi[SourceID] = 0; pi[SourceID] = 0; arcRecord[SourceID] = 0; priorityQ.push(make_pair(SourceID, 0)); while (priorityQ.empty() != true) { u = priorityQ.top().first; priorityQ.pop(); // S.insert(u); rArc *rparc; rparc = rNodeList[u].firstArc; while (rparc != NULL) { v = rparc->DestinationID; if (Pi[v] > Pi[u] + rparc->cost) { Pi[v] = Pi[u] + rparc->cost; pi[v] = u; arcRecord[v] = rparc->ArcID; priorityQ.push(make_pair(v, Pi[v])); } rparc = rparc->nextArc; } } } void recordDirectRouteForNode(vector<int>& route, int *path, int src, int dst, bool pop_back = true) { stack<int> S; while (dst != src) { S.push(dst); dst = path[dst]; } S.push(dst); while (S.empty() != true) { route.push_back(S.top()); S.pop(); } if (pop_back == true) route.pop_back(); } void recordDirectRouteForArc(vector<int>& route, int *path, int* pathHelper, int src, int dst) { stack<int> S; int temp; while (dst != src) { temp = dst; dst = path[dst]; S.push(dst); dst = pathHelper[temp]; } while (S.empty() != true) { route.push_back(S.top()); S.pop(); } }
然后是论文中提到的Correct labels of v_i^k successors using backward star form:
void correctLabels(Node *NodeList, int *Pi, int v_i_k, vector<bool>& isRemoved) { unsigned int i; queue<int> Q; Q.push(v_i_k); int x; do { x = Q.front(); Q.pop(); int priv; for (i = 0; i < NodeList[x].previous.size(); i++) { priv = NodeList[x].previous[i]->SourceID; if (isRemoved[priv] == true) // skip removed node to save time continue; if (Pi[priv] > Pi[x] + NodeList[x].previous[i]->cost) { Pi[priv] = Pi[x] + NodeList[x].previous[i]->cost; // pi_j <-- pi_i + c_{ji} Q.push(priv); // list <-- list U {j} } } } while (Q.empty() != true); } void correctLabelsReverseGraph(rNode *rNodeList, int *Pi, int v_i_k, vector<bool>& isRemoved) { queue<int> Q; Q.push(v_i_k); int x; do { x = Q.front(); Q.pop(); rArc *rparc = rNodeList[x].firstArc; while (rparc != NULL) { if (isRemoved[rparc->DestinationID] == true) // skip removed node to save time { rparc = rparc->nextArc; continue; } if (Pi[rparc->DestinationID] > Pi[x] + rparc->cost) { Pi[rparc->DestinationID] = Pi[x] + rparc->cost; // pi_j <-- pi_i + c_{ji} Q.push(rparc->DestinationID); // list <-- list U {j} } rparc = rparc->nextArc; } } while (Q.empty() != true); }
最后的重头戏来啦,直接可调用的函数,输入正向图和反向图以及源点和终点,路径就存放在传入的path中啦:
void kShortestPath(Path& path, Node *NodeList, rNode *rNodeList, int NodeNum, int src, int dst, map<int, int>& arcID2Cost, unsigned int K = 1) { dijkstra(NodeList, NodeNum, src); recordNodeRouteForYen(path, pi, src, dst); recordArcRouteForYen(path, arcRecord, pi, src, dst, arcID2Cost); if (K == 1) return; unsigned int i, j; unsigned int k = 1; int *Pi = new int[NodeNum]; vector<bool> isRemoved; vector<bool> is_in_d; for (i = 0; i < (unsigned int)NodeNum; i++) { isRemoved.push_back(false); is_in_d.push_back(false); } bool (*path_func_ptr)(const Path*, const Path*) = pathCmp; set<Path*, bool (*)(const Path*, const Path*)> X(path_func_ptr); X.insert(&path); set<int> *determinedStartingArcs = new set<int>[NodeNum]; vector<int> nodePathTemp; vector<int> nodePath_t_2_v_i_k, arcPath_t_2_v_i_k;; // store path Tt(v_i^k) vector<Path*> erasedPath; for (i = 0; i < K-1; i++) erasedPath.push_back(new Path(src)); Path path_k(src); int *d_p_k = new int[K]; int deviationFrom; // =============== While (X != empty and k < K) Do =============== while (X.empty() != true && k < K) { d_p_k[k-1] = (*X.begin())->deviation; // d_p_k[0] is src, for "d.insert(make_pair(&path, src));" path_k = **X.begin(); // p_k <-- shortest loopless path in X X.erase(X.begin()); // X <-- X - {p_k} *erasedPath[k-1] = path_k; // erasedPath[k-1] is the k Shortest Path path = path_k; // ========== Remove loopless path pk, except node t, from the network ========== Arc *parc = NULL, *qarc = NULL; rArc *rparc = NULL, *rqarc = NULL; for (i = 0; i < path_k.nodePath.size()-1; i++) // size()-1 due to except node t { int removing = path_k.nodePath[i]; // forward graph, out-degree parc = NodeList[removing].firstArc; while (parc != NULL) { parc->cost = INF; // equivalent to remove from the network parc = parc->nextArc; } // forward graph, in-degree for (j = 0; j < NodeList[removing].previous.size(); j++) { qarc = NodeList[removing].previous[j]; qarc->cost = INF; } // reverse graph, out-degree rparc = rNodeList[removing].firstArc; while (rparc != NULL) { rparc->cost = INF; // equivalent to remove from the network rparc = rparc->nextArc; } // reverse graph, in-degree for (j = 0; j < rNodeList[removing].previous.size(); j++) { rqarc = rNodeList[removing].previous[j]; rqarc->cost = INF; } isRemoved[removing] = true; } // ========== Remove arc(d(p_k), i), i in N, of p1, ..., p_{k-1} from the network ========== if (k != 1) { for (i = 0; i < k-1; i++) { for (j = 0; j < erasedPath[i]->nodePath.size() && j < path_k.nodePath.size(); j++) { if (erasedPath[i]->nodePath[j] != path_k.nodePath[j]) { deviationFrom = path_k.nodePath[j-1]; break; } } if (j != erasedPath[i]->nodePath.size() && j != path_k.nodePath.size()) determinedStartingArcs[deviationFrom].insert(erasedPath[i]->nodePath[j]); } } // ========== Tt <-- shortest tree rooted at t in the network ========== rDijkstra(rNodeList, Pi, NodeNum, dst); // ========== For (v_i^k in { v_{l_k-1}^k, ..., d(p_k)}) Do ========== nodePathTemp.assign(path_k.nodePath.begin(), path_k.nodePath.end()); int v_i_k, v_ip1_k; // v_ip1_k is v_{i+1}^k v_ip1_k = nodePathTemp.back(); nodePathTemp.pop_back(); do { v_i_k = nodePathTemp.back(); nodePathTemp.pop_back(); // ========== Restore node v_i^k in the network ========== // forward graph, out-degree parc = NodeList[v_i_k].firstArc; while (parc != NULL) { // when analysing d(p_k) the arcs starting in d(p_k) which belong to other loopless paths already determined should not be used. if (is_in_d[v_i_k] == true && determinedStartingArcs[v_i_k].count(parc->DestinationID) != 0) // this arc has belong to another loopless path { parc = parc->nextArc; continue; } if (parc->DestinationID != v_ip1_k && isRemoved[parc->DestinationID] == false) // arc (v_i^k, v_{i+1}^k) restore after { parc->cost = arcID2Cost.at(parc->ArcID); // equivalent to restore in the network if (Pi[v_i_k] > Pi[parc->DestinationID] + parc->cost) Pi[v_i_k] = Pi[parc->DestinationID] + parc->cost; // Calculate pi_{v_i^k} using forward star form } parc = parc->nextArc; } // forward graph, in-degree for (i = 0; i < NodeList[v_i_k].previous.size(); i++) { qarc = NodeList[v_i_k].previous[i]; if (isRemoved[qarc->SourceID] == false) qarc->cost = arcID2Cost.at(qarc->ArcID); // equivalent to restore in the network } // reverse graph, out-degree rparc = rNodeList[v_i_k].firstArc; while (rparc != NULL) { if (isRemoved[rparc->DestinationID] == false) rparc->cost = arcID2Cost.at(rparc->ArcID); rparc = rparc->nextArc; } // reverse graph, in-degree for (i = 0; i < rNodeList[v_i_k].previous.size(); i++) { rqarc = rNodeList[v_i_k].previous[i]; // when analysing d(p_k) the arcs starting in d(p_k) which belong to other loopless paths already determined should not be used. if (is_in_d[v_i_k] == true && determinedStartingArcs[v_i_k].count(rqarc->SourceID) != 0) // this arc has belong to another loopless path continue; if (rqarc->SourceID != v_ip1_k && isRemoved[rqarc->SourceID] == false) // arc (v_i^k, v_{i+1}^k) restore after { rqarc->cost = arcID2Cost.at(rqarc->ArcID); // equivalent to restore in the network if (Pi[v_i_k] > Pi[rqarc->SourceID] + rqarc->cost) Pi[v_i_k] = Pi[rqarc->SourceID] + rqarc->cost; // Calculate pi_{v_i^k} using forward star form } } isRemoved[v_i_k] = false; // Restore node v_i^k in the network // ========== If (pi_{v_i^k} is defined) Then ========== if (Pi[v_i_k] != INF) { // ===== Correct labels of v_i^k successors using backward star form ===== correctLabels(NodeList, Pi, v_i_k, isRemoved); correctLabelsReverseGraph(rNodeList, Pi, v_i_k, isRemoved); // ===== p <-- sub_{p_k}(s, v_i^k) square T_t(v_i^k) ===== rDijkstra(rNodeList, Pi, NodeNum, dst); recordDirectRouteForNode(nodePath_t_2_v_i_k, pi, dst, v_i_k, false); recordDirectRouteForArc(arcPath_t_2_v_i_k, arcRecord, pi, dst, v_i_k); for (i = path.nodePath.size()-1; ; i--) { if (path.nodePath[i] == v_i_k) break; path.nodePath.pop_back(); path.arcPath.pop_back(); } for (i = nodePath_t_2_v_i_k.size()-1; i > 0; i--) { path.nodePath.push_back(nodePath_t_2_v_i_k[i-1]); path.arcPath.push_back(arcPath_t_2_v_i_k[i-1]); } // ===== d(p) <-- v_i^k ===== is_in_d[v_i_k] = true; nodePath_t_2_v_i_k.clear(); arcPath_t_2_v_i_k.clear(); // X <-- X U {p} path.costSum = 0; for (i = 0; i < path.arcPath.size(); i++) path.costSum += arcID2Cost.at(path.arcPath[i]); path.deviation = v_i_k; path.deviationFromWhichPath = k - 1; Path *newPath = new Path(src); *newPath = path; X.insert(newPath); } // EndIf // ========== Restore (v_i^k, v_{i+1}^k) in the network ========== // forward graph, out-degree parc = NodeList[v_i_k].firstArc; while (parc != NULL) { if (parc->DestinationID == v_ip1_k) // now, arc (v_i^k, v_{i+1}^k) restore { parc->cost = arcID2Cost.at(parc->ArcID); // equivalent to restore in the network if (Pi[v_i_k] > Pi[v_ip1_k] + parc->cost) Pi[v_i_k] = Pi[v_ip1_k] + parc->cost; // Calculate pi_{v_i^k} using forward star form // Correct labels of v_i^k successors using backward star form correctLabels(NodeList, Pi, v_i_k, isRemoved); break; } parc = parc->nextArc; } // reverse graph, in-degree for (i = 0; i < rNodeList[v_i_k].previous.size(); i++) { rqarc = rNodeList[v_i_k].previous[i]; if (rqarc->SourceID == v_ip1_k) // now, arc (v_i^k, v_{i+1}^k) restore { rqarc->cost = arcID2Cost.at(rqarc->ArcID); // equivalent to restore in the network correctLabelsReverseGraph(rNodeList, Pi, v_i_k, isRemoved); break; } } v_ip1_k = v_i_k; } while (v_i_k != d_p_k[k-1]); // EndFor // ========== Restore nodes and arcs of path p_k from s to d(p_k) in the network ========== for (i = 0; i < path_k.nodePath.size(); i++) // sub_{p_k}(s, v_i^k), because d[k] = v_i^k { int restoring = path_k.nodePath[i]; if (isRemoved[restoring] == false) // if it has been restored, skip it to save time continue; // forward graph, out-degree parc = NodeList[restoring].firstArc; while (parc != NULL) { parc->cost = arcID2Cost.at(parc->ArcID); // equivalent to restore in the network parc = parc->nextArc; } // forward graph, in-degree for (j = 0; j < NodeList[restoring].previous.size(); j++) { qarc = NodeList[restoring].previous[j]; qarc->cost = arcID2Cost.at(qarc->ArcID); // equivalent to restore in the network } // reverse graph, out-degree rparc = rNodeList[restoring].firstArc; while (rparc != NULL) { rparc->cost = arcID2Cost.at(rparc->ArcID); // equivalent to restore in the network rparc = rparc->nextArc; } // reverse graph, in-degree for (j = 0; j < rNodeList[restoring].previous.size(); j++) { rqarc = rNodeList[restoring].previous[j]; rqarc->cost = arcID2Cost.at(rqarc->ArcID); // equivalent to restore in the network } isRemoved[restoring] = false; } // ========== Restore arcs (d(p_k), i), i in N, of p_1, ..., p_{k-1} in the network ========== if (k != 1) { for (i = 0; i < k-1; i++) { for (j = 0; j < erasedPath[i]->nodePath.size() && j < path_k.nodePath.size(); j++) { if (erasedPath[i]->nodePath[j] != path_k.nodePath[j]) { deviationFrom = path_k.nodePath[j-1]; break; } } if (j != erasedPath[i]->nodePath.size() && j != path_k.nodePath.size()) determinedStartingArcs[deviationFrom].erase(erasedPath[i]->nodePath[j]); } } // ========== Restore arcs starting in d(p_k) which belong to other loopless paths in the network ========== nodePathTemp.assign(path_k.nodePath.begin(), path_k.nodePath.end()); v_ip1_k = nodePathTemp.back(); nodePathTemp.pop_back(); do { v_i_k = nodePathTemp.back(); nodePathTemp.pop_back(); // forward graph, out-degree parc = NodeList[v_i_k].firstArc; while (parc != NULL) { if (parc->DestinationID != v_ip1_k && isRemoved[parc->DestinationID] == false) // arc (v_i^k, v_{i+1}^k) restore after parc->cost = arcID2Cost.at(parc->ArcID); // equivalent to restore in the network parc = parc->nextArc; } // reverse graph, in-degree for (i = 0; i < rNodeList[v_i_k].previous.size(); i++) { rqarc = rNodeList[v_i_k].previous[i]; if (rqarc->SourceID != v_ip1_k && isRemoved[rqarc->SourceID] == false) // arc (v_i^k, v_{i+1}^k) restore after rqarc->cost = arcID2Cost.at(rqarc->ArcID); // equivalent to restore in the network } v_ip1_k = v_i_k; } while (v_i_k != d_p_k[k-1]); // EndFor k++; } path = **X.begin(); for (set<Path*, bool (*)(const Path*, const Path*)>::iterator its = X.begin(); its != X.end(); ++its) delete *its; X.erase(X.begin(), X.end()); for (vector<Path*>::iterator itv = erasedPath.begin(); itv != erasedPath.end(); ++itv) delete *itv; erasedPath.clear(); delete []determinedStartingArcs; delete []d_p_k; delete []Pi; }
针对下图的测试结果如下:
至于正向图和方向图的构建在此给个参考,知道其中值的存储形式就很容易进行移植:
// 3. ========== construct graph ========== sort(AllUnion.begin(), AllUnion.end(), sourceCompare); Node *NodeList = new Node[NodeNum]; unsigned int NodeId = 0; unsigned int rank = 2; Arc *firstarc = new Arc, *parc, *qarc; firstarc->ArcID = AllUnion[0].second.first; firstarc->SourceID = AllUnion[0].first.first; firstarc->DestinationID = AllUnion[0].first.second; firstarc->cost = AllUnion[0].second.second; firstarc->nextArc = NULL; NodeId = AllUnion[0].first.first; // node index start with 0 arcIsAdded[NodeId] = true; NodeList[NodeId].ID = nodeIndex2ID.at(firstarc->SourceID); NodeList[NodeId].firstArc = firstarc; NodeList[firstarc->DestinationID].previous.push_back(firstarc); for (i = 1; i < ArcNum; i++) { if (AllUnion[i].first.first == AllUnion[i-1].first.first) { if (rank % 2 == 0) { parc = new Arc; parc->ArcID = AllUnion[i].second.first; parc->SourceID = AllUnion[i].first.first; parc->DestinationID = AllUnion[i].first.second; parc->cost = AllUnion[i].second.second; parc->nextArc = NULL; if (rank == 2) firstarc->nextArc = parc; else qarc->nextArc = parc; NodeList[parc->DestinationID].previous.push_back(parc); } else { qarc = new Arc; qarc->ArcID = AllUnion[i].second.first; qarc->SourceID = AllUnion[i].first.first; qarc->DestinationID = AllUnion[i].first.second; qarc->cost = AllUnion[i].second.second; qarc->nextArc = NULL; parc->nextArc = qarc; NodeList[qarc->DestinationID].previous.push_back(qarc); } rank++; } else { rank = 2; NodeId = AllUnion[i].first.first; arcIsAdded[NodeId] = true; firstarc = new Arc; firstarc->ArcID = AllUnion[i].second.first; firstarc->SourceID = NodeId; firstarc->DestinationID = AllUnion[i].first.second; firstarc->cost = AllUnion[i].second.second; firstarc->nextArc = NULL; NodeList[NodeId].ID = nodeIndex2ID.at(NodeId); NodeList[NodeId].firstArc = firstarc; NodeList[firstarc->DestinationID].previous.push_back(firstarc); } } for (i = 0; i < NodeNum; i++) { if (arcIsAdded[i] == false) { NodeList[i].ID = nodeIndex2ID.at(i); NodeList[i].firstArc = NULL; } } // 4. ========== construct reverse graph ========== sort(AllUnion.begin(), AllUnion.end(), terminalCompare); for (i = 0; i < NodeNum; i++) arcIsAdded[i] = false; NodeId = 0; rank = 2; rNode *rNodeList = new rNode[NodeNum]; rArc *rfirstarc = new rArc, *rparc, *rqarc; rfirstarc->ArcID = AllUnion[0].second.first; rfirstarc->SourceID = AllUnion[0].first.second; rfirstarc->DestinationID = AllUnion[0].first.first; rfirstarc->cost = AllUnion[0].second.second; rfirstarc->nextArc = NULL; NodeId = AllUnion[0].first.second; // NodeId ?= 0, maybe not arcIsAdded[NodeId] = true; rNodeList[NodeId].ID = nodeIndex2ID.at(rfirstarc->SourceID); rNodeList[NodeId].firstArc = rfirstarc; rNodeList[rfirstarc->DestinationID].previous.push_back(rfirstarc); for (i = 1; i < ArcNum; i++) { if (AllUnion[i].first.second == AllUnion[i-1].first.second) { if (rank % 2 == 0) { rparc = new rArc; rparc->ArcID = AllUnion[i].second.first; rparc->SourceID = AllUnion[i].first.second; rparc->DestinationID = AllUnion[i].first.first; rparc->cost = AllUnion[i].second.second; rparc->nextArc = NULL; if (rank == 2) rfirstarc->nextArc = rparc; else rqarc->nextArc = rparc; rNodeList[rparc->DestinationID].previous.push_back(rparc); } else { rqarc = new rArc; rqarc->ArcID = AllUnion[i].second.first; rqarc->SourceID = AllUnion[i].first.second; rqarc->DestinationID = AllUnion[i].first.first; rqarc->cost = AllUnion[i].second.second; rqarc->nextArc = NULL; rparc->nextArc = rqarc; rNodeList[rqarc->DestinationID].previous.push_back(rqarc); } rank++; } else { rank = 2; NodeId = AllUnion[i].first.second; arcIsAdded[NodeId] = true; rfirstarc = new rArc; rfirstarc->ArcID = AllUnion[i].second.first; rfirstarc->SourceID = NodeId; rfirstarc->DestinationID = AllUnion[i].first.first; rfirstarc->cost = AllUnion[i].second.second; rfirstarc->nextArc = NULL; rNodeList[NodeId].ID = nodeIndex2ID.at(NodeId); rNodeList[NodeId].firstArc = rfirstarc; rNodeList[rfirstarc->DestinationID].previous.push_back(rfirstarc); } } for (i = 0; i < NodeNum; i++) { if (arcIsAdded[i] == false) { rNodeList[i].ID = nodeIndex2ID.at(i); rNodeList[i].firstArc = NULL; } }