城市间最短路径问题——R和Rcpp实现
这里的最短路径问题也叫做相识问题，具体问题来自
刚好在学Rcpp，于是写了R版本和Rcpp版本的函数来解决这个问题，顺便比较一下二者的运行速度。
Rcpp包为在R中使用C++提供了很便捷的方式，这让R代码的运行速度有了很大的提升。使用Rcpp对C++有个基本了解就行，不需要对C++很熟练。学习Rcpp包的资源有很多，这里推荐一本对像我这种C++小白比较友好的电子书：https://teuder.github.io/rcpp4everyone_en/。Rcpp包的安装和配置等细节这里就不提了，感兴趣的读者可以自行搜索，网上资源很多。
-
R版本，R_shortest()函数定义如下：
R_shortest<-function(M){
n<-max(M)
A<-matrix(Inf,n,n)
for(i in 1:nrow(M)){A[M[i,1],M[i,2]]<-A[M[i,2],M[i,1]]<-1}
diag(A)<-0
while(TRUE){
B<-A
for(i in 1:n){
for(j in 1:n){
for(k in 1:n){
if(A[i,j]>A[i,k]+A[k,j]){
A[i,j]<-A[i,k]+A[k,j]
}
}
}
}
if(identical(B,A)){break}else{B<-A}
}
return(A)
}
-
Rcpp版本，Rcpp_shortest()函数定义如下：
首先是创建一个名为shortest的cpp文件，文件内容如下：
#include <Rcpp.h> using namespace Rcpp; // [[Rcpp::export]] NumericMatrix Rcpp_shortest(NumericMatrix M){ M=M-1; //注意，c++下标从0开始，而R下标从1开始，因此要减去1 int n=max(M(_,1))+1; //由于减了1，因此新矩阵的行数应该是现在M的最大值加1 NumericMatrix A(n,n); A.fill(R_PosInf); A.fill_diag(0); for(int i=0;i<M.nrow();i++){ A(M(i,0),M(i,1))=1; A(M(i,1),M(i,0))=1; } while(true){ NumericMatrix B=clone(A); for(int i=0;i<n;i++){ for(int j=0;j<n;j++){ for(int k=0;k<n;k++){ if(A(i,j)>A(i,k)+A(k,j)){ A(i,j)=A(i,k)+A(k,j); } } } } if(sum(A!=B)==0){break;}else{B=clone(A);} } return A; }
然后是对这个文件在R中调用，调用之后会得到一个名为Rcpp_shortest()的函数：
library(Rcpp)
sourceCpp("shortest.cpp")
-
速度比较
M<-matrix(c(1,2,2,3,3,4,4,5,5,8,7,9,
9,21,10,11,11,12,13,21,21,100),
ncol=2,byrow=TRUE) #生成包含100个城市的矩阵M；
identical(R_shortest(M),Rcpp_shortest(M)) #考虑返回结果是否一致
library(rbenchmark) #速度比较包
benchmark(R_shortest(M),Rcpp_shortest(M),
columns=c("test","replications","elapsed","relative")) #速度比较
> identical(R_shortest(M),Rcpp_shortest(M))
[1] TRUE
两个函数运行结果一致，速度比较结果如下：
test replications elapsed relative
1 R_shortest(M) 100 110.70 307.5
2 Rcpp_shortest(M) 100 0.36 1.0
在本人电脑上，R版本用时是Rcpp版本的307倍，Rcpp版本的速度远快于R版本，加速效果不错！所以说，懂得一点C++的话，利用Rcpp来编写循环和四则运算的话，还是可以达到对R加速的目的。此外，由于Rcpp的存在，Rcpp版本的函数和R版本的函数其实很像，因此将R版本的函数改写为Rcpp版本的函数难度也不会很大。
- 附带一个Python版本吧，权当练练手，不过这就无关主题了。
import pandas as pd import numpy as np def Py_shortest(M): M=M-1 n=M[:,1].max()+1 A=np.full(shape=(n,n),fill_value=float("inf")) np.fill_diagonal(a=A,val=0) for i in range(0,M.shape[0]): A[M[i,0],M[i,1]]=1 A[M[i,1],M[i,0]]=1 while(True): B=A.copy() for i in range(0,n): for j in range(0,n): for k in range(0,n): if A[i,j]>A[i,k]+A[k,j]: A[i,j]=A[i,k]+A[k,j] if((A==B).all()): break else: B=A.copy() return A M=np.array([[1,2],[2,3],[3,4],[4,5],[5,8],[7,9], [9,21],[10,11],[11,12],[13,21],[21,100]]) import time
time_start=time.time() Py_shortest(M) time_end=time.time() print(time_end-time_start)
本人电脑运行一次Py_shortest(M)的时间在3-4秒内，Python下标也从0开始，写程序时要注意和R、Rcpp的区别和联系。
