1

I have been experimenting with the RcppArrayFire Package, mostly rewriting some cost functions from RcppArmadillo and can't seem to get over "no viable conversion from 'af::array' to 'float'. I have also been getting some backend errors, the example below seems free of these.

This cov-var example is written poorly just to use all relevant coding pieces from my actual cost function. As of now it is the only addition in a package generated by, "RcppArrayFire.package.skeleton".

#include "RcppArrayFire.h"
#include <Rcpp.h>
// [[Rcpp::depends(RcppArrayFire)]]
// [[Rcpp::export]]
float example_ols(const RcppArrayFire::typed_array<f32>& X_vect,  const RcppArrayFire::typed_array<f32>& Y_vect){

  int Len = X_vect.dims()[0];
  int Len_Y = Y_vect.dims()[0];

  while( Len_Y < Len){
    Len --;
  }

  float mean_X = af::sum(X_vect)/Len;
  float mean_Y = af::sum(Y_vect)/Len;

  RcppArrayFire::typed_array<f32> temp(Len);
  RcppArrayFire::typed_array<f32> temp_x(Len);

  for( int f = 0; f < Len; f++){
    temp(f) = (X_vect(f) - mean_X)*(Y_vect(f) - mean_Y);
    temp_x(f) = af::pow(X_vect(f) -mean_X, 2);
  }

  return af::sum(temp)/af::sum(temp_x);
}


/*** R
  X <- 1:10
  Y <- 2*X +rnorm(10, mean = 0, sd = 1)
  example_ols(X, Y)
*/
Ralf Stubner
  • 26,263
  • 3
  • 40
  • 75
skatz
  • 115
  • 7

1 Answers1

3

The first thing to consider is the af::sum function, which comes in different forms: An sf::sum(af::array) that returns an af::array in device memory and a templated af::sum<T>(af::array) that returns a T in host memory. So the minimal change to your example would be using af::sum<float>:

#include "RcppArrayFire.h"
#include <Rcpp.h>
// [[Rcpp::depends(RcppArrayFire)]]
// [[Rcpp::export]]
float example_ols(const RcppArrayFire::typed_array<f32>& X_vect,
                  const RcppArrayFire::typed_array<f32>& Y_vect){

    int Len = X_vect.dims()[0];
    int Len_Y = Y_vect.dims()[0];

    while( Len_Y < Len){
        Len --;
    }

    float mean_X = af::sum<float>(X_vect)/Len;
    float mean_Y = af::sum<float>(Y_vect)/Len;

    RcppArrayFire::typed_array<f32> temp(Len);
    RcppArrayFire::typed_array<f32> temp_x(Len);

    for( int f = 0; f < Len; f++){
        temp(f) = (X_vect(f) - mean_X)*(Y_vect(f) - mean_Y);
        temp_x(f) = af::pow(X_vect(f) -mean_X, 2);
    }

    return af::sum<float>(temp)/af::sum<float>(temp_x);
}


/*** R
set.seed(1)
X <- 1:10
Y <- 2*X +rnorm(10, mean = 0, sd = 1)
example_ols(X, Y)
*/

However, there are more things one can improve. In no particular order:

  • You don't need to include Rcpp.h.
  • There is an af::mean function for computing the mean of an af::array.
  • In general RcppArrayFire::typed_array<T> is only needed for getting arrays from R into C++. Within C++ and for the way back you can use af::array.
  • Even when your device does not support double, you can still use double values on the host.
  • In order to get good performance, you should avoid for loops and use vectorized functions, just like in R. You have to impose equal dimensions for X and Y, though.

Interestingly I get a different result when I use vectorized functions. Right now I am not sure why this is the case, but the following form makes more sense to me. You should verify that the result is what you want to get:

#include <RcppArrayFire.h>
// [[Rcpp::depends(RcppArrayFire)]]
// [[Rcpp::export]]
double example_ols(const RcppArrayFire::typed_array<f32>& X_vect,
                   const RcppArrayFire::typed_array<f32>& Y_vect){

    double mean_X = af::mean<double>(X_vect);
    double mean_Y = af::mean<double>(Y_vect);

    af::array temp = (X_vect - mean_X) * (Y_vect - mean_Y);
    af::array temp_x = af::pow(X_vect - mean_X, 2.0);

    return af::sum<double>(temp)/af::sum<double>(temp_x);
}

/*** R
set.seed(1)
X <- 1:10
Y <- 2*X +rnorm(10, mean = 0, sd = 1)
example_ols(X, Y)
*/ 

BTW, an even shorter version would be:

#include <RcppArrayFire.h>
// [[Rcpp::depends(RcppArrayFire)]]
// [[Rcpp::export]]
af::array example_ols(const RcppArrayFire::typed_array<f32>& X_vect,
                      const RcppArrayFire::typed_array<f32>& Y_vect){

    return af::cov(X_vect, Y_vect) / af::var(X_vect);
}

Generally it is a good idea to use the in-build functions as much as possible.

Ralf Stubner
  • 26,263
  • 3
  • 40
  • 75
  • 1
    I will add this bit incase anyone else has the same problem. Not necessarily needed for this example but useful. To extract an element from an af::array *array_name(index_number).host() i.e. in this example something like *X_vect(0).host() . This is needed for stuff like logical comparisons – skatz Dec 06 '19 at 21:53