Сессия R прервана с использованием RcppParallel - PullRequest
Новая
       308

Сессия R прервана с использованием RcppParallel

0 голосов
/ 13 июня 2019

Я пытаюсь создать пакет R, реализующий алгоритм Дейкстры с использованием Rcpp и RcppParallel. Вы можете увидеть мою работу здесь . Теперь я хочу добавить новую функцию и появляется странное поведение. Когда я компилирую эту функцию с помощью функции 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 :: векторы непрерывны в памяти, верно?

Я заметил, что если я удаляю оператор 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, а не в пакете?

EDIT : Согласно @thc, моя реализация не является поточно-ориентированной, так есть ли способ переписать этот алгоритм безопасным способом? Зная, что приоритетные очереди здесь обязательны.

Любая помощь приветствуется!

...