1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
#include <Rcpp.h>
using namespace Rcpp;

//' fast Euclidean distance matrix
//'
//' @param x matrix with sample rows for which the distance matrix is computed (to use with vectors, use \code{as.matrix(x)})
//' @examples
//' #require(microbenchmark)
//' #x = rnorm(100)
//' #microbenchmark(fastdist(as.matrix(x)),as.matrix(dist(x)))
//' @export
// [[Rcpp::export]]
NumericMatrix fastdist (const NumericMatrix & x){
  unsigned int outrows = x.nrow(), i = 0, j = 0;
  double d;
  Rcpp::NumericMatrix out(outrows,outrows);

  for (i = 0; i < outrows - 1; i++){
    Rcpp::NumericVector v1 = x.row(i);
    for (j = i + 1; j < outrows ; j ++){
      d = sqrt(sum(pow(v1-x.row(j), 2.0)));
      out(j,i)=d;
      out(i,j)=d;
    }
  }

  return out;
}

//' double center a symmetric matrix
//'
//' @param x symmetric matrix
//' @param normalize boolean. If \code{TRUE} the matrix will be normalized to mean 1.
//' @keywords internal
// [[Rcpp::export]]
NumericMatrix doubleCenterSymMat(const NumericMatrix & x, bool & normalize) {
  int i, j;
  int N = x.nrow();
  NumericVector colmeans(N);
  NumericMatrix out(N, N);
  double fullmean;
  double tmp;

  for (i=0; i<N; i++) {
    colmeans(i) = sum(x(i,_))/(double)(N);
  }
  fullmean = sum(colmeans)/N;

  if ( (fullmean == 0) | !normalize) {
    if (fullmean == 0) warning("It seems that one variable is constant. Constants are always independent. \n");
    for (i=0; i<N; i++)
      for (j=i; j<N; j++) {
        tmp = - x(i, j) + colmeans(i) + colmeans(j) - fullmean;
        out(j, i) = tmp;
        out(i, j) = tmp;
      }
  } else {
    for (i=0; i<N; i++)
      for (j=i; j<N; j++) {
        tmp = (- x(i, j) + colmeans(i) + colmeans(j) - fullmean)/fullmean;
        out(j, i) = tmp;
        out(i, j) = tmp;
      }
  }

  return out;
}


//' fast centered Euclidean distance matrix
//'
//' @param x matrix with sample rows for which the distance matrix is computed (to use with vectors, use \code{as.matrix(x)})
//' @param normalize boolean. If \code{TRUE} the matrix will be normalized to mean 1.
//' @export
// [[Rcpp::export]]
NumericMatrix fastEuclideanCdm (const NumericMatrix & x, bool & normalize){
  unsigned const int N = x.nrow();
  unsigned int i = 0, j = 0;
  NumericMatrix out(N,N);
  NumericVector colmeans(N);
  double tmp, m;


  for (i = 0; i < N - 1; i++){ // row
    NumericVector v1 = x.row(i);
    for (j = i + 1; j < N ; j ++){ // column
      tmp = sqrt(sum(pow(v1-x.row(j), 2.0)));
      out(i,j) = tmp;
      colmeans(i) += tmp;
      colmeans(j) += tmp;
    }
  }

  colmeans = colmeans/N;
  m = sum(colmeans)/(double) N;

  if ( (m == 0) | !normalize) {
    if (m == 0) warning("It seems that one variable is constant. Constants are always independent. \n");
    for (i = 0; i < N-1; i++){ // row
      for (j = i+1 ; j < N ; j++){ // column
        tmp = -out(i,j) + colmeans(i) + colmeans(j) - m;
        out(i,j) = tmp;
        out(j,i) = tmp;

      }
    }


    for (i = 0; i < N ; i++){ // diag
      out(i,i) =  2*colmeans(i) - m;
    }
  } else {
    for (i = 0; i < N-1; i++){ // row
      for (j = i+1 ; j < N ; j++){ // column
        tmp = (-out(i,j) + colmeans(i) + colmeans(j) - m)/m;
        out(i,j) = tmp;
        out(j,i) = tmp;
      }
    }


    for (i = 0; i < N ; i++){ // diag
      out(i,i) =  (2*colmeans(i) - m)/m;
    }
  }
  /* */
   return out;
}


//' for the fast detection of the full dependence structure
//'
//' Returns the row indicies of matrix A which match with B
//'
//' @param A matrix
//' @param B matrix whose rows are subset of A
//'
//' @examples
//' # A = t(utils::combn(10,3))
//' # B = A[sort(sample.int(nrow(A),10)),]
//' # match_rows(A,B)
//'
//' @keywords internal
// [[Rcpp::export]]

NumericVector match_rows(NumericMatrix & A,NumericMatrix &B){
  int i = 0, k;
  NumericVector res (B.nrow());

  for (k = 0; k < B.nrow(); k++) {
    while( is_true(any(A.row(i) != B.row(k)))) {
      i++;
    }
    res(k) = i;
  }
  return res+1;
}



/*** R
*/