require(methods)

require(RcppParallel)

require(testthat)

RCPP <- Sys.getenv( "RCPP" )

if( RCPP == "Rcpp" ){

message( "testing against Rcpp" )

require(Rcpp)

} else if( RCPP == "Rcpp11" ){

message( "testing against Rcpp11" )

require(attributes)

} else {

stop( "Rcpp implementation not setup, please set the $RCPP environment variable" )

}

test_dir("testthat")

#include <Rcpp.h>

using namespace Rcpp;

#include <cmath>

#include <algorithm>

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;

}

inline double average(double val1, double val2) {

return (val1 + val2) / 2;

}

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;

}

#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));

}

}

}

};

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;

}

#include <Rcpp.h>

using namespace Rcpp;

#include <algorithm>

double innerProduct(NumericVector x, NumericVector y) {

return std::inner_product(x.begin(), x.end(), y.begin(), 0.0);

}

#include <RcppParallel.h>

using namespace RcppParallel;

struct InnerProduct : public Worker

{

// source vectors

const RVector<double> x;

const RVector<double> y;

// product that I have accumulated

double product;

// constructors

InnerProduct(const NumericVector x, const NumericVector y)

: x(x), y(y), product(0) {}

InnerProduct(const InnerProduct& innerProduct, Split)

: x(innerProduct.x), y(innerProduct.y), product(0) {}

// process just the elements of the range I have been asked to

void operator()(std::size_t begin, std::size_t end) {

product += std::inner_product(x.begin() + begin,

x.begin() + end,

y.begin() + begin,

0.0);

}

// join my value with that of another InnerProduct

void join(const InnerProduct& rhs) {

product += rhs.product;

}

};

double parallelInnerProduct(NumericVector x, NumericVector y) {

// declare the InnerProduct instance that takes a pointer to the vector data

InnerProduct innerProduct(x, y);

// call paralleReduce to start the work

parallelReduce(0, x.length(), innerProduct);

// return the computed product

return innerProduct.product;

}

#include <Rcpp.h>

#include <RcppParallel.h>

using namespace RcppParallel;

using namespace Rcpp;

struct Sum : public Worker

{

// source vector

const RVector<double> input;

// accumulated value

double value;

// constructors

Sum(const NumericVector input) : input(input), value(0) {}

Sum(const Sum& sum, Split) : input(sum.input), value(0) {}

// accumulate just the element of the range I have been asked to

void operator()(std::size_t begin, std::size_t end) {

value += std::accumulate(input.begin() + begin, input.begin() + end, 0.0);

}

// join my value with that of another Sum

void join(const Sum& rhs) {

value += rhs.value;

}

};

double parallelVectorSum(NumericVector x) {

// declare the SumBody instance

Sum sum(x);

// call parallel_reduce to start the work

parallelReduce(0, x.length(), sum);

// return the computed sum

return sum.value;

}

double vectorSum(NumericVector x) {

return std::accumulate(x.begin(), x.end(), 0.0);

}

#include <Rcpp.h>

using namespace Rcpp;

#include <cmath>

#include <algorithm>

NumericMatrix matrixSqrt(NumericMatrix orig) {

// allocate the matrix we will return

NumericMatrix mat(orig.nrow(), orig.ncol());

// transform it

std::transform(orig.begin(), orig.end(), mat.begin(), ::sqrt);

// return the new matrix

return mat;

}

#include <RcppParallel.h>

using namespace RcppParallel;

struct SquareRoot : public Worker

{

// source matrix

const RMatrix<double> input;

// destination matrix

RMatrix<double> output;

// initialize with source and destination

SquareRoot(const NumericMatrix input, NumericMatrix output)

: input(input), output(output) {}

// take the square root of the range of elements requested

void operator()(std::size_t begin, std::size_t end) {

std::transform(input.begin() + begin,

input.begin() + end,

output.begin() + begin,

::sqrt);

}

};

NumericMatrix parallelMatrixSqrt(NumericMatrix x) {

// allocate the output matrix

NumericMatrix output(x.nrow(), x.ncol());

// SquareRoot functor (pass input and output matrixes)

SquareRoot squareRoot(x, output);

// call parallelFor to do the work

parallelFor(0, x.nrow() * x.ncol(), squareRoot);

// return the output matrix

return output;

}

context( "distance" )

test_that( "distance works with Rcpp", {

sourceCpp( "cpp/distance.cpp" )

n = 1000

m = matrix(runif(n*10), ncol = 10)

m = m/rowSums(m)

expect_equal(

rcpp_js_distance(m),

rcpp_parallel_js_distance(m)

)

})

context( "inner product" )

test_that( "parallelInnerProduct works with Rcpp", {

sourceCpp( "cpp/innerproduct.cpp" )

x <- runif(1000000)

y <- runif(1000000)

expect_equal(innerProduct(x, y), parallelInnerProduct(x, y))

})