Skip to content

Commit

Permalink
Fix RcppPlugin request from CRAN + new version
Browse files Browse the repository at this point in the history
  • Loading branch information
calbertsen committed Nov 20, 2017
1 parent 9cf4d21 commit 30180a7
Show file tree
Hide file tree
Showing 10 changed files with 229 additions and 7 deletions.
4 changes: 2 additions & 2 deletions 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 <cmoe@aqua.dtu.dk>
Expand Down
5 changes: 5 additions & 0 deletions 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:
Expand Down
2 changes: 1 addition & 1 deletion R/RcppPlugin.R
@@ -1,6 +1,6 @@


inlineCxxPlugin <- Rcpp:::Rcpp.plugin.maker(
inlineCxxPlugin <- Rcpp::Rcpp.plugin.maker(
include.before = paste("#include<RcppEigen.h>",
"#include<covafill/Tree>",
"",
Expand Down
4 changes: 2 additions & 2 deletions inst/examples/ggplot/ggplotexample.R
Expand Up @@ -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)"))
41 changes: 41 additions & 0 deletions 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)
38 changes: 38 additions & 0 deletions inst/examples/interpolate/interpolate.cpp
@@ -0,0 +1,38 @@

#include <TMB.hpp>
#include <covafill/TMB>

template<class Type>
Type objective_function<Type>::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<Type> cf(coord,covObs,h,p);
covatree<Type> ct(d,&cf);
Type nll = 0.0;
for(int i = 0; i < obs.size(); ++i){
vector<Type> xtmp = x.col(i);
vector<Type> ctOut = ct(xtmp);
nll -= dnorm(obs(i),ctOut(0),exp(logSd),true);
}

matrix<Type> predXY(repx.size(),repy.size());
predXY.setZero();
for(int i = 0; i < repx.size(); ++i)
for(int j = 0; j < repx.size(); ++j){
vector<Type> tmp(2);
tmp(0) = repx(i); tmp(1) = repy(j);
predXY(i,j) = cf(tmp)(0);
}

REPORT(predXY);
return nll;
}
51 changes: 51 additions & 0 deletions 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)
38 changes: 38 additions & 0 deletions inst/examples/interpolate_ssm/interpolate_ssm.cpp
@@ -0,0 +1,38 @@

#include <TMB.hpp>
#include <covafill/TMB>

template<class Type>
Type objective_function<Type>::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<Type> cf(coord,covObs,h,p);
covatree<Type> ct(d,&cf);
Type nll = 0.0;
for(int i = 0; i < obs.size(); ++i){
vector<Type> xtmp = x.col(i);
vector<Type> ctOut = cf(xtmp);
nll -= dnorm(obs(i),ctOut(0),exp(logSd),true);
}

matrix<Type> predXY(repx.size(),repy.size());
predXY.setZero();
for(int i = 0; i < repx.size(); ++i)
for(int j = 0; j < repx.size(); ++j){
vector<Type> tmp(2);
tmp(0) = repx(i); tmp(1) = repy(j);
predXY(i,j) = cf(tmp)(0);
}

REPORT(predXY);
return nll;
}
48 changes: 47 additions & 1 deletion inst/examples/tmbtest/tmbtest.R
Expand Up @@ -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))
Expand All @@ -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))

Expand Down
5 changes: 4 additions & 1 deletion inst/examples/tmbtest/tmbtest.cpp
Expand Up @@ -10,8 +10,11 @@ Type objective_function<Type>::operator() ()
DATA_VECTOR(covObs);
DATA_INTEGER(p);
DATA_VECTOR(h);
DATA_SCALAR(d);
PARAMETER_VECTOR(x);
covafill<Type> cf(coord,covObs,h,p);
Type val = evalFill((CppAD::vector<Type>)x, cf)[0];
covatree<Type> ct(d,&cf);
//Type val = evalFill((CppAD::vector<Type>)x, cf)[0];
Type val = ct(x)(0);
return val;
}

0 comments on commit 30180a7

Please sign in to comment.