城市间最短路径问题—

　　这里的最短路径问题也叫做相识问题，具体问题来自

https://www.math.pku.edu.cn/teachers/lidf/docs/Rbook/html/_Rbook/examples.html#examples-scicomp-citylink

　　刚好在学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的区别和联系。

巴特西

城市间最短路径问题——R和Rcpp实现

R版本，R_shortest()函数定义如下：

Rcpp版本，Rcpp_shortest()函数定义如下：

速度比较

最新文章

热门文章