0

我正在尝试使用 Rcpp 和 RcppParallel 构建一个实现 Dijkstra 算法的 R 包。你可以在这里看到我的作品。现在我想添加一个新功能并且出现了一个奇怪的行为。当我通过 sourceCpp 函数编译这个函数并尝试它时，它运行良好，但是当我用这个函数重建包时，当我调用该函数时，R 几乎立即崩溃。当我通过 Rstudio（清理和重建）构建包时，我没有错误。
该算法是双向的 Dijkstra，因此它使用两个图（正向和反向）和来自 STL 的两个优先级队列。图只是成对向量的向量。

这是代码（对不起，我不知道如何简化它）：

// [[Rcpp::depends(RcppParallel)]]
// [[Rcpp::plugins(cpp11)]]

#include <queue>
#include <vector>
#include <limits>
#include <functional>
#include <RcppParallel.h>
#include <Rcpp.h>
using namespace RcppParallel;

struct Pardijkstra : public Worker
{
  //input
  const std::vector<std::vector<std::pair<int, double> > > m_graph;
  const std::vector<std::vector<std::pair<int, double> > > m_graphr;
  RVector<int> m_dep;
  RVector<int> m_arr;
  const int m_nbnodes;

  //output
  RcppParallel::RVector<double> m_result;

  //constructor
  Pardijkstra(const std::vector<std::vector<std::pair<int, double> > >  &graph,
              const std::vector<std::vector<std::pair<int, double> > >  &graphr,
              Rcpp::IntegerVector & dep,
              Rcpp::IntegerVector & arr,
              const int nbnodes,
              Rcpp::NumericVector result) : m_graph(graph),m_graphr(graphr),m_dep(dep), m_arr(arr),m_nbnodes(nbnodes),m_result(result)
  {

  }

  //overload () operator
  void operator()(std::size_t begin, std::size_t end){
    struct comp{

      //Custom comparator included in priority queues

      bool operator()(const std::pair<int, double> &a, const std::pair<int, double> &b){
        return a.second > b.second;
      }
    };

    //

    for (std::size_t k=begin; k!=end;k++){

      //Here is the algorithm (bidirectional search)

      int StartNode=m_dep[k];
      int EndNode=m_arr[k];

      std::vector<double> Distances(m_nbnodes, std::numeric_limits<double>::max()); 
      std::vector<double> Distances2(m_nbnodes, std::numeric_limits<double>::max()); 

      //Distances are stored here
      Distances[StartNode] = 0.0;  
      Distances2[EndNode] = 0.0;

      std::vector <int> Visited(m_nbnodes,0);

      std::priority_queue<std::pair<int, double>, std::vector<std::pair<int, double> >, comp > Q;
      std::priority_queue<std::pair<int, double>, std::vector<std::pair<int, double> >, comp > Qr;
      Q.push(std::make_pair(StartNode, 0.0)); 
      Qr.push(std::make_pair(EndNode, 0.0));
      Visited[StartNode]+=1;
      Visited[EndNode]+=1;

      double mu=std::numeric_limits<double>::max();

      while (!Q.empty() && !Qr.empty()) {  

        //Stopping criterion

        if (Q.top().second+Qr.top().second >= mu){
          break;
        }  

        //Forward search
        if (!Q.empty()){
          int v=Q.top().first;
          int w=Q.top().second;
          Q.pop();

          if (w <= Distances[v]) {
            for (int i=0; i< m_graph[v].size(); i++){
              int v2 = m_graph[v][i].first;                                                     
              double w2 = m_graph[v][i].second;

              if (Distances[v] + w2 < Distances[v2]) {                               
                Distances[v2] = Distances[v] + w2;                                   

                Q.push(std::make_pair(v2, Distances[v2]));
                Visited[v2]+=1;

              }
            }
          }
          if ((Visited[v]>1)  && (Distances[v]+Distances2[v]) < mu){

            mu=Distances[v]+Distances2[v];

          }
        }

        //Backward search
        if (!Qr.empty()){
          int vv=Qr.top().first;
          int ww=Qr.top().second;
          Qr.pop();

          Visited[vv]+=1;

          if (ww <= Distances2[vv]) {
            for (int i=0; i< m_graphr[vv].size(); i++){
              int vv2 = m_graphr[vv][i].first;                                                     
              double ww2 = m_graphr[vv][i].second;

              if (Distances2[vv] + ww2 < Distances2[vv2]) {                               
                Distances2[vv2] = Distances2[vv] + ww2;                                  

                Qr.push(std::make_pair(vv2, Distances2[vv2]));
                Visited[vv2]+=1;
              }
            }
          }
          if ((Visited[vv]> 1) && (Distances[vv]+Distances2[vv]) < mu){
            mu=Distances[vv]+Distances2[vv];
          }
        }
      }

      if (mu >= std::numeric_limits<double>::max()){
        m_result[k] = Rcpp::NumericVector::get_na();
      }
      else {
        m_result[k]=mu;
      }

    }

  }

};


//Fonction exported in R

// [[Rcpp::export]]
Rcpp::NumericVector Bidir_par(Rcpp::IntegerVector dep, Rcpp::IntegerVector arr,Rcpp::IntegerVector gfrom,Rcpp::IntegerVector gto,Rcpp::NumericVector gw,int NbNodes){


  //dep : origin nodes
  //arr: destination nodes
  //gfrom: first nodes of graphs edges
  //gto: last nodes of graphs edges
  //gw: weights of graphs edges
  //NbNodes: number of nodes in the graph

  Rcpp::NumericVector result(dep.size());


  //Construction of the two graphs
  int NbEdges=gfrom.size();

  std::vector<std::vector<std::pair<int, double> > > G(NbNodes);   
  std::vector<std::vector<std::pair<int, double> > > Gr(NbNodes);
  for (int i = 0; i != NbEdges; ++i) {

    G[gfrom[i]].push_back(std::make_pair(gto[i], gw[i]));
    Gr[gto[i]].push_back(std::make_pair(gfrom[i], gw[i]));

  }


  Pardijkstra Dijfunc(G,Gr,dep,arr,NbNodes,result);
  parallelFor(0,dep.length(),Dijfunc);

  return Rcpp::wrap(result);



}

我知道输入应该是 RVector 或 RMatrix，但不幸的是，有向图不能以这种形式简化而不会降低效率。std::vectors 在内存中是连续的，对吗？

我注意到，如果我删除一个 if 语句（例如所有向后搜索），它就会起作用。

如果您想查看包中实现的其他功能，可以查看我的 github 存储库。所有其他功能都运行良好。
这是 Makevars 文件： 

CXX_STD = CXX11
PKG_LIBS += $(shell ${R_HOME}/bin/Rscript -e "RcppParallel::RcppParallelLibs()")

Makevars.win 文件： 

CXX_STD = CXX11
PKG_CXXFLAGS += -DRCPP_PARALLEL_USE_TBB=1
PKG_LIBS += $(shell "${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe" -e "RcppParallel::RcppParallelLibs()")

和说明文件：

Imports: Rcpp (>= 1.0.1), RcppParallel (>= 4.4.2)
LinkingTo: Rcpp, RcppParallel
SystemRequirements: GNU make

那么我做错了什么？
为什么该函数与 sourceCpp 函数一起使用而不是在包中？

编辑：根据@thc，我的实现不是线程安全的，所以有没有办法以安全的方式重写这个算法？知道优先级队列在这里是强制性的。

任何帮助表示赞赏！

4

0 回答 0