#include <Rcpp.h>
+#include <R.h> // Rprintf()
+//#include <R_ext/Utils.h> // user interrupts
+#include <Rdefines.h>
+#include <Rinternals.h>
+
+#include <stdio.h>
+
using namespace Rcpp;
//' filter
//'
//' @return The filtered CWT
// [[Rcpp::export]]
-NumericMatrix epclustFilter(NumericMatrix cwt)
+RcppExport SEXP epclustFilter(SEXP cwt_)
{
- int L = cwt.nrow(),
- D = cwt.ncol();
- NumericMatrix fcwt(L, D); //fill with 0... TODO: back to SEXP C-style?
- double *cwt_c = cwt.begin(),
- *fcwt_c = fcwt.begin();
+ int L = INTEGER(Rf_getAttrib(cwt_, R_DimSymbol))[0],
+ D = INTEGER(Rf_getAttrib(cwt_, R_DimSymbol))[1];
+ double *cwt = REAL(cwt_);
+ SEXP fcwt_;
+ PROTECT(fcwt_ = Rf_allocMatrix(REALSXP, L, D));
+ double* fcwt = REAL(fcwt_); //(double*)malloc(L*D*sizeof(double));
//TODO: coding style is terrible... no time for now.
for (int col=0; col<D; col++)
{
- double v1 = cwt_c[0];
- double ma = v1 + cwt[1] + cwt_c[2];
+ double v1 = cwt[0];
+ double ma = v1 + cwt[1] + cwt[2];
for (int i=1; i<L-2; i++)
{
- fcwt_c[i] = ma / 3.;
- ma = ma - v1 + cwt_c[i+2];
- v1 = cwt_c[i];
+ fcwt[i] = ma / 3.;
+ ma = ma - v1 + cwt[i+2];
+ v1 = cwt[i];
}
- fcwt_c[0] = cwt_c[0];
- fcwt_c[L-2] = ma / 3.;
- fcwt_c[L-1] = cwt_c[L-1];
- cwt_c = cwt_c + L;
- fcwt_c = fcwt_c + L;
+ fcwt[0] = cwt[0];
+ fcwt[L-2] = ma / 3.;
+ fcwt[L-1] = cwt[L-1];
+ cwt = cwt + L;
+ fcwt = fcwt + L;
}
- return fcwt;
+// REAL(fcwt_) = fcwt;
+ UNPROTECT(1);
+
+ return fcwt_;
}