save state; test clustering not OK, all others OK

[epclust.git] / epclust / src / filter.cpp

diff --git a/epclust/src/filter.cpp b/epclust/src/filter.cpp
index ee24af6..cb5a8a0 100644 (file)
--- a/epclust/src/filter.cpp
+++ b/epclust/src/filter.cpp
@@ -1,5 +1,12 @@
 #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
@@ -10,31 +17,35 @@ using namespace Rcpp;
 //'
 //' @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_;
 }