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
|
/**
* @title Parallel Distance Matrix Calculation with RcppParallel
* @author JJ Allaire and Jim Bullard
* @license GPL (>= 2)
*/
#include <Rcpp.h>
using namespace Rcpp;
#include <cmath>
#include <algorithm>
// generic function for kl_divergence
template <typename InputIterator1, typename InputIterator2>
inline double kl_divergence(InputIterator1 begin1, InputIterator1 end1,
InputIterator2 begin2) {
// value to return
double rval = 0;
// set iterators to beginning of ranges
InputIterator1 it1 = begin1;
InputIterator2 it2 = begin2;
// for each input item
while (it1 != end1) {
// take the value and increment the iterator
double d1 = *it1++;
double d2 = *it2++;
// accumulate if appropirate
if (d1 > 0 && d2 > 0)
rval += std::log(d1 / d2) * d1;
}
return rval;
}
// helper function for taking the average of two numbers
inline double average(double val1, double val2) {
return (val1 + val2) / 2;
}
// [[Rcpp::export]]
NumericMatrix rcpp_js_distance(NumericMatrix mat) {
// allocate the matrix we will return
NumericMatrix rmat(mat.nrow(), mat.nrow());
for (int i = 0; i < rmat.nrow(); i++) {
for (int j = 0; j < i; j++) {
// rows we will operate on
NumericMatrix::Row row1 = mat.row(i);
NumericMatrix::Row row2 = mat.row(j);
// compute the average using std::tranform from the STL
std::vector<double> avg(row1.size());
std::transform(row1.begin(), row1.end(), // input range 1
row2.begin(), // input range 2
avg.begin(), // output range
average); // function to apply
// calculate divergences
double d1 = kl_divergence(row1.begin(), row1.end(), avg.begin());
double d2 = kl_divergence(row2.begin(), row2.end(), avg.begin());
// write to output matrix
rmat(i,j) = std::sqrt(.5 * (d1 + d2));
}
}
return rmat;
}
// [[Rcpp::depends(RcppParallel)]]
#include <RcppParallel.h>
using namespace RcppParallel;
struct JsDistance : public Worker {
// input matrix to read from
const RMatrix<double> mat;
// output matrix to write to
RMatrix<double> rmat;
// initialize from Rcpp input and output matrixes (the RMatrix class
// can be automatically converted to from the Rcpp matrix type)
JsDistance(const NumericMatrix mat, NumericMatrix rmat)
: mat(mat), rmat(rmat) {}
// function call operator that work for the specified range (begin/end)
void operator()(std::size_t begin, std::size_t end) {
for (std::size_t i = begin; i < end; i++) {
for (std::size_t j = 0; j < i; j++) {
// rows we will operate on
RMatrix<double>::Row row1 = mat.row(i);
RMatrix<double>::Row row2 = mat.row(j);
// compute the average using std::tranform from the STL
std::vector<double> avg(row1.length());
std::transform(row1.begin(), row1.end(), // input range 1
row2.begin(), // input range 2
avg.begin(), // output range
average); // function to apply
// calculate divergences
double d1 = kl_divergence(row1.begin(), row1.end(), avg.begin());
double d2 = kl_divergence(row2.begin(), row2.end(), avg.begin());
// write to output matrix
rmat(i,j) = sqrt(.5 * (d1 + d2));
}
}
}
};
// [[Rcpp::export]]
NumericMatrix rcpp_parallel_js_distance(NumericMatrix mat) {
// allocate the matrix we will return
NumericMatrix rmat(mat.nrow(), mat.nrow());
// create the worker
JsDistance jsDistance(mat, rmat);
// call it with parallelFor
parallelFor(0, mat.nrow(), jsDistance);
return rmat;
}
|
