1
#include <Rcpp.h>
2
// #include <cstddef>
3
using namespace Rcpp;
4

5
// typedef int_least64_t linex;
6

7 1
double AssignTopDown(int v, 
8
                  NumericVector &lef, 
9
                  IntegerVector &visited,
10
                  std::vector<std::vector<int> > &ImSucc
11
                  ){
12 1
  visited[v]=1;
13 1
  double e=0;
14 1
  for(std::vector<int>::size_type i = 0; i!=ImSucc[v].size(); ++i){
15
  // for(int i=0; i<ImSucc[v].size(); ++i){
16 1
    int vPrime=ImSucc[v][i];
17 1
    if(vPrime==0){
18 1
      e+=1;
19 1
      lef[vPrime]=1;
20
    }
21
    else{
22 1
      if(visited[vPrime]==0){
23 1
        e=e+AssignTopDown(vPrime,lef,visited,ImSucc);
24
      }
25
      else{
26 1
        e=e+lef[vPrime];
27
      }
28
    }
29
  }
30 1
  lef[v]=e;
31 1
  return e;
32
}
33

34 1
void AssignBottomUp(int nElem,
35
                    NumericVector &lei, 
36
                    IntegerVector &visited,
37
                    std::vector<std::vector<int> > &ImSucc){
38 1
  std::vector<int> Q;
39

40 1
  lei(nElem)=1;
41

42 1
  for(std::vector<int>::size_type i = 0; i!=ImSucc[nElem].size();++i){
43
  // for(int i=0; i<ImSucc[nElem].size();++i){
44 1
    int vPrime=ImSucc[nElem][i];
45 1
    Q.push_back(vPrime);
46 1
    lei[vPrime]=1;
47
  }
48 1
  while(Q.size()!=0){
49 1
    Rcpp::checkUserInterrupt();
50 1
    int v=Q[0];
51 1
    Q.erase(Q.begin());
52

53 1
    for(std::vector<int>::size_type j = 0; j!=ImSucc[v].size();++j){
54
    // for(int j=0;j<ImSucc[v].size();++j){
55 1
      int vPrime=ImSucc[v][j];
56 1
      lei[vPrime]+=lei[v];
57 1
      if(visited[vPrime]==0){
58 1
        Q.push_back(vPrime);
59 1
        visited[vPrime]=1;
60
      }
61
    }
62
  }
63
}
64

65 1
void ComputeRankProb(int v,int h, NumericMatrix &rp,
66
                     std::vector<std::vector<int> > &ImSucc,
67
                     std::vector<std::vector<int> > &ideals,
68
                     IntegerVector &visited,
69
                     NumericVector &lei,
70
                     NumericVector &lef,
71
                     double &e){
72 1
  visited[v]=1;
73
  // std::vector<int>::size_type i = 0; i!=child[0].size(); ++i
74 1
  for(std::vector<int>::size_type j = 0;j!=ImSucc[v].size();++j){
75
    // for(int j=0;j<ImSucc[v].size();++j){
76 1
    int vPrime=ImSucc[v][j];
77
    int x;
78 1
    std::set_difference(ideals[vPrime].begin(), ideals[vPrime].end(),
79 1
                        ideals[v].begin(), ideals[v].end(), &x);
80 1
    rp(x,h)=rp(x,h)+double(lei[v])*double(lef[vPrime])/double(e);
81 1
    if((vPrime!=0) & (visited[vPrime]==0)){
82 1
      ComputeRankProb(vPrime,h+1,rp,ImSucc,ideals,visited,lei,lef,e);
83
    }
84
  }
85
  
86
}
87

88 1
void ComputeMutualRankProb(int v,int h, int &nElem,
89
                           NumericMatrix &mrp,
90
                           std::vector<std::vector<int> > &ImSucc,
91
                           std::vector<std::vector<int> > &ideals,
92
                           IntegerVector &visited,
93
                           IntegerVector &visitedElem,
94
                           NumericVector &lei,
95
                           NumericVector &lef,
96
                           double &e){
97 1
  visited[v]=1;
98 1
  for(std::vector<int>::size_type j = 0;j!=ImSucc[v].size();++j){
99
  // for(int j=0;j<ImSucc[v].size();++j){
100 1
    int vPrime=ImSucc[v][j];
101 1
    for(int y=0;y<nElem; ++y){
102 1
      if(visitedElem[y]==1){
103
        int x;
104 1
        std::set_difference(ideals[vPrime].begin(), ideals[vPrime].end(),
105 1
                            ideals[v].begin(), ideals[v].end(), &x);
106 1
        mrp(x,y)=mrp(x,y)+double(lei[v])*double(lef[vPrime])/double(e);
107
      }
108 1
      if((vPrime!=0) & (visited[vPrime]==0)){
109
        int x;
110 1
        std::set_difference(ideals[vPrime].begin(), ideals[vPrime].end(),
111 1
                            ideals[v].begin(), ideals[v].end(), &x);
112 1
        visitedElem[x]=1;
113 1
        ComputeMutualRankProb(vPrime,h+1,nElem,mrp,ImSucc,ideals,visited,visitedElem,lei,lef,e);
114 1
        visitedElem[x]=0;
115
      }
116
        
117
    }
118
  }
119
  
120
  
121
}
122

123
// [[Rcpp::export]]
124 1
Rcpp::List rankprobs(std::vector<std::vector<int> > ImPred,
125
                     std::vector<std::vector<int> > ideals,
126
                     int nElem,
127
                     int nIdeals){
128
  double e;
129 1
  NumericVector lei(nIdeals);
130 1
  NumericVector lef(nIdeals);
131 1
  IntegerVector visited(nIdeals);
132 1
  IntegerVector visitedElem(nElem);
133
  
134
  /* Turn ImPred to ImSucc*/
135 1
  std::vector<std::vector<int> > ImSucc(nIdeals);
136 1
  for(int i=0; i<nIdeals;++i){
137 1
    for(std::vector<int>::size_type j=0;j!=ImPred[i].size();++j){
138
    // for(int j=0;j<ImPred[i].size();++j){
139 1
      int idx=ImPred[i][j];
140 1
      ImSucc[idx].push_back(i);
141
    }
142
  }
143
  /* Sort increasingly*/
144 1
  for(int i=0;i<nIdeals;++i){
145 1
    std::sort(ImSucc[i].begin(), ImSucc[i].end());
146
  }
147
  /*calculate number of path*/
148 1
  e=AssignTopDown(nElem, lef,visited,ImSucc);
149 1
  std::fill(visited.begin(), visited.end(), 0);
150 1
  AssignBottomUp(nElem,lei,visited,ImSucc);
151

152
  /*rank probabilities*/
153 1
  std::fill(visited.begin(), visited.end(), 0);
154 1
  NumericMatrix rp(nElem,nElem);
155 1
  ComputeRankProb(nElem,0,rp,ImSucc,ideals,visited,lei,lef,e);
156
  
157
  /*mutual rank probabilities*/
158 1
  std::fill(visited.begin(), visited.end(), 0);
159 1
  NumericMatrix mrp(nElem,nElem);
160 1
  ComputeMutualRankProb(nElem,1,nElem,mrp,ImSucc,ideals,visited,visitedElem,lei,lef,e);
161 1
  return Rcpp::List::create(Rcpp::Named("linext") = e, 
162 1
                            Rcpp::Named("rp")=rp,
163 1
                            Rcpp::Named("mrp")=mrp);
164
}
165

166

Read our documentation on viewing source code .

Loading