#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
//'
//' Filter time-series
//'
//' @param cwt Continuous wavelets transform
//'
//' @return The filtered CWT
// [[Rcpp::export]]
RcppExport SEXP epclustFilter(SEXP cwt_)
{
	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[0];
		double ma = v1 + cwt[1] + cwt[2];
		for (int i=1; i<L-2; i++)
		{
			fcwt[i] = ma / 3.;
			ma = ma - v1 + cwt[i+2];
			v1 = cwt[i];
		}
		fcwt[0] = cwt[0];
		fcwt[L-2] = ma / 3.;
		fcwt[L-1] = cwt[L-1];
		cwt = cwt + L;
		fcwt = fcwt + L;
	}

//	REAL(fcwt_) = fcwt;
	UNPROTECT(1);

	return fcwt_;
}