setwd("C:/Users/kenrice/Desktop/SISG-Advanced") # before starting, download nnfind.c, to your working directory (and read it) file.show("nnfind.c") file.show("bin_heap.c") # # at the Command Line; # cd C:\Users\kenrice\Desktop\SISG-Advanced # remove old versions del nnfind.o bin_heap.o nnfind.dll # notice that there was a bunch of path changes at Rtools installation # (run as Administrator, on Windows) # but we need one more; path C:\Program Files\R\R-3.4.1\bin\x64;%PATH% R CMD SHLIB nnfind.c bin_heap.c # May also need to force 64-bit architecture # R --arch x64 CMD SHLIB nnfind.c bin_heap.c # back to R setwd("C:/Users/kenrice/Desktop/SISG-Advanced") dyn.load("nnfind.dll") nnfind<-function(x, scale=FALSE){ if (scale) x<-scale(x) else x<-as.matrix(x) n<-nrow(x) p<-ncol(x) rval<-.C("within_neighbours",as.double(x), as.integer(n), as.integer(p), integer(n), double(n)) rval[[4]]+1 # the fourth element, plus one to "translate" the indexing } data(faithful) # built-in geyser dataset faithful dim(faithful) names(faithful) plot(faithful) par(mfrow=c(1,2)) nnoutput <- nnfind(faithful) plot(faithful) with(faithful, segments( x0=eruptions, y0=waiting, x1=eruptions[nnoutput], y1=waiting[nnoutput]) ) nnoutputs <- nnfind(faithful, scale=TRUE) plot(faithful) with(faithful, segments( x0=eruptions, y0=waiting, x1=eruptions[nnoutputs], y1=waiting[nnoutputs]) ) between <- function(x,y){ rval<-.C("between_neighbours",as.double(x), nrow(x),as.double(y),nrow(y), ncol(x), neighbour=integer(nrow(y)), double(nrow(y))) rval$neighbour+1 } # a "fun" example; drawing lines from where you click to the nearest point f2<-scale(faithful) par(mfrow=c(1,1)) plot(f2, main="Finding nearest point; escape to exit") repeat({ click <- t(as.matrix(unlist(locator(1)))) nearest <- between(as.matrix(f2),click) segments(click[,1],click[,2],f2[nearest,1],f2[nearest,2],col="red") }) dyn.unload("nnfind.dll") # another example; uniform points in space # how far are they from their nearest neighbor? # v similar function to above extracts the distances dyn.load("nnfind.dll") # if needed dists.to.neeb <-function(x, scale=FALSE){ if (scale) x<-scale(x) else x<-as.matrix(x) n<-nrow(x) p<-ncol(x) rval<-.C("within_neighbours",as.double(x), as.integer(n), as.integer(p), integer(n), double(n)) rval[[5]] } set.seed(1314) par(mfrow=c(1,2)) p <- 2 n <- 1000 X2 <- matrix(runif(n*p), n, p) plot(X2, asp=1, col="#0000FF10") d2 <- dists.to.neeb(X2) plot( density(d2) ) lines( density(dists.to.neeb(matrix(runif(n*p), n, p)) ) , lty=2) lines( density(dists.to.neeb(matrix(runif(n*p), n, p)) ) , lty=2) lines( density(dists.to.neeb(matrix(runif(n*p), n, p)) ) , lty=2) p <- 3 n <- 1000 X3 <- matrix(runif(n*p), n, p) library("scatterplot3d") scatterplot3d(X3[,1] ,X3[,2], X3[,3]) # yuk! d3 <- dists.to.neeb(X3) plot( density(d2), xlim=c(0,0.1)) lines(density(d3), col="red") # further over to the right, less peaked lines( density(dists.to.neeb(matrix(runif(n*p), n, p)) ) , lty=2, col=2) lines( density(dists.to.neeb(matrix(runif(n*p), n, p)) ) , lty=2, col=2) lines( density(dists.to.neeb(matrix(runif(n*p), n, p)) ) , lty=2, col=2) p <- 100 n <- 1000 X100 <- matrix(runif(n*p), n, p) d100 <- dists.to.neeb(X100) par(mfrow=c(1,1)) plot(density(d100)) # distances way bigger lines( density(dists.to.neeb(matrix(runif(n*p), n, p)) ) , lty=2) plot(density(scale(d100)), xlim=c(-4,4)) lines(density(scale(d3)), col="red") lines(density(scale(d2)), col="green") qqnorm(scale(d2), asp=1, ylim=range(c(scale(d2),scale(d3),scale(d100)))) points(qqnorm(scale(d3), plot.it=F), col="red") points(qqnorm(scale(d100), plot.it=F), col="green") abline(0,1,lty=2) dyn.unload("nnfind.dll") ### ###Q2 ### setwd("C:/Users/kenrice/Desktop/SISG-ADV") # First, do everything in R # The code is from the question; boxcar.R <- function(Y, X, radius, n=length(Y)){ y.smooth <- rep(0,n) x <- 0 for (i in 1:n){ count <- 0 x <- X[i] for (j in 1:n){ if(abs(X[j]-x) \n void boxcar ( int * n, double * x, double * y, double * radius, double * ysmooth ) {", codeSq,"}") file.show("boxcar6.c") Sys.getenv("path") opath <- Sys.getenv("path") opath npath <- paste("C:\\Program Files\\R\\R-3.4.1\\bin\\x64;", opath, sep="") Sys.setenv("path"=npath) system("R CMD SHLIB boxcar6.c", show.output.on.console = TRUE) dyn.load("boxcar6.dll") time3 <- system.time({ boxcar6.out <- .C("boxcar", n=length(x.test), x=x.test, y=y.test, radius=0.5 , ysmooth=rep(0, n) ) }) time1/time3 # even blimier str(boxcar6.out) plot(y.test~x.test, pch=19, col="gray70") lines(x.test, y=boxcar6.out$ysmooth, col="green", lwd=8, lty=1) lines(x.test, bc.test, col="red", lwd=2) lines(x.test, y=bc.test2, col="blue", lwd=2, lty=2) dyn.unload("boxcar6.dll") # particularly if you want to re-make it file.remove("boxcar6.dll") ## Rcpp versions # first the convolve example - as used in their vignette vignette("Rcpp-introduction") src <- " Rcpp::NumericVector xa(a); Rcpp::NumericVector xb(b); int n_xa = xa.size(), n_xb = xb.size(); Rcpp::NumericVector xab(n_xa + n_xb - 1); for (int i = 0; i < n_xa; i++) for (int j = 0; j < n_xb; j++) xab[i + j] += xa[i] * xb[j]; return xab; " fun <- cxxfunction(signature(a = "numeric", b = "numeric"), src, plugin = "Rcpp") fun(1:3, 1:4) plot( x=2:12, y=fun(rep(1/6,6), rep(1/6,6)), type="h") # same ideas but for the boxcar filter sigSq <- signature(x="numeric", y="numeric", radius="numeric") src <- " Rcpp::NumericVector zx(x); Rcpp::NumericVector zy(y); Rcpp::NumericVector zradius(radius); int n_x = zx.size(); Rcpp::NumericVector ysmooth(n_x); int count = 0; for (int i = 0; i < n_x; i++){ count = 0; ysmooth[i] = 0.0; for (int j = 0; j < n_x; j++){ if (fabs(zx[i]-zx[j]) < zradius[1]){ count++; ysmooth[i]+=zy[j]; } } ysmooth[i] = ysmooth[i]/count; } return ysmooth; } " fun <- cxxfunction(sigSq , src, plugin = "Rcpp")