I try to implement RCPP code for the R function below. I'm very new to C++. Below I put R function and the RCPP function. The errors occur in the C++ script lines 43,44,50,51,55.
#include <Rcpp.h>
using namespace Rcpp;
// Helper function to extract a row slice from a NumericMatrix
NumericVector extract_row_slice(const NumericMatrix& data, int row, int start, int end) {
if (end >= data.ncol()) end = data.ncol() - 1;
NumericVector slice(end - start + 1);
for (int i = start; i <= end; i++) {
slice[i - start] = data(row, i);
}
return slice;
}
// Main function to compute critical values
// [[Rcpp::export]]
NumericVector crit_val_cpp(int reps, int k, int n, double Quantile, double eps) {
NumericVector W(reps);
int sublength = n - 1;
double n_squared = n * n;
for (int j = 0; j < reps; j++) {
NumericMatrix data(k, n, Rcpp::rnorm(n * k).begin());
NumericVector substat(sublength);
NumericVector scale(sublength);
// Precompute the scaling factor for each t
for (int i = 0; i < sublength; i++) {
scale[i] = (i + 1) * (sublength - i) / n_squared;
}
for (int t = 0; t < sublength; t++) {
NumericVector mean1(k), mean2(k);
// Calculate means for each segment
for (int i = 0; i < k; i++) {
mean1[i] = mean(extract_row_slice(data, i, 0, t));
mean2[i] = mean(extract_row_slice(data, i, t + 1, n - 1));
}
// Calculate inter1 and inter2
NumericMatrix inter1(k, t + 1), inter2(k, n - t - 1);
for (int h = 0; h < k; h++) {
NumericVector temp1 = cumsum(extract_row_slice(data, h, 0, t)) - (seq_len(t + 1) - 1) * mean1[h];
NumericVector temp2 = cumsum(extract_row_slice(data, h, t + 1, n - 1)) - (seq_len(n - t - 1) - 1) * mean2[h];
std::copy(temp1.begin(), temp1.end(), inter1(h, _).begin());
std::copy(temp2.begin(), temp2.end(), inter2(h, _).begin());
}
// Compute M1 and M2 matrices
NumericMatrix M1 = inter1.t() * inter1 / n_squared;
NumericMatrix M2 = inter2.t() * inter2 / n_squared;
// Calculate substat[t] using the quadratic form
NumericVector diff = scale[t] * (mean1 - mean2);
substat[t] = n * crossprod(diff, solve(M1 + M2, diff));
}
// Store the maximum of the substatistics adjusted for eps
int lower_idx = floor(n * eps);
int upper_idx = ceil(n * (1 - eps)) - 1;
W[j] = max(abs(substat[Range(lower_idx, upper_idx)]));
}
// Return W
return W;
}
This RCPP function returns W. I would use W to get the critical values as in the R script below. Here is the R script;
library(Rcpp)
crit.val <- function(reps, k, n, Quantile, eps){
W <- rep(0,reps)
sublength <- n-1
scale <- (1:(n-1))/n^2*((n-1):1)
for (j in 1:reps){
data <- matrix( rnorm(n*k, 0, 1), k, n )
substat <- rep(0, sublength)
for (t in 1:(n-1)){
if (k==1){
mean1 <- mean( data[1, 1:t] )
mean2 <- mean( data[1, (1+t):n] )
}
if (k>1){
mean1 <- apply( matrix(data[, 1:t], ncol = t), 1, mean )
mean2 <- apply( matrix(data[, (t+1):n], ncol = n-t), 1, mean )
}
inter1 <- NULL
inter2 <- NULL
for (h in 1:k){
inter1 <- rbind(inter1, cumsum( data[h, 1:t ] ) - (1:t)*mean1[h])
inter2 <- rbind(inter2, cumsum( data[h, n:(t+1)] ) - (1:(n-t))*mean2[h])
}
M1 <- inter1 %*% t(inter1)/n^2
M2 <- inter2 %*% t(inter2)/n^2
substat[t] <- n * matrix(scale[t]*(mean1 - mean2), 1, k) %*% solve(M1+M2) %*%
matrix(scale[t]*(mean1 - mean2), k, 1)
}
W[j] <- max( abs(substat)[ floor(n*eps):ceiling(n*(1-eps)) ] )
print(j)
print(Sys.time())
}
return( quantile( W, Quantile ) )
}
crit.val(reps=100, k=1, n=100, Quantile = c(0.9,0.95,0.99,0.995,0.999), eps=0.01)
Here are the errors in RCPP script;
-Line 43 invalid operands to binary expression ('sugar::Cumsum<14, true, Vector<14>>' and 'typename traits::enable_if<traits::is_convertible<typename traits::remove_const_and_reference<double>::type, typename traits::storage_type<13>::type>::value, sugar::Times_Vector_Primitive<13, true, Minus_Vector_Primitive<13, false, SeqLen>>>::type' (aka 'Rcpp::sugar::Times_Vector_Primitive<13, true, Rcpp::sugar::Minus_Vector_Primitive<13, false, Rcpp::sugar::SeqLen>>'))
Line 44 invalid operands to binary expression ('sugar::Cumsum<14, true, Vector<14>>' and 'typename traits::enable_if<traits::is_convertible<typename traits::remove_const_and_reference<double>::type, typename traits::storage_type<13>::type>::value, sugar::Times_Vector_Primitive<13, true, Minus_Vector_Primitive<13, false, SeqLen>>>::type' (aka 'Rcpp::sugar::Times_Vector_Primitive<13, true, Rcpp::sugar::Minus_Vector_Primitive<13, false, Rcpp::sugar::SeqLen>>'))
-Line 50 no member named 't' in 'Rcpp::Matrix<14>'
-Line 51no member named 't' in 'Rcpp::Matrix<14>'
-Line 55use of undeclared identifier 'solve'
I initially created a row_slicer function in Rcpp to handle specific errors successfully. However, I've encountered subsequent errors that I'm struggling to resolve, despite attempting to use various troubleshooting methods, including AI tools. My background is primarily in R, and I'm relatively new to C++, which is why I'm finding it challenging to address these issues independently.
I would greatly appreciate any guidance, comments, or assistance in troubleshooting and integrating my Rcpp code into a functional R function.
Thank you in advance for any help!
Regarding lines 43,44: This is because the operator '-' is defined for NumericVector, but not for the type outputted by the cumsum function (Rcpp::sugar::Cumsum<14, true, Rcpp::NumericVector>). In this case it isn't converted implicitly, so you have to do it, like this:
NumericVector temp1 = as<NumericVector>(wrap(cumsum(extract_row_slice(data, h, 0, t)))) - as<NumericVector>(wrap(seq_len(t + 1) - 1)) * mean1[h];
NumericVector temp2 = as<NumericVector>(wrap(cumsum(extract_row_slice(data, h, t + 1, n - 1)))) - as<NumericVector>(wrap(seq_len(n - t - 1) - 1)) * mean2[h];
Bit messy but there you go, implicit conversion has limits.
Regarding your other issues, it looks like you are trying to solve problems better suited to RcppArmadillo. You can transpose and solve matrixes using this library. See:
https://arma.sourceforge.net/docs.html#solve https://github.com/RcppCore/RcppArmadillo
Rcpp matrixes simply do not have this functionality. I hope this helps.