mcmc_rank.cpp at d47ca11 ⋅ schochastics/netrankr

#include <Rcpp.h>

using namespace Rcpp;

// [[Rcpp::export]]
List mcmc_rank(IntegerMatrix P, 
               IntegerVector init_rank,
               int rp) {
  int t = 1;
  int n = P.rows();
  NumericMatrix rrp(n,n);
  NumericVector expected(n);
  IntegerVector elements=seq(0,n-1);

  //rrp
  for(int i=0; i<(n-1); ++i){
    int j=init_rank[i];
    IntegerVector tmp=init_rank[elements>i];
    for(int k=0;k<tmp.length();++k){
      rrp(j,tmp[k])=1;
    }
  }
  //expected
  for(int i=0;i<n; ++i){
      expected[init_rank[i]] =i;
  }
  
  //MC
  while(t<=rp){
    int p = floor(R::runif(0,1)*(n-1));//rand() % (n-1);
    int c = round(R::runif(0,1));//rand() % 2;
    int a = init_rank[p];
    int b = init_rank[p+1];
    t+=1;
    if((c==1) & (P(a,b)!=1)){
      init_rank[p]=b;
      init_rank[p+1]=a;
      //expected update
      for(int i=0;i<n; ++i){
        expected[init_rank[i]]=double(expected[init_rank[i]]*(t-1)+i)/double(t);
      }
      //rrp update
      for(int i=0; i<(n-1); ++i){
        int j=init_rank[i];
        IntegerVector tmp=init_rank[elements>i];
        for(int k=0;k<tmp.length();++k){
          rrp(j,tmp[k])=double(rrp(j,tmp[k])*(t-1)+1)/double(t);
          rrp(tmp[k],j)=1-rrp(j,tmp[k]);
        }
      }
    }
  }
  
  return Rcpp::List::create(Rcpp::Named("expected") = expected,
                            Rcpp::Named("rrp")=rrp);
}


Read our documentation on viewing source code .

1		#include <Rcpp.h>
2
3		using namespace Rcpp;
4
5		// [[Rcpp::export]]
6	1	List mcmc_rank(IntegerMatrix P,
7		IntegerVector init_rank,
8		int rp) {
9	1	int t = 1;
10	1	int n = P.rows();
11	1	NumericMatrix rrp(n,n);
12	1	NumericVector expected(n);
13	1	IntegerVector elements=seq(0,n-1);
14
15		//rrp
16	1	for(int i=0; i<(n-1); ++i){
17	1	int j=init_rank[i];
18	1	IntegerVector tmp=init_rank[elements>i];
19	1	for(int k=0;k<tmp.length();++k){
20	1	rrp(j,tmp[k])=1;
21		}
22		}
23		//expected
24	1	for(int i=0;i<n; ++i){
25	1	expected[init_rank[i]] =i;
26		}
27
28		//MC
29	1	while(t<=rp){
30	1	int p = floor(R::runif(0,1)*(n-1));//rand() % (n-1);
31	1	int c = round(R::runif(0,1));//rand() % 2;
32	1	int a = init_rank[p];
33	1	int b = init_rank[p+1];
34	1	t+=1;
35	1	if((c==1) & (P(a,b)!=1)){
36	0	init_rank[p]=b;
37	0	init_rank[p+1]=a;
38		//expected update
39	0	for(int i=0;i<n; ++i){
40	0	expected[init_rank[i]]=double(expected[init_rank[i]]*(t-1)+i)/double(t);
41		}
42		//rrp update
43	0	for(int i=0; i<(n-1); ++i){
44	0	int j=init_rank[i];
45	0	IntegerVector tmp=init_rank[elements>i];
46	0	for(int k=0;k<tmp.length();++k){
47	0	rrp(j,tmp[k])=double(rrp(j,tmp[k])*(t-1)+1)/double(t);
48	0	rrp(tmp[k],j)=1-rrp(j,tmp[k]);
49		}
50		}
51		}
52		}
53
54	1	return Rcpp::List::create(Rcpp::Named("expected") = expected,
55	1	Rcpp::Named("rrp")=rrp);
56		}
57

Navigation	Overlay
`t` Navigate files	`h` Toggle hits
`y` Change url to tip of branch	`m` Toggle misses
`b / v` Jump to prev/next hit line	`p` Toggle partial
`z / x` Jump to prev/next missed or partial line	`1..9` Toggle flags
`shift + o` Open current page in GitHub	`a` Toggle all on
`/ or ?` Show keyboard shortcuts dialog	`c` Toggle context lines or commits