Permalink
Browse files

Profile script changes

    collect more details
    Makefile: option to profile all examples
    New example ar1_4D
  • Loading branch information...
kaskr committed Feb 19, 2014
1 parent 4e7c2ad commit ce4520040f5714e4b965c4306b6bff18ad493b5e
View
@@ -2,6 +2,7 @@ rfiles = $(basename $(wildcard *.R))
cppfiles = $(basename $(wildcard *.cpp))
intersection= $(filter $(rfiles), $(cppfiles))
outputfiles = $(intersection:=.output.RData)
+profiletargets = $(intersection:=.profile)
%.output.RData : %.R %.cpp
unset MAKEFLAGS; example=$(basename $<) R --vanilla < unittest.R
@@ -14,3 +15,6 @@ clean :
%.profile : %.R %.cpp
example=$(basename $<) R --vanilla < profiler.R
+
+profile_all: $(outputfiles) $(profiletargets)
+
View
@@ -0,0 +1,54 @@
+require(TMB)
+compile("ar1_4D.cpp")
+set.seed(123)
+n <- 8 ## Size of problem = n*n
+
+## ======================= Simulate separable 2D GMRF
+## - With exponential correlation in both directions
+## - phi1 = 1-lag correlation in 1st direction
+## - phi2 = 1-lag correlation in 2nd direction
+ar1corr <- function(n,phi){
+ phi^abs(outer(1:n,1:n,"-"))
+}
+simgmrf4 <- function(n,phi){
+ dim <- c(n,n,n,n)
+ u <- array(rnorm(prod(dim)),dim)
+ L <- t(chol(ar1corr(n,phi)))
+ ar2mat <- function(x)matrix(x,nrow(x))
+ for(i in 1:4){
+ u[] <- L%*%ar2mat(u)
+ u <- aperm(u,c(4,1,2,3))
+ }
+ u
+}
+
+## ======================= Simulate data
+phi=exp(-1/(.2*n)) ## Correlation range=20% of grid size second dimension
+eta <- simgmrf4(n,phi)
+N <- rpois(length(eta),exp(eta))
+##d <- expand.grid(x=factor(1:n),y=factor(1:n))
+##d$N <- N
+
+## ======================= Parameterization of phi
+f <- function(x) 2/(1 + exp(-2 * x)) - 1
+invf <- function(y) -0.5 * log(2/(y + 1) - 1)
+
+## ======================= Fit model
+dyn.load(dynlib("ar1_4D"))
+obj <- MakeADFun(data=list(N=N),
+ parameters=list(
+ eta=array(0,c(n,n,n,n)),
+ transf_phi=invf(0.5)
+ ),
+ random=c("eta")
+ )
+runSymbolicAnalysis(obj)
+obj$control <- list(trace=1,parscale=c(1)*1e-2,REPORT=1,reltol=1e-12)
+obj$hessian <- TRUE
+##Rprof();opt <- do.call("optim",obj);Rprof(NULL)
+newtonOption(smartsearch=FALSE)
+system.time(opt <- do.call("optim",obj))
+##obj$gr()
+#c(phi1,phi2)
+#f(opt$par)
+rep <- sdreport(obj)
View
@@ -0,0 +1,27 @@
+// Separable covariance on 4D lattice with AR1 structure in each direction.
+#include <TMB.hpp>
+
+/* Parameter transform */
+template <class Type>
+Type f(Type x){return Type(2)/(Type(1) + exp(-Type(2) * x)) - Type(1);}
+
+template<class Type>
+Type objective_function<Type>::operator() ()
+{
+ DATA_VECTOR(N)
+ PARAMETER_ARRAY(eta);
+ PARAMETER(transf_phi); /* fastest running dim */
+ Type phi=f(transf_phi);
+ ADREPORT(phi);
+
+ using namespace density;
+ Type res=0;
+
+ res+=AR1(phi,AR1(phi,AR1(phi,AR1(phi))))(eta);
+
+ // logdpois = N log lam - lam
+ for(int i=0;i<N.size();i++)res-=N[i]*eta[i]-exp(eta[i]);
+
+ return res;
+
+}
View
@@ -6,9 +6,40 @@ example <- Sys.getenv("example")
## ===============================================================
Rexe <- paste(Sys.getenv("R_HOME"),"bin/exec/R",sep="/")
+
cmd <- paste("amplxe-cl -collect hotspots -result-dir ",example,".profile"," -- ",Rexe," --vanilla < ",example,".R",sep="")
system(cmd)
-cmd <- paste("amplxe-cl -report hotspots -result-dir ",example,".profile"," > ",example,".profile.txt",sep="")
-system(cmd)
-df <- read.table(paste0(example,".profile.txt"),header=TRUE,sep="\t")
-as.matrix(xtabs(CPU.Time~Module,data=df))
+
+## Read total duration of profile
+cmd <- paste("amplxe-cl -report summary -result-dir ",example,".profile"," | grep Elapsed",sep="")
+out <- system(cmd,TRUE,FALSE,TRUE)
+endTime <- as.numeric(gsub("[^0-9.]*","",out))
+## Bisect to find relevant profile time interval
+f <- function(a,b){
+ time.filter <- paste(a,b,sep=":")
+ cmd <- paste("amplxe-cl -report top-down -result-dir ",example,".profile"," -time-filter ",time.filter,sep="")
+ system(cmd,TRUE,FALSE,TRUE)
+}
+g <- function(x,regexp,tol=0.05,direction=1){
+ empty <- length(grep(regexp,f(x[1],x[2])))==0
+ dx <- x[2]-x[1]
+ print(x)
+ if(dx<tol)return(x)
+ if(empty){
+ x <- x+dx*direction
+ } else {
+ if(direction>0)
+ x <- c(min(x),min(x)+dx/2)
+ else
+ x <- c(min(x)+dx/2,max(x))
+ }
+ return(g(x,regexp,tol,direction))
+}
+if(length(grep("CppAD",f(0,endTime)))!=0){
+ begin <- min(g(c(0,2),"CppAD"))
+ end <- max(g(c(endTime-1,endTime),"CppAD",direction=-1))
+ time.filter <- paste(begin,end,sep=":")
+ cat("Profiling time interval",time.filter,"\n")
+ cmd <- paste("amplxe-cl -report top-down -result-dir ",example,".profile"," -time-filter ",time.filter," > ",example,".profile.txt",sep="")
+ system(cmd)
+}
View
0 tmb_examples/sdv_multi.R 100755 → 100644
No changes.
View
0 tmb_examples/sdv_multi_compact.R 100755 → 100644
No changes.

0 comments on commit ce45200

Please sign in to comment.