af::array conversion to float or double

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)
*/

Solution

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.