#include <Rcpp.h>

using namespace Rcpp;

//' filter
//'
//' Filter [time-]series by replacing all values by the moving average of previous, current and
//' next value. Possible extensions: larger window, gaussian kernel... (but would be costly!).
//' Note: border values are unchanged.
//'
//' @param cwt Continuous wavelets transform (in columns): a matrix of size LxD
//'
//' @return The filtered CWT, in a matrix of same size (LxD)
// [[Rcpp::export]]
RcppExport SEXP epclustFilter(SEXP cwt_)
{
	// NOTE: C-style for better efficiency (this must be as fast as possible)
	int L = INTEGER(Rf_getAttrib(cwt_, R_DimSymbol))[0],
		D = INTEGER(Rf_getAttrib(cwt_, R_DimSymbol))[1];
	double *cwt = REAL(cwt_);

	SEXP fcwt_; //the filtered result
	PROTECT(fcwt_ = Rf_allocMatrix(REALSXP, L, D));
	double* fcwt = REAL(fcwt_); //pointer to the encapsulated vector

	// NOTE: unused loop parameter: shifting at the end of the loop is more efficient
	for (int col=D-1; col>=0; col++)
	{
		double v1 = cwt[0]; //first value
		double ms = v1 + cwt[1] + cwt[2]; //moving sum at second value
		for (int i=1; i<L-2; i++)
		{
			fcwt[i] = ms / 3.; //ms / 3: moving average at current index i
			ms = ms - v1 + cwt[i+2]; //update ms: remove oldest items, add next
			v1 = cwt[i]; //update first value for next iteration
		}

		// Fill a few border values not computed in the loop
		fcwt[0] = cwt[0];
		fcwt[L-2] = ms / 3.;
		fcwt[L-1] = cwt[L-1];

		// Shift by L == increment column index by 1 (storage per column)
		cwt = cwt + L;
		fcwt = fcwt + L;
	}

	UNPROTECT(1);
	return fcwt_;
}