added Thorson experiments

TMB_experiments/.gitignore

@@ -0,0 +1,3 @@
+*.o
+*.dll
+*.so

TMB_experiments/MakeADFun_windows_debug/MakeADFun_windows_debug.R

@@ -0,0 +1,54 @@
+
+MakeADFun_windows_debug = function( cpp_name, data, parameters, random=NULL, dir=getwd(), overwrite=FALSE, recompile=TRUE, ... ){
+
+  # Set local working directory
+  orig_dir = getwd()
+  setwd( dir )
+  on.exit( setwd(orig_dir) )
+
+  # Recompile
+  if( recompile==TRUE ){
+    unlink( dynlib(cpp_name) )
+    compile( paste0(cpp_name,".cpp"), "-O1 -g", DLLFLAGS="")
+    dyn.load( dynlib(cpp_name) )
+    message( "Recompiling .cpp file with appropriate compiler flags" )
+  }
+
+  # Get and save inputs
+  Other_inputs = list(...)
+  All_inputs = list( "data"=data, "parameters"=parameters, "random"=random, "Other_inputs"=Other_inputs )
+  save( All_inputs, file="All_inputs.RData")
+
+  # Write file to source
+  if( (paste0(cpp_name,".R") %in% list.files()) & overwrite==FALSE ){
+    message( "By default, we don't overwrite existing file ",cpp_name,".R")
+    Return = "Didn't attempt running"
+  }else{
+    sink( paste0(cpp_name,".R") )
+    cat("
+      library( TMB )
+      dyn.load( dynlib('",cpp_name,"') )
+      setwd('",dir,"')
+      load( 'All_inputs.RData' )
+      Obj = MakeADFun(data=All_inputs[['data']], parameters=All_inputs[['parameters']], random=All_inputs[['random']], All_inputs[['Other_inputs']])
+      save( Obj, file='Obj.RData')
+    ",fill=FALSE,sep="")
+    sink()
+
+    # Try running
+    Bdg_output = gdbsource(paste0(cpp_name,".R"))
+
+    # Sort out outcomes
+    if( length(grep("#0",Bdg_output))==0 ){
+      Return = MakeADFun(data=All_inputs[['data']], parameters=All_inputs[['parameters']], random=All_inputs[['random']], All_inputs[['Other_inputs']])
+      message("Compiled fine, and returning output from MakeADFun")
+    }else{
+      Return = Bdg_output
+      message("Did not compile, and returning output from bdbsource")
+      #print( Bdg_output )
+    }
+  }
+
+  # Return
+  return( Return )
+}

TMB_experiments/MakeADFun_windows_debug/test_MakeADFun_windows_debug.R

@@ -0,0 +1,31 @@
+
+File = "C:/Users/James.Thorson/Desktop/Project_git/TMB_experiments/windows_debugger"
+setwd( File )
+source( "MakeADFun_windows_debug.R" )
+
+#####################
+# Simulate data
+######################
+
+Factor = rep( 1:10, each=10)
+Z = rnorm( length(unique(Factor)), mean=0, sd=1)
+
+X0 = 0
+
+Y = Z[Factor] + X0 + rnorm( length(Factor), mean=0, sd=1)
+
+######################
+# Experiment with debugger
+######################
+
+library(TMB)
+cpp_name = "problem"
+#compile( paste0(cpp_name,".cpp"), "-O1 -g", DLLFLAGS="")
+compile( paste0(cpp_name,".cpp") )
+
+data = list( "n_data"=length(Y), "n_factors"=length(unique(Factor)), "Factor"=Factor-1, "Y"=Y)
+parameters = list( "X0"=-10, "log_SD0"=2, "log_SDZ"=2, "Z"=rep(0,data$n_factor) )
+random = c("Z")
+
+#dyn.load( dynlib(cpp_name) )
+Output = MakeADFun_windows_debug( data=data, cpp_name=cpp_name, dir=File, parameters=parameters, random=random, overwrite=TRUE )

TMB_experiments/README.md

@@ -0,0 +1,2 @@
+# TMB_experiments
+Experiments involving TMB

TMB_experiments/experiment_with_SCALE/experiment_involving_SCALE.R

@@ -0,0 +1,54 @@
+
+library(TMB)
+
+mu_b = 2
+sigma_b = 2
+sigma_c = 1
+
+b_g = rnorm(10, mean=mu_b, sd=sigma_b)
+g_i = rep(1:10,each=10)
+y_i = b_g[g_i] + rnorm(length(g_i), mean=0, sd=sigma_c)
+
+# sanity check
+library(lme4)
+Lm = lmer( y_i ~ 1 + 1|g_i )
+summary(Lm)
+
+# compile both models
+compile( "linear_model.cpp" )
+compile( "linear_model_using_scale.cpp" )
+
+# Build inputs (same for both models)
+Data = list( "g_i"=g_i-1, "y_i"=y_i)
+Parameters = list( "mu"=-10, "logsigma_y"=2, "logsigma_g"=2, "b_g"=rep(0,length(unique(g_i))) )
+Random = c("b_g")
+
+
+################
+# In TMB, without using scale
+# WORKS CORRECTLY
+################
+
+# Build object
+dyn.load( dynlib("linear_model") )
+Obj = MakeADFun(data=Data, parameters=Parameters, random=Random)
+
+# Optimize
+Opt = nlminb( start=Obj$par, objective=Obj$fn, gradient=Obj$gr, control=list("trace"=1) )
+Report = Obj$report()
+exp( unlist(Report[c('logsigma_g','logsigma_y')]) )
+
+################
+# In TMB, using scale
+# DOES NOT WORK CORRECTLY
+################
+
+# Build object
+dyn.load( dynlib("linear_model_using_scale") )
+Obj = MakeADFun(data=Data, parameters=Parameters, random=Random)
+
+# Optimize
+Opt = nlminb( start=Obj$par, objective=Obj$fn, gradient=Obj$gr, control=list("trace"=1) )
+Report = Obj$report()
+exp( unlist(Report[c('logsigma_g','logsigma_y')]) )
+

TMB_experiments/experiment_with_SCALE/linear_model.cpp

@@ -0,0 +1,40 @@
+// Space time
+#include <TMB.hpp>
+template<class Type>
+Type objective_function<Type>::operator() ()
+{
+  using namespace density;
+
+  // Data
+  DATA_VECTOR( y_i );
+  DATA_FACTOR( g_i );
+
+  // Parameters
+  PARAMETER( mu );
+  PARAMETER( logsigma_y );
+  PARAMETER( logsigma_g );
+  PARAMETER_VECTOR( b_g );
+
+  int n_data = y_i.size();
+  int n_factors = b_g.size();
+
+  // Objective funcction
+  Type jnll = 0;
+
+  // Probability of random coefficients
+  matrix<Type> Cov_gg(n_factors, n_factors);
+  Cov_gg.setIdentity();
+  Cov_gg = exp(2*logsigma_g) * Cov_gg;
+  MVNORM_t<Type> nll_mvnorm(Cov_gg);
+  jnll += nll_mvnorm( b_g );
+
+  // Probability of data conditional on fixed and random effect values
+  for( int i=0; i<n_data; i++){
+    jnll -= dnorm( y_i(i), mu + b_g(g_i(i)), exp(logsigma_y), true );
+  }
+
+  // Reporting
+  REPORT( logsigma_y );
+  REPORT( logsigma_g );
+  return jnll;
+}

TMB_experiments/experiment_with_SCALE/linear_model_using_scale.cpp

@@ -0,0 +1,43 @@
+// Space time
+#include <TMB.hpp>
+template<class Type>
+Type objective_function<Type>::operator() ()
+{
+  using namespace density;
+
+  // Data
+  DATA_VECTOR( y_i );
+  DATA_FACTOR( g_i );
+
+  // Parameters
+  PARAMETER( mu );
+  PARAMETER( logsigma_y );
+  PARAMETER( logsigma_g );
+  PARAMETER_VECTOR( b_g );
+
+  int n_data = y_i.size();
+  int n_factors = b_g.size();
+
+  // Objective funcction
+  Type jnll = 0;
+
+  // Rescale
+  vector<Type> beta_g(n_factors);
+  beta_g = b_g * exp(logsigma_g);
+
+  // Probability of random coefficients
+  matrix<Type> Cov_gg(n_factors, n_factors);
+  Cov_gg.setIdentity();
+  MVNORM_t<Type> nll_mvnorm(Cov_gg);
+  jnll += SCALE(nll_mvnorm, exp(logsigma_g))( beta_g );
+
+  // Probability of data conditional on fixed and random effect values
+  for( int i=0; i<n_data; i++){
+    jnll -= dnorm( y_i(i), mu + beta_g(g_i(i)), exp(logsigma_y), true );
+  }
+
+  // Reporting
+  REPORT( logsigma_y );
+  REPORT( logsigma_g );
+  return jnll;
+}

TMB_experiments/experiment_with_oneStepPredict/linear_mixed_model.R

@@ -0,0 +1,45 @@
+
+setwd( "C:/Users/James.Thorson/Desktop/Project_git/TMB_experiments/experiment_with_oneStepPredict" )
+
+######################
+# Simulate data
+######################
+
+Factor = rep( 1:10, each=10)
+Z = rnorm( length(unique(Factor)), mean=0, sd=1)
+
+X0 = 0
+
+Y = Z[Factor] + X0 + rnorm( length(Factor), mean=0, sd=1)
+
+######################
+# Run in TMB
+######################
+
+library(TMB)
+
+Version = c( "linear_mixed_model" )[1]
+compile( paste0(Version,".cpp") )
+
+Data = list( "n_data"=length(Y), "n_factors"=length(unique(Factor)), "Factor"=Factor-1, "Y"=Y)
+Parameters = list( "X0"=-10, "log_SD0"=2, "log_SDZ"=2, "Z"=rep(0,Data$n_factor) )
+Random = c("Z")
+
+dyn.load( dynlib("linear_mixed_model") )
+Obj = MakeADFun(data=Data, parameters=Parameters, random=Random)
+
+# Prove that function and gradient calls work
+Obj$fn( Obj$par )
+Obj$gr( Obj$par )
+
+# Optimize
+start_time = Sys.time()
+Opt = nlminb( start=Obj$par, objective=Obj$fn, gradient=Obj$gr, control=list("trace"=1) )
+  Opt[["final_gradient"]] = Obj$gr( Opt$par )
+  Opt[["total_time"]] = Sys.time() - start_time
+  Report = Obj$report()
+  SD = sdreport( Obj, bias.correct=TRUE )
+
+# One-step-predict
+oneStepPredict( Obj, observation.name="Y", method="fullGaussian" )
+oneStepPredict( Obj, observation.name="Y", method="oneStepGeneric", data.term.indicator="keep" )

TMB_experiments/experiment_with_oneStepPredict/linear_mixed_model.cpp

@@ -0,0 +1,55 @@
+// Space time
+#include <TMB.hpp>
+template<class Type>
+Type objective_function<Type>::operator() ()
+{
+  // Data
+  DATA_INTEGER( n_data );
+  DATA_INTEGER( n_factors );
+  DATA_FACTOR( Factor );
+  DATA_VECTOR( Y );                
+ DATA_VECTOR_INDICATOR(keep, Y);
+
+  // Parameters
+  PARAMETER( X0 );
+  PARAMETER( log_SD0 );
+  PARAMETER( log_SDZ );
+  PARAMETER_VECTOR( Z );
+
+  // Objective funcction
+  Type jnll = 0;
+
+  // Probability of data conditional on fixed and random effect values
+  for( int i=0; i<n_data; i++){
+    jnll -= dnorm( Y(i), X0 + Z(Factor(i)), exp(log_SD0), true );
+  }
+
+  // Probability of random coefficients
+  for( int i=0; i<n_factors; i++){
+    jnll -= dnorm( Z(i), Type(0.0), exp(log_SDZ), true );
+  }
+
+  // Reporting
+  Type SDZ = exp(log_SDZ);
+  Type SD0 = exp(log_SD0);
+  ADREPORT( SDZ );
+  REPORT( SDZ );
+  ADREPORT( SD0 );
+  REPORT( SD0 );
+  ADREPORT( Z );
+  REPORT( Z );
+  ADREPORT( X0 );
+  REPORT( X0 );
+
+  // bias-correction testing
+  Type MeanZ = Z.sum() / Z.size();
+  Type SampleVarZ = ( (Z-MeanZ) * (Z-MeanZ) ).sum();
+  Type SampleSDZ = pow( SampleVarZ + 1e-20, 0.5);
+  REPORT( SampleVarZ );
+  REPORT( SampleSDZ );  
+  ADREPORT( SampleVarZ );
+  ADREPORT( SampleSDZ );  
+
+
+  return jnll;
+}

TMB_experiments/experiment_with_parallel_accumulator/linreg_parallel.R

@@ -0,0 +1,39 @@
+
+setwd("C:/Users/James.Thorson/Desktop/Project_git/TMB_experiments/experiment_with_parallel_accumulator")
+
+library(TMB)
+set.seed(123)
+x <- seq(0,10,length=50001)
+data <- list(Y=rnorm(length(x))+x,x=x)
+parameters <- list(a=0,b=0,logSigma=0)
+compile( "linreg_parallel.cpp" )
+dyn.load(dynlib("linreg_parallel"))
+obj <- MakeADFun(data,parameters,DLL="linreg_parallel")
+obj$hessian <- TRUE
+opt <- do.call("optim",obj)
+
+#########################
+# Experiment with parallelization
+#########################
+
+# Parallelized experimentations
+#ben <- benchmark(obj, n=1e3, cores=1:4)
+#plot(ben)
+
+# Parallel
+ben <- benchmark(obj, n=1000, cores=1:4, expr=expression(do.call("optim",obj)))
+png( file="Benchmark.png", width=6, height=6, res=200, units="in")
+  plot(ben)
+dev.off()
+
+##########################
+# Help file version
+##########################
+
+runExample("linreg_parallel",thisR=TRUE) ## Create obj
+ben <- benchmark(obj,n=100,cores=1:4)
+plot(ben)
+ben <- benchmark(obj,n=10,cores=1:4,expr=expression(do.call("optim",obj)))
+plot(ben)
+
+

TMB_experiments/experiment_with_parallel_accumulator/linreg_parallel.cpp

@@ -0,0 +1,13 @@
+#include <TMB.hpp>
+template<class Type>
+Type objective_function<Type>::operator() ()
+{
+  DATA_VECTOR(Y);
+  DATA_VECTOR(x);
+  PARAMETER(a);
+  PARAMETER(b);
+  PARAMETER(logSigma);
+  parallel_accumulator<Type> nll(this);
+  for(int i=0;i<x.size();i++)nll-=dnorm(Y[i],a+b*x[i],exp(logSigma),true);
+  return nll;
+}