diff --git a/DESCRIPTION b/DESCRIPTION index 026f154..2b2481f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: covafillr Title: Local Polynomial Regression of State Dependent Covariates in State-Space Models -Version: 0.4.1-0001 -Date: 2017-05-04 +Version: 0.4.2 +Date: 2017-11-20 Authors@R: person(c("Christoffer","Moesgaard"), "Albertsen", email = "cmoe@aqua.dtu.dk", role = c("aut", "cre")) Maintainer: Christoffer Moesgaard Albertsen diff --git a/NEWS b/NEWS index 7727f41..df09cd9 100644 --- a/NEWS +++ b/NEWS @@ -1,5 +1,10 @@ Covafillr News +CHANGES IN VERSION 0.4.2 + BUG FIXES: + - Rcpp:::Rcpp.plugin.maker changed to Rcpp::Rcpp.plugin.maker + - Fixed index bug in tricubic interpolation + CHANGES IN VERSION 0.4.1 NEW FEATURES: diff --git a/R/RcppPlugin.R b/R/RcppPlugin.R index bcc8423..85c22b3 100644 --- a/R/RcppPlugin.R +++ b/R/RcppPlugin.R @@ -1,6 +1,6 @@ -inlineCxxPlugin <- Rcpp:::Rcpp.plugin.maker( +inlineCxxPlugin <- Rcpp::Rcpp.plugin.maker( include.before = paste("#include", "#include", "", diff --git a/inst/examples/ggplot/ggplotexample.R b/inst/examples/ggplot/ggplotexample.R index 2c742d6..3781665 100644 --- a/inst/examples/ggplot/ggplotexample.R +++ b/inst/examples/ggplot/ggplotexample.R @@ -17,5 +17,5 @@ cf$predict(0,TRUE) ggplot(aes(x=x,y=y),data=data) + geom_point() + - stat_covafill(bandwith=1,polyDegree=3L) + - geom_line(aes(x=x,y=y),data=data.frame(x=seq(-2,2,len=1000),y=pol(seq(-2,2,len=1000))),col="red") + stat_covafill(bandwith=1,polyDegree=3L)+ +geom_label(aes(x=x,y=y,label=label,hjust=0),data=data.frame(x=0,y=-25,label="ggplot(aes(x=x,y=y),data=data) +\n\tgeom_point() +\n\tstat_covafill(bandwith=1,polyDegree=3L)")) diff --git a/inst/examples/interpolate/interpolate.R b/inst/examples/interpolate/interpolate.R new file mode 100644 index 0000000..81a02a4 --- /dev/null +++ b/inst/examples/interpolate/interpolate.R @@ -0,0 +1,41 @@ +library(TMB) +library(covafillr) + +compile("interpolate.cpp",CXXFLAGS = paste(cxxFlags(),"-std=c++14")) + +dyn.load(dynlib("interpolate")) + +coord <- as.matrix(expand.grid(seq(-1,1,len=5),seq(-1,1,len=5))) +ftrue <- function(x)sum(x^2) +covObsTrue <- apply(coord,1,function(x)ftrue(x) )##+ rnorm(1,0,0.0)) + +n <- 5000 +X <- rbind(runif(n,-1,1), + runif(n,-1,1)) +obs <- apply(X,2,ftrue) + rnorm(n,0,0.1) + +dat <- list(coord = coord, + p = 2, + h = 0, + d = 2, + obs = obs, + x = X, + repx = seq(-1,1,0.05), + repy = seq(-1,1,0.05)) +dat$h <- covafillr:::suggestBandwith(coord,dat$p) +pars <- list(covObs = rep(0,nrow(coord)), + logSd = 0) + +obj <- MakeADFun(data = dat, + parameters = pars, + DLL = "interpolate") + +opt <- nlminb(obj$par,obj$fn,obj$gr,obj$he,control=list(iter.max=1000,eval.max=1000)) + +opt + + +rp <- obj$report() + +contour(dat$repx,dat$repy,rp$predXY) +contour(dat$repx,dat$repy,outer(dat$repx,dat$repy,Vectorize(function(x,y)ftrue(c(x,y)))),col="red",add=TRUE) diff --git a/inst/examples/interpolate/interpolate.cpp b/inst/examples/interpolate/interpolate.cpp new file mode 100644 index 0000000..05519ca --- /dev/null +++ b/inst/examples/interpolate/interpolate.cpp @@ -0,0 +1,38 @@ + +#include +#include + +template +Type objective_function::operator() () +{ + DATA_MATRIX(coord); + DATA_INTEGER(p); + DATA_VECTOR(h); + DATA_SCALAR(d); + DATA_VECTOR(obs); + DATA_MATRIX(x); + DATA_VECTOR(repx); + DATA_VECTOR(repy); + PARAMETER_VECTOR(covObs); + PARAMETER(logSd); + covafill cf(coord,covObs,h,p); + covatree ct(d,&cf); + Type nll = 0.0; + for(int i = 0; i < obs.size(); ++i){ + vector xtmp = x.col(i); + vector ctOut = ct(xtmp); + nll -= dnorm(obs(i),ctOut(0),exp(logSd),true); + } + + matrix predXY(repx.size(),repy.size()); + predXY.setZero(); + for(int i = 0; i < repx.size(); ++i) + for(int j = 0; j < repx.size(); ++j){ + vector tmp(2); + tmp(0) = repx(i); tmp(1) = repy(j); + predXY(i,j) = cf(tmp)(0); + } + + REPORT(predXY); + return nll; +} diff --git a/inst/examples/interpolate_ssm/interpolate_ssm.R b/inst/examples/interpolate_ssm/interpolate_ssm.R new file mode 100644 index 0000000..a38de57 --- /dev/null +++ b/inst/examples/interpolate_ssm/interpolate_ssm.R @@ -0,0 +1,51 @@ +library(TMB) +library(covafillr) + +compile("interpolate_ssm.cpp",CXXFLAGS = paste(cxxFlags(),"-std=c++14")) + +dyn.load(dynlib("interpolate_ssm")) + + +coord <- as.matrix(expand.grid(seq(-1,1,len=5),seq(-1,1,len=5))) +ftrue <- function(x) sum(diag(1,2) %*% (0.5 * x * x - c(1,-1) * x)) +covObsTrue <- apply(coord,1,function(x)ftrue(x) )##+ rnorm(1,0,0.0)) + +n <- 1000 +dt <- 1/100 +Xt <- matrix(0,2,n*(1/dt)) +Xt[,1] <- 0.5 +for(i in 2:ncol(Xt)){ + v <- numDeriv::grad(ftrue,Xt[,i-1]) + Xt[,i] <- Xt[,i-1] + -v*dt + rnorm(2,0,0.1 * sqrt(dt)) +} +X <- Xt[,(1:ncol(Xt))*dt == round((1:ncol(Xt))*dt)] +plot(t(X),type="l") +plot(X[1,],type="l") + +obs <- apply(X,2,ftrue) + rnorm(n,0,0.1) + +dat <- list(coord = coord, + p = 2, + h = 0, + d = 2, + obs = obs, + x = X, + repx = seq(-1,1,0.05), + repy = seq(-1,1,0.05)) +dat$h <- covafillr:::suggestBandwith(coord,dat$p) +pars <- list(covObs = rep(0,nrow(coord)), + logSd = 0) + +obj <- MakeADFun(data = dat, + parameters = pars, + DLL = "interpolate") + +opt <- nlminb(obj$par,obj$fn,obj$gr,obj$he,control=list(iter.max=1000,eval.max=1000)) + +opt + + +rp <- obj$report() + +contour(dat$repx,dat$repy,rp$predXY) +contour(dat$repx,dat$repy,outer(dat$repx,dat$repy,Vectorize(function(x,y)ftrue(c(x,y)))),col="red",add=TRUE) diff --git a/inst/examples/interpolate_ssm/interpolate_ssm.cpp b/inst/examples/interpolate_ssm/interpolate_ssm.cpp new file mode 100644 index 0000000..86f5bb4 --- /dev/null +++ b/inst/examples/interpolate_ssm/interpolate_ssm.cpp @@ -0,0 +1,38 @@ + +#include +#include + +template +Type objective_function::operator() () +{ + DATA_MATRIX(coord); + DATA_INTEGER(p); + DATA_VECTOR(h); + DATA_SCALAR(d); + DATA_VECTOR(obs); + DATA_MATRIX(x); + DATA_VECTOR(repx); + DATA_VECTOR(repy); + PARAMETER_VECTOR(covObs); + PARAMETER(logSd); + covafill cf(coord,covObs,h,p); + covatree ct(d,&cf); + Type nll = 0.0; + for(int i = 0; i < obs.size(); ++i){ + vector xtmp = x.col(i); + vector ctOut = cf(xtmp); + nll -= dnorm(obs(i),ctOut(0),exp(logSd),true); + } + + matrix predXY(repx.size(),repy.size()); + predXY.setZero(); + for(int i = 0; i < repx.size(); ++i) + for(int j = 0; j < repx.size(); ++j){ + vector tmp(2); + tmp(0) = repx(i); tmp(1) = repy(j); + predXY(i,j) = cf(tmp)(0); + } + + REPORT(predXY); + return nll; +} diff --git a/inst/examples/tmbtest/tmbtest.R b/inst/examples/tmbtest/tmbtest.R index df44720..1b54e9e 100644 --- a/inst/examples/tmbtest/tmbtest.R +++ b/inst/examples/tmbtest/tmbtest.R @@ -2,11 +2,56 @@ library(numDeriv) library(TMB) library(covafillr) -system("touch tmbtest.cpp") TMB::compile("tmbtest.cpp",CXXFLAGS=paste(cxxFlags(),"-std=c++14")) dyn.load(dynlib("tmbtest")) + + + +## Test one dim +coord <- as.matrix(expand.grid(seq(-1,1,0.001))) +ftrue <- function(x)sum(x^3) +covObs <- apply(coord,1,function(x)ftrue(x) + rnorm(1,0,0.01)) + +plot(coord,covObs) + +dat <- list(coord = coord, + covObs = covObs, + p = 3, + h = 0, + d = 100) +dat$h <- suggestBandwith(dat$coord,dat$p) + +obj <- MakeADFun(data = dat, + parameters = list(x = c(0)), + DLL = "tmbtest") + +fdiff <- Vectorize(function(x) obj$f(x) - x^3) +fdiff(seq(-1,1,0.01)) + + +gdiff <- Vectorize(function(x) obj$gr(x) - 3*x^2) +gdiff(seq(-1,1,0.01)) + +cf <- covafill(dat$coord,dat$covObs,h=dat$h,p=as.integer(dat$p)) +x0 <- seq(-1,1,0.1) +pr <- cf$predict(x0) + +pr[,1]-x0^3 +pr[,1]-2*x0^2 + + +g1 <- 3*x0^2 +g2 <- sapply(x0,obj$gr) +g3 <- pr[,2] + +plot(x0,g1) +lines(x0,g2,col="red") +lines(x0,g3,col="blue") + +## Test two dim + coord <- as.matrix(expand.grid(seq(-10,10,0.2),seq(-10,10,0.2))) ftrue <- function(x)sum(x^3) covObs <- apply(coord,1,function(x)ftrue(x) + rnorm(1,0,0.01)) @@ -27,6 +72,7 @@ obj$fn(c(0,1)) obj$gr() obj$he() + numDeriv::grad(obj$fn,c(0,1)) numDeriv::hessian(obj$fn,c(0,1)) diff --git a/inst/examples/tmbtest/tmbtest.cpp b/inst/examples/tmbtest/tmbtest.cpp index 7a9b116..2a1ded7 100644 --- a/inst/examples/tmbtest/tmbtest.cpp +++ b/inst/examples/tmbtest/tmbtest.cpp @@ -10,8 +10,11 @@ Type objective_function::operator() () DATA_VECTOR(covObs); DATA_INTEGER(p); DATA_VECTOR(h); + DATA_SCALAR(d); PARAMETER_VECTOR(x); covafill cf(coord,covObs,h,p); - Type val = evalFill((CppAD::vector)x, cf)[0]; + covatree ct(d,&cf); + //Type val = evalFill((CppAD::vector)x, cf)[0]; + Type val = ct(x)(0); return val; }