Fixing (partially) the zero weights issue in 1-D and moving all lwls1…

…d to Rlwls1d. locfit dependancy dropped.

R/RcppExports.R

@@ -0,0 +1,23 @@
+# This file was generated by Rcpp::compileAttributes
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+dropZeroElementsXYWin <- function(win, xin, yin) {
+    .Call('tPACE_dropZeroElementsXYWin', PACKAGE = 'tPACE', win, xin, yin)
+}
+
+interp2lin <- function(xin, yin, zin, xou, you) {
+    .Call('tPACE_interp2lin', PACKAGE = 'tPACE', xin, yin, zin, xou, you)
+}
+
+Rlwls1d <- function(bw, kernel_type, win, xin, yin, xout, npoly = 1L, nder = 0L) {
+    .Call('tPACE_Rlwls1d', PACKAGE = 'tPACE', bw, kernel_type, win, xin, yin, xout, npoly, nder)
+}
+
+Rmullwlsk <- function(bw, kernel_type, tPairs, cxxn, win, xgrid, ygrid, bwCheck) {
+    .Call('tPACE_Rmullwlsk', PACKAGE = 'tPACE', bw, kernel_type, tPairs, cxxn, win, xgrid, ygrid, bwCheck)
+}
+
+Rrotatedmullwlsk <- function(bw, kernel_type, tPairs, cxxn, win, xygrid, npoly, bwCheck) {
+    .Call('tPACE_Rrotatedmullwlsk', PACKAGE = 'tPACE', bw, kernel_type, tPairs, cxxn, win, xygrid, npoly, bwCheck)
+}
+

R/cvlwls1d.R

@@ -63,13 +63,13 @@ cvlwls1d <- function(yy, t, kernel, npoly, nder, dataType ){
 
 	}
 
-        # yin = xxn[order(ttn)]
-        # xin = sort(ttn)    
         if (any(win==0)){
-           mu = lwls1d(bw= bw[j], kern=kernel, npoly=npoly, nder= nder, xin = ttn, yin= xxn, xout=out, win = win) 
+          nz = c(win != 0)
+          mu = Rlwls1d(bw= bw[j], kern=kernel, npoly=npoly, nder= nder, xin = ttn[nz], yin= xxn[nz], xout=out, win = win[nz])       
         } else {
-          mu = Rlwls1d(bw= bw[j], kern=kernel, npoly=npoly, nder= nder, xin = ttn, yin= xxn, xout=out, win = win)       
+          mu = Rlwls1d(bw= bw[j], kern=kernel, npoly=npoly, nder= nder, xin = ttn, yin= xxn, xout=out, win = win)
         }
+
         # if invalid==0 {
         cv[j]=cv[j]+t(obs-mu)%*%(obs-mu);
         count[j]=count[j]+1;

R/pc_covE.R

@@ -63,8 +63,13 @@ pc_covE = function(obsGrid, regGrid, bw_userCov, rotationCut=c(0, 1), kernel = '
   win2 = rep(1, nrow(rcovdiag))
 
   # yvar is the smoothed variance function along the diagonal line
-  yvar = lwls1d(bw = bw_userCov[1], kern = kernel, xin = rcovdiag[,1],
-    yin = rcovdiag[,2], win = win2, xout = rcutGrid, returnFit = FALSE)
+  # yvar = lwls1d(bw = bw_userCov[1], kern = kernel, xin = rcovdiag[,1],
+  #  yin = rcovdiag[,2], win = win2, xout = rcutGrid, returnFit = FALSE)
+    yvar = Rlwls1d(bw = bw_userCov[1], kern = kernel, xin = rcovdiag[,1],
+                  yin = rcovdiag[,2], win = win2, xout = rcutGrid)
+
+
+
 
   # Estimate variance of measurement error term
   # use quadratic form on diagonal to estimate Var(x(t))

src/Rlwls1d.cpp

@@ -7,7 +7,6 @@
 // [[Rcpp::depends(RcppEigen)]]
 // [[Rcpp::export]]
 
-
 Eigen::VectorXd Rlwls1d( const double & bw, const std::string kernel_type, const Eigen::Map<Eigen::VectorXd> & win, const Eigen::Map<Eigen::VectorXd> & xin, const Eigen::Map<Eigen::VectorXd> & yin, const Eigen::Map<Eigen::VectorXd> & xout, const int & npoly = 1, const int & nder = 0){
 
 

src/dropZeroElementsXYWin.cpp

@@ -0,0 +1,48 @@
+#include <RcppEigen.h>
+#include <map>          // to map kernels to integers for the switch
+#include <string>       // to read in the kernel name
+#include <vector>       // to use vectors
+#include <algorithm>    // to get the intersect and sort
+
+// [[Rcpp::depends(RcppEigen)]]
+// [[Rcpp::export]]
+
+
+Eigen::MatrixXd dropZeroElementsXYWin( const Eigen::Map<Eigen::VectorXd> & win, const Eigen::Map<Eigen::VectorXd> & xin, const Eigen::Map<Eigen::VectorXd> & yin){
+
+  const unsigned int nXGrid = xin.size();
+
+  // Check that we have equal number of readings
+  if( nXGrid != yin.size()){
+    Rcpp::stop("The input Y-grid does not have the same number of points as input X-grid.");
+  }
+
+  if( nXGrid != win.size()){
+    Rcpp::stop("The input weight vector does not have the same number of points as input X-grid.");
+  }
+
+  unsigned int nZeroElements = std::count(&win[0], &win[nXGrid], 0.);
+
+  Eigen::MatrixXd Q(nZeroElements,3);
+
+  // Check that we do not have zero weights // Should do a try-catch here
+  if( nZeroElements != 0 ){  //
+    Eigen::MatrixXd Q(nXGrid - nZeroElements,3);
+    unsigned int q = 0;
+    for( unsigned int i = 0; i != nXGrid; ++i){
+      if ( win[i] != 0 ) {
+          Q(q,0) = xin[i];
+          Q(q,1) = yin[i];
+          Q(q,2) = win[i];
+          ++q;
+      }
+    }
+    return( Q );         
+  } else {
+    Eigen::MatrixXd Q(nXGrid,3);
+    Q.col(0) = xin;
+    Q.col(1) = yin;
+    Q.col(2) = win;
+    return( Q );  
+  }
+}