TMB Documentation  v1.9.11
tmb_core.hpp
Go to the documentation of this file.
1 // Copyright (C) 2013-2015 Kasper Kristensen
2 // License: GPL-2
3 
8 /*
9  Call to external C++ code can potentially result in exeptions that
10  will crash R. However, we do not want R to crash on failed memory
11  allocations. Therefore:
12 
13  * All interface functions (those called with .Call from R) must have
14  TMB_TRY wrapped around CppAD/Eigen code that allocates memory.
15 
16  * Special attention must be payed to parallel code, as each thread
17  is responsible for catching its own exceptions.
18 */
19 
20 #ifndef TMB_TRY
21 #define TMB_TRY try
22 #endif
23 // By default we only accept 'bad_alloc' as a valid exception. Everything else => debugger !
24 // Behaviour can be changed by re-defining this macro.
25 #ifndef TMB_CATCH
26 #define TMB_CATCH catch(std::bad_alloc& excpt)
27 #endif
28 // Inside the TMB_CATCH comes 'cleanup code' followed by this error
29 // call (allowed to depend on the exception 'excpt')
30 // Error message can be changed by re-defining this macro.
31 #ifndef TMB_ERROR_BAD_ALLOC
32 #define TMB_ERROR_BAD_ALLOC \
33 Rf_error("Caught exception '%s' in function '%s'\n", \
34  excpt.what(), \
35  __FUNCTION__)
36 #endif
37 // Error call comes outside TMB_CATCH in OpenMP case (so *cannot*
38 // depend on exception e.g. 'excpt')
39 // Error message can be changed by re-defining this macro.
40 #ifndef TMB_ERROR_BAD_THREAD_ALLOC
41 #define TMB_ERROR_BAD_THREAD_ALLOC \
42 Rf_error("Caught exception '%s' in function '%s'\n", \
43  bad_thread_alloc, \
44  __FUNCTION__)
45 #endif
46 
47 /* Memory manager:
48  Count the number of external pointers alive.
49  When total number is zero it is safe to dyn.unload
50  the library.
51 */
52 #include <set>
53 extern "C" void finalizeDoubleFun(SEXP x);
54 extern "C" void finalizeADFun(SEXP x);
55 extern "C" void finalizeparallelADFun(SEXP x);
56 extern "C" SEXP FreeADFunObject(SEXP f) CSKIP ({
57  SEXP tag = R_ExternalPtrTag(f);
58  if (tag == Rf_install("DoubleFun")) {
59  finalizeDoubleFun(f);
60  }
61  else if (tag == Rf_install("ADFun")) {
62  finalizeADFun(f);
63  }
64  else if (tag == Rf_install("parallelADFun")) {
65  finalizeparallelADFun(f);
66  }
67  else {
68  Rf_error("Unknown external ptr type");
69  }
70  R_ClearExternalPtr(f); // Set pointer to 'nil'
71  return R_NilValue;
72 })
74 struct memory_manager_struct {
75  int counter;
77  std::set<SEXP> alive;
79  void RegisterCFinalizer(SEXP list);
81  void CallCFinalizer(SEXP x);
83  void clear();
84  memory_manager_struct();
85 };
86 #ifndef WITH_LIBTMB
87 void memory_manager_struct::RegisterCFinalizer(SEXP x) {
88  counter++;
89  alive.insert(x);
90 }
91 void memory_manager_struct::CallCFinalizer(SEXP x){
92  counter--;
93  alive.erase(x);
94 }
95 void memory_manager_struct::clear(){
96  std::set<SEXP>::iterator it;
97  while (alive.size() > 0) {
98  FreeADFunObject(*alive.begin());
99  }
100 }
101 memory_manager_struct::memory_manager_struct(){
102  counter=0;
103 }
104 #endif
105 TMB_EXTERN memory_manager_struct memory_manager;
106 
117 #ifdef WITH_LIBTMB
118 SEXP ptrList(SEXP x);
119 #else
120 SEXP ptrList(SEXP x)
121 {
122  SEXP ans,names;
123  PROTECT(ans=Rf_allocVector(VECSXP,1));
124  PROTECT(names=Rf_allocVector(STRSXP,1));
125  SET_VECTOR_ELT(ans,0,x);
126  SET_STRING_ELT(names,0,Rf_mkChar("ptr"));
127  Rf_setAttrib(ans,R_NamesSymbol,names);
128  memory_manager.RegisterCFinalizer(x);
129  UNPROTECT(2);
130  return ans;
131 }
132 #endif
133 
134 extern "C"{
135 #ifdef LIB_UNLOAD
136 #include <R_ext/Rdynload.h>
137  void LIB_UNLOAD(DllInfo *dll)
138  {
139  if(memory_manager.counter>0)Rprintf("Warning: %d external pointers will be removed\n",memory_manager.counter);
140  memory_manager.clear();
141  for(int i=0;i<1000;i++){ // 122 seems to be sufficient.
142  if(memory_manager.counter>0){
143  R_gc();
144  R_RunExitFinalizers();
145  } else break;
146  }
147  if(memory_manager.counter>0)Rf_error("Failed to clean. Please manually clean up before unloading\n");
148  }
149 #endif
150 }
151 
152 #ifdef _OPENMP
153 TMB_EXTERN bool _openmp CSKIP( =true; )
154 #else
155 TMB_EXTERN bool _openmp CSKIP( =false; )
156 #endif
157 
159 template<class ADFunPointer>
160 void optimizeTape(ADFunPointer pf){
161  if(!config.optimize.instantly){
162  /* Drop out */
163  return;
164  }
165  if (!config.optimize.parallel){
166 #ifdef _OPENMP
167 #pragma omp critical
168 #endif
169  { /* Avoid multiple tape optimizations at the same time (to reduce memory) */
170  if(config.trace.optimize)std::cout << "Optimizing tape... ";
171  pf->optimize();
172  if(config.trace.optimize)std::cout << "Done\n";
173  }
174  }
175  else
176  { /* Allow multiple tape optimizations at the same time */
177  if(config.trace.optimize)std::cout << "Optimizing tape... ";
178  pf->optimize();
179  if(config.trace.optimize)std::cout << "Done\n";
180  }
181 }
182 
183 /* Macros to obtain data and parameters from R */
184 
208 #define TMB_OBJECTIVE_PTR \
209 this
210 
213 #define PARAMETER_MATRIX(name) \
214 tmbutils::matrix<Type> name(TMB_OBJECTIVE_PTR -> fillShape( \
215 asMatrix<Type> ( TMB_OBJECTIVE_PTR -> getShape( #name, &Rf_isMatrix) ), \
216 #name) );
217 
220 #define PARAMETER_VECTOR(name) \
221 vector<Type> name(TMB_OBJECTIVE_PTR -> fillShape( \
222 asVector<Type>(TMB_OBJECTIVE_PTR -> getShape(#name, &Rf_isReal )), \
223 #name));
224 
227 #define PARAMETER(name) \
228 Type name(TMB_OBJECTIVE_PTR -> fillShape( \
229 asVector<Type>(TMB_OBJECTIVE_PTR -> getShape(#name,&isNumericScalar)), \
230 #name)[0]);
231 
236 #define DATA_VECTOR(name) \
237 vector<Type> name; \
238 if (!Rf_isNull(getListElement(TMB_OBJECTIVE_PTR -> parameters,#name))){ \
239  name = TMB_OBJECTIVE_PTR -> fillShape(asVector<Type>( \
240  TMB_OBJECTIVE_PTR -> getShape(#name, &Rf_isReal )), #name); \
241 } else { \
242  name = asVector<Type>(getListElement( \
243  TMB_OBJECTIVE_PTR -> data,#name,&Rf_isReal )); \
244 }
245 
248 #define DATA_MATRIX(name) \
249 matrix<Type> name(asMatrix<Type>( \
250 getListElement(TMB_OBJECTIVE_PTR -> data, #name, &Rf_isMatrix)));
251 
254 #define DATA_SCALAR(name) \
255 Type name(asVector<Type>(getListElement(TMB_OBJECTIVE_PTR -> data, \
256 #name,&isNumericScalar))[0]);
257 
260 #define DATA_INTEGER(name) int name(CppAD::Integer(asVector<Type>( \
261 getListElement(TMB_OBJECTIVE_PTR -> data, \
262 #name, &isNumericScalar))[0]));
263 
281 #define DATA_FACTOR(name) vector<int> name(asVector<int>( \
282 getListElement(TMB_OBJECTIVE_PTR -> data, #name, &Rf_isReal )));
283 
287 #define DATA_IVECTOR(name) vector<int> name(asVector<int>( \
288 getListElement(TMB_OBJECTIVE_PTR -> data, #name, &Rf_isReal )));
289 
292 #define NLEVELS(name) \
293 LENGTH(Rf_getAttrib(getListElement(TMB_OBJECTIVE_PTR -> data, #name), \
294 Rf_install("levels")))
295 
299 #define DATA_SPARSE_MATRIX(name) \
300 Eigen::SparseMatrix<Type> name(tmbutils::asSparseMatrix<Type>( \
301 getListElement(TMB_OBJECTIVE_PTR -> data, \
302 #name, &isValidSparseMatrix)));
303 
304 // NOTE: REPORT() constructs new SEXP so never report in parallel!
313 #define REPORT(name) \
314 if( isDouble<Type>::value && \
315  TMB_OBJECTIVE_PTR -> current_parallel_region<0 ) \
316 { \
317  SEXP _TMB_temporary_sexp_; \
318  PROTECT( _TMB_temporary_sexp_ = asSEXP(name) ); \
319  Rf_defineVar(Rf_install(#name), \
320  _TMB_temporary_sexp_, TMB_OBJECTIVE_PTR -> report); \
321  UNPROTECT(1); \
322 }
323 
329 #define SIMULATE \
330 if(isDouble<Type>::value && TMB_OBJECTIVE_PTR -> do_simulate)
331 
342 #define ADREPORT(name) \
343 TMB_OBJECTIVE_PTR -> reportvector.push(name, #name);
344 
345 #define PARALLEL_REGION \
346 if( TMB_OBJECTIVE_PTR -> parallel_region() )
347 
352 #define DATA_ARRAY(name) \
353 tmbutils::array<Type> name; \
354 if (!Rf_isNull(getListElement(TMB_OBJECTIVE_PTR -> parameters,#name))){ \
355  name = TMB_OBJECTIVE_PTR -> fillShape(tmbutils::asArray<Type>( \
356  TMB_OBJECTIVE_PTR -> getShape(#name, &Rf_isArray)), #name); \
357 } else { \
358  name = tmbutils::asArray<Type>(getListElement( \
359  TMB_OBJECTIVE_PTR -> data, #name, &Rf_isArray)); \
360 }
361 
364 #define PARAMETER_ARRAY(name) \
365 tmbutils::array<Type> name(TMB_OBJECTIVE_PTR -> fillShape( \
366 tmbutils::asArray<Type>(TMB_OBJECTIVE_PTR -> getShape( \
367 #name, &Rf_isArray)), #name));
368 
371 #define DATA_IMATRIX(name) \
372 matrix<int> name(asMatrix<int>( \
373 getListElement(TMB_OBJECTIVE_PTR -> data,#name, &Rf_isMatrix)));
374 
377 #define DATA_IARRAY(name) \
378 tmbutils::array<int> name(tmbutils::asArray<int>( \
379 getListElement(TMB_OBJECTIVE_PTR -> data, #name, &Rf_isArray)));
380 
394 #define DATA_STRING(name) \
395 std::string name = \
396  CHAR(STRING_ELT(getListElement(TMB_OBJECTIVE_PTR -> data, #name), 0));
397 
433 #define DATA_STRUCT(name, struct) \
434 struct<Type> name(getListElement(TMB_OBJECTIVE_PTR -> data, #name));
435 
440 template<class VT, class Type = typename VT::Scalar>
441 struct data_indicator : VT{
443  VT cdf_lower;
445  VT cdf_upper;
447  vector<int> ord;
449  bool osa_flag;
451  data_indicator() { osa_flag = false; }
456  data_indicator(VT obs, bool init_one = false){
457  VT::operator=(obs);
458  if (init_one) VT::fill(Type(1.0));
459  cdf_lower = obs; cdf_lower.setZero();
460  cdf_upper = obs; cdf_upper.setZero();
461  osa_flag = false;
462  }
464  void fill(vector<Type> p, SEXP ord_){
465  int n = (*this).size();
466  if(p.size() >= n ) VT::operator=(p.segment(0, n));
467  if(p.size() >= 2*n) cdf_lower = p.segment(n, n);
468  if(p.size() >= 3*n) cdf_upper = p.segment(2 * n, n);
469  if(!Rf_isNull(ord_)) {
470  this->ord = asVector<int>(ord_);
471  }
472  for (int i=0; i<p.size(); i++) {
473  osa_flag |= CppAD::Variable(p[i]);
474  }
475  }
478  data_indicator segment(int pos, int n) {
479  data_indicator ans ( VT::segment(pos, n) );
480  ans.cdf_lower = cdf_lower.segment(pos, n);
481  ans.cdf_upper = cdf_upper.segment(pos, n);
482  if (ord.size() != 0) {
483  ans.ord = ord.segment(pos, n);
484  }
485  ans.osa_flag = osa_flag;
486  return ans;
487  }
489  vector<int> order() {
490  int n = this->size();
491  vector<int> ans(n);
492  if (ord.size() == 0) {
493  for (int i=0; i<n; i++)
494  ans(i) = i;
495  } else {
496  if (ord.size() != n) Rf_error("Unexpected 'ord.size() != n'");
497  std::vector<std::pair<int, int> > y(n);
498  for (int i=0; i<n; i++) {
499  y[i].first = ord[i];
500  y[i].second = i;
501  }
502  std::sort(y.begin(), y.end()); // sort inplace
503  for (int i=0; i<n; i++) {
504  ans[i] = y[i].second;
505  }
506  }
507  return ans;
508  }
510  bool osa_active() { return osa_flag; }
511 };
512 
518 #define DATA_ARRAY_INDICATOR(name, obs) \
519 data_indicator<tmbutils::array<Type> > name(obs, true); \
520 if (!Rf_isNull(getListElement(TMB_OBJECTIVE_PTR -> parameters,#name))){ \
521  name.fill( TMB_OBJECTIVE_PTR -> fillShape(asVector<Type>( \
522  TMB_OBJECTIVE_PTR -> getShape(#name, &Rf_isReal )), \
523  #name), \
524  Rf_getAttrib( \
525  TMB_OBJECTIVE_PTR -> getShape(#name, &Rf_isReal ), \
526  Rf_install("ord")) ); \
527 }
528 
534 #define DATA_VECTOR_INDICATOR(name, obs) \
535 data_indicator<tmbutils::vector<Type> > name(obs, true); \
536 if (!Rf_isNull(getListElement(TMB_OBJECTIVE_PTR -> parameters,#name))){ \
537  name.fill( TMB_OBJECTIVE_PTR -> fillShape(asVector<Type>( \
538  TMB_OBJECTIVE_PTR -> getShape(#name, &Rf_isReal )), \
539  #name), \
540  Rf_getAttrib( \
541  TMB_OBJECTIVE_PTR -> getShape(#name, &Rf_isReal ), \
542  Rf_install("ord")) ); \
543 }
544 
545 // kasper: Not sure used anywhere
549 template<class Type>
550 matrix<int> HessianSparsityPattern(ADFun<Type> *pf){
551  int n=pf->Domain();
552  vector<bool> Px(n * n);
553  for(int i = 0; i < n; i++)
554  {
555  for(int j = 0; j < n; j++)
556  Px[ i * n + j ] = false;
557  Px[ i * n + i ] = true;
558  }
559  pf->ForSparseJac(n, Px);
560  vector<bool> Py(1); Py[0]=true;
561  vector<int> tmp = (pf->RevSparseHes(n,Py)).template cast<int>();
562  return asMatrix(tmp, n, n);
563 }
564 
566 void Independent(vector<double> x)CSKIP({})
567 
569 template <class Type>
570 struct report_stack{
571  std::vector<const char*> names;
572  std::vector<vector<int> > namedim;
573  std::vector<Type> result;
574  void clear(){
575  names.resize(0);
576  namedim.resize(0);
577  result.resize(0);
578  }
579  // Get dimension of various object types
580  vector<int> getDim(const matrix<Type> &x) {
581  vector<int> dim(2);
582  dim << x.rows(), x.cols();
583  return dim;
584  }
585  vector<int> getDim(const tmbutils::array<Type> &x) {
586  return x.dim;
587  }
588  template<class Other> // i.e. vector or expression
589  vector<int> getDim(const Other &x) {
590  vector<int> dim(1);
591  dim << x.size();
592  return dim;
593  }
594  // push vector, matrix or array
595  template<class Vector_Matrix_Or_Array>
596  void push(Vector_Matrix_Or_Array x, const char* name) {
597  names.push_back(name);
598  namedim.push_back(getDim(x));
599  Eigen::Array<Type, Eigen::Dynamic, Eigen::Dynamic> xa(x);
600  result.insert(result.end(), xa.data(), xa.data() + x.size());
601  }
602  // push scalar (convert to vector case)
603  void push(Type x, const char* name){
604  vector<Type> xvec(1);
605  xvec[0] = x;
606  push(xvec, name);
607  }
608  // Eval: cast to vector<Type>
609  vector<Type> operator()() {
610  return result;
611  }
612  /* Get names (with replicates) to R */
613  SEXP reportnames()
614  {
615  int n = result.size();
616  SEXP nam;
617  PROTECT( nam = Rf_allocVector(STRSXP, n) );
618  int k = 0;
619  for(size_t i = 0; i < names.size(); i++) {
620  int namelength = namedim[i].prod();
621  for(int j = 0; j < namelength; j++) {
622  SET_STRING_ELT(nam, k, Rf_mkChar(names[i]) );
623  k++;
624  }
625  }
626  UNPROTECT(1);
627  return nam;
628  }
629  /* Get AD reported object dims */
630  SEXP reportdims() {
631  SEXP ans, nam;
632  typedef vector<vector<int> > VVI;
633  PROTECT( ans = asSEXP(VVI(namedim)) );
634  PROTECT( nam = Rf_allocVector(STRSXP, names.size()) );
635  for(size_t i = 0; i < names.size(); i++) {
636  SET_STRING_ELT(nam, i, Rf_mkChar(names[i]));
637  }
638  Rf_setAttrib(ans, R_NamesSymbol, nam);
639  UNPROTECT(2);
640  return ans;
641  }
642  EIGEN_DEFAULT_DENSE_INDEX_TYPE size(){return result.size();}
643 }; // report_stack
644 
645 extern "C" {
646  void GetRNGstate(void);
647  void PutRNGstate(void);
648 }
649 
651 template <class Type>
652 class objective_function
653 {
654 // private:
655 public:
656  SEXP data;
657  SEXP parameters;
658  SEXP report;
659 
660  int index;
661  vector<Type> theta;
662  vector<const char*> thetanames;
663  report_stack<Type> reportvector;
664  bool reversefill; // used to find the parameter order in user template (not anymore - use pushParname instead)
665  vector<const char*> parnames;
668  void pushParname(const char* x){
669  parnames.conservativeResize(parnames.size()+1);
670  parnames[parnames.size()-1]=x;
671  }
672 
673  /* ================== For parallel Hessian computation
674  Need three different parallel evaluation modes:
675  (1) *Parallel mode* where a parallel region is evaluated iff
676  current_parallel_region == selected_parallel_region
677  (2) *Serial mode* where all parallel region tests are evaluated
678  to TRUE so that "PARALLEL_REGION" tests are effectively removed.
679  A negative value of "current_parallel_region" or "selected_parallel_region"
680  is used to select this mode (the default).
681  (3) *Count region mode* where statements inside "PARALLEL_REGION{...}"
682  are *ignored* and "current_parallel_region" is increased by one each
683  time a parallel region is visited.
684  NOTE: The macro "PARALLEL_REGION" is supposed to be defined as
685  #define PARALLEL_REGION if(this->parallel_region())
686  where the function "parallel_region" does the book keeping.
687  */
688  bool parallel_ignore_statements;
689  int current_parallel_region; /* Identifier of a code-fragment of user template */
690  int selected_parallel_region; /* Consider _this_ code-fragment */
691  int max_parallel_regions; /* Max number of parallel region identifiers,
692  e.g. max_parallel_regions=config.nthreads;
693  probably best in most cases. */
694  bool parallel_region(){ /* Is this the selected parallel region ? */
695  bool ans;
696  if(config.autopar || current_parallel_region<0 || selected_parallel_region<0)return true; /* Serial mode */
697  ans = (selected_parallel_region==current_parallel_region) && (!parallel_ignore_statements);
698  current_parallel_region++;
699  if(max_parallel_regions>0)current_parallel_region=current_parallel_region % max_parallel_regions;
700  return ans;
701  }
702  /* Note: Some other functions rely on "count_parallel_regions" to run through the users code (!) */
703  int count_parallel_regions(){
704  current_parallel_region=0; /* reset counter */
705  selected_parallel_region=0;
706  parallel_ignore_statements=true; /* Do not evaluate stuff inside PARALLEL_REGION{...} */
707  this->operator()(); /* Run through users code */
708  if (config.autopar) return 0;
709  if(max_parallel_regions>0)return max_parallel_regions;
710  else
711  return current_parallel_region;
712  }
713  void set_parallel_region(int i){ /* Select parallel region (from within openmp loop) */
714  current_parallel_region=0;
715  selected_parallel_region=i;
716  parallel_ignore_statements=false;
717  }
718 
719  bool do_simulate;
720  void set_simulate(bool do_simulate_) {
721  do_simulate = do_simulate_;
722  }
723 
724  /* data_ and parameters_ are R-lists containing R-vectors or R-matrices.
725  report_ is an R-environment.
726  The elements of the vector "unlist(parameters_)" are filled into "theta"
727  which contains the default parameter-values. This happens during the
728  *construction* of the objective_function object.
729  The user defined template "objective_function::operator()" is called
730  from "MakeADFunObject" which tapes the operations and creates the final
731  ADFun-object.
732  */
734  objective_function(SEXP data, SEXP parameters, SEXP report) :
735  data(data), parameters(parameters), report(report), index(0)
736  {
737  /* Fill theta with the default parameters.
738  Pass R-matrices column major. */
739  theta.resize(nparms(parameters));
740  int length_parlist = Rf_length(parameters);
741  for(int i = 0, counter = 0; i < length_parlist; i++) {
742  // x = parameters[[i]]
743  SEXP x = VECTOR_ELT(parameters, i);
744  int nx = Rf_length(x);
745  double* px = REAL(x);
746  for(int j = 0; j < nx; j++) {
747  theta[counter++] = Type( px[j] );
748  }
749  }
750  thetanames.resize(theta.size());
751  for(int i=0;i<thetanames.size();i++)thetanames[i]="";
752  current_parallel_region=-1;
753  selected_parallel_region=-1;
754  max_parallel_regions=-1;
755 #ifdef _OPENMP
756  max_parallel_regions = config.nthreads;
757 #endif
758  reversefill=false;
759  do_simulate = false;
760  GetRNGstate(); /* Read random seed from R. Note: by default we do
761  not write the seed back to R *after*
762  simulation. This ensures that multiple tapes for
763  one model object get the same seed. When in
764  simulation mode (enabled when calling
765  obj$simulate() from R) we *do* write the seed
766  back after simulation in order to get varying
767  replicates. */
768  }
769 
771  void sync_data() {
772  SEXP env = ENCLOS(this->report);
773  this->data = Rf_findVar(Rf_install("data"), env);
774  }
775 
777  SEXP defaultpar()
778  {
779  int n=theta.size();
780  SEXP res;
781  SEXP nam;
782  PROTECT(res=Rf_allocVector(REALSXP,n));
783  PROTECT(nam=Rf_allocVector(STRSXP,n));
784  for(int i=0;i<n;i++){
785  //REAL(res)[i]=CppAD::Value(theta[i]);
786  REAL(res)[i]=value(theta[i]);
787  SET_STRING_ELT(nam,i,Rf_mkChar(thetanames[i]));
788  }
789  Rf_setAttrib(res,R_NamesSymbol,nam);
790  UNPROTECT(2);
791  return res;
792  }
793 
795  SEXP parNames()
796  {
797  int n=parnames.size();
798  SEXP nam;
799  PROTECT(nam=Rf_allocVector(STRSXP,n));
800  for(int i=0;i<n;i++){
801  SET_STRING_ELT(nam,i,Rf_mkChar(parnames[i]));
802  }
803  UNPROTECT(1);
804  return nam;
805  }
806 
807  /* FIXME: "Value" should be "var2par" I guess
808  kasper: Why not use asDouble defined previously? */
818  double value(double x){return x;}
819  double value(AD<double> x){return CppAD::Value(x);}
820  double value(AD<AD<double> > x){return CppAD::Value(CppAD::Value(x));}
821  double value(AD<AD<AD<double> > > x){return CppAD::Value(CppAD::Value(CppAD::Value(x)));}
822 #ifdef TMBAD_FRAMEWORK
823  double value(TMBad::ad_aug x){return x.Value();}
824 #endif
825 
828  int nparms(SEXP obj)
829  {
830  int count=0;
831  for(int i=0;i<Rf_length(obj);i++){
832  if(!Rf_isReal(VECTOR_ELT(obj,i)))Rf_error("PARAMETER COMPONENT NOT A VECTOR!");
833  count+=Rf_length(VECTOR_ELT(obj,i));
834  }
835  return count;
836  }
837 
838  /* The "fill functions" are all used to populate parameter vectors,
839  arrays, matrices etc with the values of the parameter vector theta. */
840  void fill(vector<Type> &x, const char *nam)
841  {
842  pushParname(nam);
843  for(int i=0;i<x.size();i++){
844  thetanames[index]=nam;
845  if(reversefill)theta[index++]=x[i];else x[i]=theta[index++];
846  }
847  }
848  void fill(matrix<Type> &x, const char *nam)
849  {
850  pushParname(nam);
851  for(int j=0;j<x.cols();j++){
852  for(int i=0;i<x.rows();i++){
853  thetanames[index]=nam;
854  if(reversefill)theta[index++]=x(i,j);else x(i,j)=theta[index++];
855  }
856  }
857  }
858  template<class ArrayType>
859  void fill(ArrayType &x, const char *nam)
860  {
861  pushParname(nam);
862  for(int i=0;i<x.size();i++){
863  thetanames[index]=nam;
864  if(reversefill)theta[index++]=x[i];else x[i]=theta[index++];
865  }
866  }
867 
868  /* Experiment: new map feature - currently arrays only */
869  template<class ArrayType>
870  void fillmap(ArrayType &x, const char *nam)
871  {
872  pushParname(nam);
873  SEXP elm=getListElement(parameters,nam);
874  int* map=INTEGER(Rf_getAttrib(elm,Rf_install("map")));
875  int nlevels=INTEGER(Rf_getAttrib(elm,Rf_install("nlevels")))[0];
876  for(int i=0;i<x.size();i++){
877  if(map[i]>=0){
878  thetanames[index+map[i]]=nam;
879  if(reversefill)theta[index+map[i]]=x(i);else x(i)=theta[index+map[i]];
880  }
881  }
882  index+=nlevels;
883  }
884  // Auto detect whether we are in "map-mode"
885  SEXP getShape(const char *nam, RObjectTester expectedtype=NULL){
886  SEXP elm=getListElement(parameters,nam);
887  SEXP shape=Rf_getAttrib(elm,Rf_install("shape"));
888  SEXP ans;
889  if(shape==R_NilValue)ans=elm; else ans=shape;
890  RObjectTestExpectedType(ans, expectedtype, nam);
891  return ans;
892  }
893  template<class ArrayType>
894  //ArrayType fillShape(ArrayType &x, const char *nam){
895  ArrayType fillShape(ArrayType x, const char *nam){
896  SEXP elm=getListElement(parameters,nam);
897  SEXP shape=Rf_getAttrib(elm,Rf_install("shape"));
898  if(shape==R_NilValue)fill(x,nam);
899  else fillmap(x,nam);
900  return x;
901  }
902 
903  void fill(Type &x, char const *nam)
904  {
905  pushParname(nam);
906  thetanames[index]=nam;
907  if(reversefill)theta[index++]=x;else x=theta[index++];
908  }
909 
910  Type operator() ();
911 
912  Type evalUserTemplate(){
913  Type ans=this->operator()();
914  /* After evaluating the template, "index" should be equal to the length of "theta".
915  If not, we assume that the "epsilon method" has been requested from R, I.e.
916  that the un-used theta parameters are reserved for an inner product contribution
917  with the numbers reported via ADREPORT. */
918  if(index != theta.size()){
919  PARAMETER_VECTOR( TMB_epsilon_ );
920  ans += ( this->reportvector() * TMB_epsilon_ ).sum();
921  }
922  return ans;
923  }
924 
925 }; // objective_function
926 
956 template<class Type>
957 struct parallel_accumulator{
958  Type result;
959  objective_function<Type>* obj;
960  parallel_accumulator(objective_function<Type>* obj_){
961  result=Type(0);
962  obj=obj_;
963 #ifdef _OPENMP
964  obj->max_parallel_regions=config.nthreads;
965 #endif
966  }
967  inline void operator+=(Type x){
968  if(obj->parallel_region())result+=x;
969  }
970  inline void operator-=(Type x){
971  if(obj->parallel_region())result-=x;
972  }
973  operator Type(){
974  return result;
975  }
976 };
977 
978 
979 #ifndef WITH_LIBTMB
980 
981 #ifdef TMBAD_FRAMEWORK
982 template<class ADFunType>
983 SEXP EvalADFunObjectTemplate(SEXP f, SEXP theta, SEXP control)
984 {
985  if(!Rf_isNewList(control))Rf_error("'control' must be a list");
986  ADFunType* pf;
987  pf=(ADFunType*)R_ExternalPtrAddr(f);
988  int data_changed = getListInteger(control, "data_changed", 0);
989  if (data_changed) {
990  pf->force_update();
991  }
992  int set_tail = getListInteger(control, "set_tail", 0) - 1;
993  if (set_tail == -1) {
994  pf -> unset_tail();
995  } else {
996  std::vector<TMBad::Index> r(1, set_tail);
997  pf -> set_tail(r);
998  }
999  PROTECT(theta=Rf_coerceVector(theta,REALSXP));
1000  int n=pf->Domain();
1001  int m=pf->Range();
1002  if(LENGTH(theta)!=n)Rf_error("Wrong parameter length.");
1003  //R-index -> C-index
1004  int rangecomponent = getListInteger(control, "rangecomponent", 1) - 1;
1005  if(!((0<=rangecomponent)&(rangecomponent<=m-1)))
1006  Rf_error("Wrong range component.");
1007  int order = getListInteger(control, "order");
1008  if((order!=0) & (order!=1) & (order!=2) & (order!=3))
1009  Rf_error("order can be 0, 1, 2 or 3");
1010  //int sparsitypattern = getListInteger(control, "sparsitypattern");
1011  //int dumpstack = getListInteger(control, "dumpstack");
1012  SEXP hessiancols; // Hessian columns
1013  PROTECT(hessiancols=getListElement(control,"hessiancols"));
1014  int ncols=Rf_length(hessiancols);
1015  SEXP hessianrows; // Hessian rows
1016  PROTECT(hessianrows=getListElement(control,"hessianrows"));
1017  int nrows=Rf_length(hessianrows);
1018  if((nrows>0)&(nrows!=ncols))Rf_error("hessianrows and hessianrows must have same length");
1019  vector<size_t> cols(ncols);
1020  vector<size_t> cols0(ncols);
1021  vector<size_t> rows(nrows);
1022  if(ncols>0){
1023  for(int i=0;i<ncols;i++){
1024  cols[i]=INTEGER(hessiancols)[i]-1; //R-index -> C-index
1025  cols0[i]=0;
1026  if(nrows>0)rows[i]=INTEGER(hessianrows)[i]-1; //R-index -> C-index
1027  }
1028  }
1029  std::vector<double> x(REAL(theta), REAL(theta) + LENGTH(theta));
1030 
1031  SEXP res=R_NilValue;
1032  SEXP rangeweight=getListElement(control,"rangeweight");
1033  if(rangeweight!=R_NilValue){
1034  if(LENGTH(rangeweight)!=m)Rf_error("rangeweight must have length equal to range dimension");
1035  std::vector<double> w(REAL(rangeweight),
1036  REAL(rangeweight) + LENGTH(rangeweight));
1037  vector<double> ans = pf->Jacobian(x, w);
1038  res = asSEXP(ans);
1039  UNPROTECT(3);
1040  return res;
1041  }
1042  if(order==3){
1043  Rf_error("Not implemented for TMBad");
1044  // vector<double> w(1);
1045  // w[0]=1;
1046  // if((nrows!=1) | (ncols!=1))Rf_error("For 3rd order derivatives a single hessian coordinate must be specified.");
1047  // pf->ForTwo(x,rows,cols); /* Compute forward directions */
1048  // PROTECT(res=asSEXP(asMatrix(pf->Reverse(3,w),n,3)));
1049  }
1050  if(order==0){
1051  //if(dumpstack)CppAD::traceforward0sweep(1);
1052  std::vector<double> ans = pf->operator()(x);
1053  PROTECT(res=asSEXP(ans));
1054  //if(dumpstack)CppAD::traceforward0sweep(0);
1055  SEXP rangenames=Rf_getAttrib(f,Rf_install("range.names"));
1056  if(LENGTH(res)==LENGTH(rangenames)){
1057  Rf_setAttrib(res,R_NamesSymbol,rangenames);
1058  }
1059  }
1060  if(order==1){
1061  std::vector<double> jvec;
1062  SEXP keepx = getListElement(control, "keepx");
1063  if (keepx != R_NilValue && LENGTH(keepx) > 0) {
1064  SEXP keepy = getListElement(control, "keepy");
1065  std::vector<bool> keep_x(pf->Domain(), false);
1066  std::vector<bool> keep_y(pf->Range(), false);
1067  for (int i=0; i<LENGTH(keepx); i++) {
1068  keep_x[INTEGER(keepx)[i] - 1] = true;
1069  }
1070  for (int i=0; i<LENGTH(keepy); i++) {
1071  keep_y[INTEGER(keepy)[i] - 1] = true;
1072  }
1073  n = LENGTH(keepx);
1074  m = LENGTH(keepy);
1075  jvec = pf->Jacobian(x, keep_x, keep_y);
1076  } else {
1077  jvec = pf->Jacobian(x);
1078  }
1079  // if(doforward)pf->Forward(0,x);
1080  matrix<double> jac(m, n);
1081  int k=0;
1082  for (int i=0; i<m; i++) {
1083  for (int j=0; j<n; j++) {
1084  jac(i, j) = jvec[k];
1085  k++;
1086  }
1087  }
1088  PROTECT( res = asSEXP(jac) );
1089  }
1090  //if(order==2)res=asSEXP(pf->Hessian(x,0),1);
1091  if(order==2){
1092  // if(ncols==0){
1093  // if(sparsitypattern){
1094  // PROTECT(res=asSEXP(HessianSparsityPattern(pf)));
1095  // } else {
1096  // PROTECT(res=asSEXP(asMatrix(pf->Hessian(x,rangecomponent),n,n)));
1097  // }
1098  // }
1099  // else if (nrows==0){
1100  // /* Fixme: the cols0 argument should be user changeable */
1101  // PROTECT(res=asSEXP(asMatrix(pf->RevTwo(x,cols0,cols),n,ncols)));
1102  // }
1103  // else PROTECT(res=asSEXP(asMatrix(pf->ForTwo(x,rows,cols),m,ncols)));
1104  }
1105  UNPROTECT(4);
1106  return res;
1107 } // EvalADFunObjectTemplate
1108 #endif
1109 
1110 #ifdef CPPAD_FRAMEWORK
1111 
1145 template<class ADFunType>
1146 SEXP EvalADFunObjectTemplate(SEXP f, SEXP theta, SEXP control)
1147 {
1148  if(!Rf_isNewList(control))Rf_error("'control' must be a list");
1149  ADFunType* pf;
1150  pf=(ADFunType*)R_ExternalPtrAddr(f);
1151  PROTECT(theta=Rf_coerceVector(theta,REALSXP));
1152  int n=pf->Domain();
1153  int m=pf->Range();
1154  if(LENGTH(theta)!=n)Rf_error("Wrong parameter length.");
1155  // Do forwardsweep ?
1156  int doforward = getListInteger(control, "doforward", 1);
1157  //R-index -> C-index
1158  int rangecomponent = getListInteger(control, "rangecomponent", 1) - 1;
1159  if(!((0<=rangecomponent)&(rangecomponent<=m-1)))
1160  Rf_error("Wrong range component.");
1161  int order = getListInteger(control, "order");
1162  if((order!=0) & (order!=1) & (order!=2) & (order!=3))
1163  Rf_error("order can be 0, 1, 2 or 3");
1164  int sparsitypattern = getListInteger(control, "sparsitypattern");
1165  int dumpstack = getListInteger(control, "dumpstack");
1166  SEXP hessiancols; // Hessian columns
1167  PROTECT(hessiancols=getListElement(control,"hessiancols"));
1168  int ncols=Rf_length(hessiancols);
1169  SEXP hessianrows; // Hessian rows
1170  PROTECT(hessianrows=getListElement(control,"hessianrows"));
1171  int nrows=Rf_length(hessianrows);
1172  if((nrows>0)&(nrows!=ncols))Rf_error("hessianrows and hessianrows must have same length");
1173  vector<size_t> cols(ncols);
1174  vector<size_t> cols0(ncols);
1175  vector<size_t> rows(nrows);
1176  if(ncols>0){
1177  for(int i=0;i<ncols;i++){
1178  cols[i]=INTEGER(hessiancols)[i]-1; //R-index -> C-index
1179  cols0[i]=0;
1180  if(nrows>0)rows[i]=INTEGER(hessianrows)[i]-1; //R-index -> C-index
1181  }
1182  }
1183  vector<double> x = asVector<double>(theta);
1184  SEXP res=R_NilValue;
1185  SEXP rangeweight=getListElement(control,"rangeweight");
1186  if(rangeweight!=R_NilValue){
1187  if(LENGTH(rangeweight)!=m)Rf_error("rangeweight must have length equal to range dimension");
1188  if(doforward)pf->Forward(0,x);
1189  res=asSEXP(pf->Reverse(1,asVector<double>(rangeweight)));
1190  UNPROTECT(3);
1191  return res;
1192  }
1193  if(order==3){
1194  vector<double> w(1);
1195  w[0]=1;
1196  if((nrows!=1) | (ncols!=1))Rf_error("For 3rd order derivatives a single hessian coordinate must be specified.");
1197  pf->ForTwo(x,rows,cols); /* Compute forward directions */
1198  PROTECT(res=asSEXP(asMatrix(pf->Reverse(3,w),n,3)));
1199  }
1200  if(order==0){
1201  if(dumpstack)CppAD::traceforward0sweep(1);
1202  PROTECT(res=asSEXP(pf->Forward(0,x)));
1203  if(dumpstack)CppAD::traceforward0sweep(0);
1204  SEXP rangenames=Rf_getAttrib(f,Rf_install("range.names"));
1205  if(LENGTH(res)==LENGTH(rangenames)){
1206  Rf_setAttrib(res,R_NamesSymbol,rangenames);
1207  }
1208  }
1209  if(order==1){
1210  if(doforward)pf->Forward(0,x);
1211  matrix<double> jac(m, n);
1212  vector<double> u(n);
1213  vector<double> v(m);
1214  v.setZero();
1215  for(int i=0; i<m; i++) {
1216  v[i] = 1.0; u = pf->Reverse(1,v);
1217  v[i] = 0.0;
1218  jac.row(i) = u;
1219  }
1220  PROTECT( res = asSEXP(jac) );
1221  }
1222  //if(order==2)res=asSEXP(pf->Hessian(x,0),1);
1223  if(order==2){
1224  if(ncols==0){
1225  if(sparsitypattern){
1226  PROTECT(res=asSEXP(HessianSparsityPattern(pf)));
1227  } else {
1228  PROTECT(res=asSEXP(asMatrix(pf->Hessian(x,rangecomponent),n,n)));
1229  }
1230  }
1231  else if (nrows==0){
1232  /* Fixme: the cols0 argument should be user changeable */
1233  PROTECT(res=asSEXP(asMatrix(pf->RevTwo(x,cols0,cols),n,ncols)));
1234  }
1235  else PROTECT(res=asSEXP(asMatrix(pf->ForTwo(x,rows,cols),m,ncols)));
1236  }
1237  UNPROTECT(4);
1238  return res;
1239 } // EvalADFunObjectTemplate
1240 #endif
1241 
1243 template <class ADFunType>
1244 void finalize(SEXP x)
1245 {
1246  ADFunType* ptr=(ADFunType*)R_ExternalPtrAddr(x);
1247  if(ptr!=NULL)delete ptr;
1248  memory_manager.CallCFinalizer(x);
1249 }
1250 
1251 #ifdef TMBAD_FRAMEWORK
1252 
1253 TMBad::ADFun< TMBad::ad_aug >* MakeADFunObject_(SEXP data, SEXP parameters,
1254  SEXP report, SEXP control, int parallel_region=-1,
1255  SEXP &info=R_NilValue)
1256 {
1257  typedef TMBad::ad_aug ad;
1258  typedef TMBad::ADFun<ad> adfun;
1259  int returnReport = (control!=R_NilValue) && getListInteger(control, "report");
1260  /* Create objective_function "dummy"-object */
1261  objective_function< ad > F(data,parameters,report);
1262  F.set_parallel_region(parallel_region);
1263  /* Create ADFun pointer.
1264  We have the option to tape either the value returned by the
1265  objective_function template or the vector reported using the
1266  macro "ADREPORT" */
1267  adfun* pf = new adfun();
1268  pf->glob.ad_start();
1269  //TMBad::Independent(F.theta); // In both cases theta is the independent variable
1270  for (int i=0; i<F.theta.size(); i++) F.theta(i).Independent();
1271  if(!returnReport){ // Default case: no ad report - parallel run allowed
1272  vector< ad > y(1);
1273  y[0] = F.evalUserTemplate();
1274  //TMBad::Dependent(y);
1275  for (int i=0; i<y.size(); i++) y[i].Dependent();
1276  } else { // ad report case
1277  F(); // Run through user template (modifies reportvector)
1278  //TMBad::Dependent(F.reportvector.result);
1279  for (int i=0; i<F.reportvector.size(); i++) F.reportvector.result[i].Dependent();
1280  info=F.reportvector.reportnames(); // parallel run *not* allowed
1281  }
1282  pf->glob.ad_stop();
1283  return pf;
1284 }
1285 #endif
1286 
1287 #ifdef CPPAD_FRAMEWORK
1288 
1289 ADFun<double>* MakeADFunObject_(SEXP data, SEXP parameters,
1290  SEXP report, SEXP control, int parallel_region=-1,
1291  SEXP &info=R_NilValue)
1292 {
1293  int returnReport = getListInteger(control, "report");
1294  /* Create objective_function "dummy"-object */
1295  objective_function< AD<double> > F(data,parameters,report);
1296  F.set_parallel_region(parallel_region);
1297  /* Create ADFun pointer.
1298  We have the option to tape either the value returned by the
1299  objective_function template or the vector reported using the
1300  macro "ADREPORT" */
1301  Independent(F.theta); // In both cases theta is the independent variable
1302  ADFun< double >* pf;
1303  if(!returnReport){ // Default case: no ad report - parallel run allowed
1304  vector< AD<double> > y(1);
1305  y[0]=F.evalUserTemplate();
1306  pf = new ADFun< double >(F.theta,y);
1307  } else { // ad report case
1308  F(); // Run through user template (modifies reportvector)
1309  pf = new ADFun< double >(F.theta,F.reportvector());
1310  info=F.reportvector.reportnames(); // parallel run *not* allowed
1311  }
1312  return pf;
1313 }
1314 #endif
1315 
1316 extern "C"
1317 {
1318 
1319 #ifdef TMBAD_FRAMEWORK
1320 
1321  void finalizeADFun(SEXP x)
1322  {
1323  finalize<TMBad::ADFun<TMBad::ad_aug> > (x);
1324  }
1325  void finalizeparallelADFun(SEXP x)
1326  {
1327  finalize<parallelADFun<double> > (x);
1328  }
1329 #endif
1330 
1331 #ifdef CPPAD_FRAMEWORK
1332 
1333  void finalizeADFun(SEXP x)
1334  {
1335  finalize<ADFun<double> > (x);
1336  }
1337  void finalizeparallelADFun(SEXP x)
1338  {
1339  finalize<parallelADFun<double> > (x);
1340  }
1341 #endif
1342 
1343  /* --- MakeADFunObject ----------------------------------------------- */
1344 
1345 #ifdef TMBAD_FRAMEWORK
1346 
1347  SEXP MakeADFunObject(SEXP data, SEXP parameters,
1348  SEXP report, SEXP control)
1349  {
1350  typedef TMBad::ad_aug ad;
1351  typedef TMBad::ADFun<ad> adfun;
1352 
1353  adfun* pf = NULL;
1354  /* Some type checking */
1355  if(!Rf_isNewList(data))Rf_error("'data' must be a list");
1356  if(!Rf_isNewList(parameters))Rf_error("'parameters' must be a list");
1357  if(!Rf_isEnvironment(report))Rf_error("'report' must be an environment");
1358  if(!Rf_isNewList(control))Rf_error("'control' must be a list");
1359  int returnReport = getListInteger(control, "report");
1360 
1361  /* Get the default parameter vector (tiny overhead) */
1362  SEXP par,res=NULL,info;
1363  objective_function< double > F(data,parameters,report);
1364 #ifdef _OPENMP
1365  int n=F.count_parallel_regions(); // Evaluates user template
1366 #else
1367  F.count_parallel_regions(); // Evaluates user template
1368 #endif
1369  if(returnReport && F.reportvector.size()==0){
1370  /* Told to report, but no ADREPORT in template: Get out quickly */
1371  return R_NilValue;
1372  }
1373  PROTECT(par=F.defaultpar());
1374  PROTECT(info=R_NilValue); // Important
1375 
1376  if(_openmp && !returnReport){ // Parallel mode
1377 #ifdef _OPENMP
1378  if(config.trace.parallel)
1379  std::cout << n << " regions found.\n";
1380  if (n==0) n++; // No explicit parallel accumulation
1381  start_parallel(); /* FIXME: NOT NEEDED */
1382  vector< adfun* > pfvec(n);
1383  const char* bad_thread_alloc = NULL;
1384 #pragma omp parallel for num_threads(config.nthreads) if (config.tape.parallel && n>1)
1385  for(int i = 0; i < n; i++) {
1386  TMB_TRY {
1387  pfvec[i] = NULL;
1388  pfvec[i] = MakeADFunObject_(data, parameters, report, control, i, info);
1389  if (config.optimize.instantly) pfvec[i]->optimize();
1390  }
1391  TMB_CATCH {
1392  if (pfvec[i] != NULL) delete pfvec[i];
1393  bad_thread_alloc = excpt.what();
1394  }
1395  }
1396  if (bad_thread_alloc) {
1397  TMB_ERROR_BAD_THREAD_ALLOC;
1398  }
1399 
1400  // FIXME: NOT DONE YET
1401 
1402  parallelADFun<double>* ppf=new parallelADFun<double>(pfvec);
1403  /* Convert parallel ADFun pointer to R_ExternalPtr */
1404  PROTECT(res=R_MakeExternalPtr((void*) ppf,Rf_install("parallelADFun"),R_NilValue));
1405  //R_RegisterCFinalizer(res,finalizeparallelADFun);
1406 #endif
1407  } else { // Serial mode
1408  TMB_TRY{
1409  /* Actual work: tape creation */
1410  pf = NULL;
1411  pf = MakeADFunObject_(data, parameters, report, control, -1, info);
1412  if (config.optimize.instantly) pf->optimize();
1413  }
1414  TMB_CATCH {
1415  if (pf != NULL) delete pf;
1416  TMB_ERROR_BAD_ALLOC;
1417  }
1418  /* Convert ADFun pointer to R_ExternalPtr */
1419  PROTECT(res=R_MakeExternalPtr((void*) pf,Rf_install("ADFun"),R_NilValue));
1420  Rf_setAttrib(res,Rf_install("range.names"),info);
1421  }
1422 
1423  /* Return list of external pointer and default-parameter */
1424  SEXP ans;
1425  Rf_setAttrib(res,Rf_install("par"),par);
1426  PROTECT(ans=ptrList(res));
1427  UNPROTECT(4);
1428 
1429  return ans;
1430  } // MakeADFunObject
1431 #endif
1432 
1433 #ifdef CPPAD_FRAMEWORK
1434 
1435  SEXP MakeADFunObject(SEXP data, SEXP parameters,
1436  SEXP report, SEXP control)
1437  {
1438  ADFun<double>* pf = NULL;
1439  /* Some type checking */
1440  if(!Rf_isNewList(data))Rf_error("'data' must be a list");
1441  if(!Rf_isNewList(parameters))Rf_error("'parameters' must be a list");
1442  if(!Rf_isEnvironment(report))Rf_error("'report' must be an environment");
1443  if(!Rf_isNewList(control))Rf_error("'control' must be a list");
1444  int returnReport = getListInteger(control, "report");
1445 
1446  /* Get the default parameter vector (tiny overhead) */
1447  SEXP par,res=NULL,info;
1448  objective_function< double > F(data,parameters,report);
1449 #ifdef _OPENMP
1450  int n=F.count_parallel_regions(); // Evaluates user template
1451 #else
1452  F.count_parallel_regions(); // Evaluates user template
1453 #endif
1454  if(returnReport && F.reportvector.size()==0){
1455  /* Told to report, but no ADREPORT in template: Get out quickly */
1456  return R_NilValue;
1457  }
1458  PROTECT(par=F.defaultpar());
1459  PROTECT(info=R_NilValue); // Important
1460 
1461  if(_openmp && !returnReport){ // Parallel mode
1462 #ifdef _OPENMP
1463  if(config.trace.parallel)
1464  std::cout << n << " regions found.\n";
1465  if (n==0) n++; // No explicit parallel accumulation
1466  start_parallel(); /* Start threads */
1467  vector< ADFun<double>* > pfvec(n);
1468  const char* bad_thread_alloc = NULL;
1469 #pragma omp parallel for num_threads(config.nthreads) if (config.tape.parallel && n>1)
1470  for(int i=0;i<n;i++){
1471  TMB_TRY {
1472  pfvec[i] = NULL;
1473  pfvec[i] = MakeADFunObject_(data, parameters, report, control, i, info);
1474  if (config.optimize.instantly) pfvec[i]->optimize();
1475  }
1476  TMB_CATCH {
1477  if (pfvec[i] != NULL) delete pfvec[i];
1478  bad_thread_alloc = excpt.what();
1479  }
1480  }
1481  if (bad_thread_alloc) {
1482  TMB_ERROR_BAD_THREAD_ALLOC;
1483  }
1484  parallelADFun<double>* ppf=new parallelADFun<double>(pfvec);
1485  /* Convert parallel ADFun pointer to R_ExternalPtr */
1486  PROTECT(res=R_MakeExternalPtr((void*) ppf,Rf_install("parallelADFun"),R_NilValue));
1487 #endif
1488  } else { // Serial mode
1489  TMB_TRY{
1490  /* Actual work: tape creation */
1491  pf = NULL;
1492  pf = MakeADFunObject_(data, parameters, report, control, -1, info);
1493  if (config.optimize.instantly) pf->optimize();
1494  }
1495  TMB_CATCH {
1496  if (pf != NULL) delete pf;
1497  TMB_ERROR_BAD_ALLOC;
1498  }
1499  /* Convert ADFun pointer to R_ExternalPtr */
1500  PROTECT(res=R_MakeExternalPtr((void*) pf,Rf_install("ADFun"),R_NilValue));
1501  Rf_setAttrib(res,Rf_install("range.names"),info);
1502  }
1503 
1504  /* Return list of external pointer and default-parameter */
1505  SEXP ans;
1506  Rf_setAttrib(res,Rf_install("par"),par);
1507  PROTECT(ans=ptrList(res));
1508  UNPROTECT(4);
1509 
1510  return ans;
1511  } // MakeADFunObject
1512 #endif
1513 
1514  /* --- TransformADFunObject ----------------------------------------------- */
1515 
1516 #ifdef TMBAD_FRAMEWORK
1517 inline int get_num_tapes(SEXP f) {
1518  if (Rf_isNull(f))
1519  return 0;
1520  SEXP tag = R_ExternalPtrTag(f);
1521  if (tag != Rf_install("parallelADFun"))
1522  return 0;
1523  return
1524  ((parallelADFun<double>*) R_ExternalPtrAddr(f))->ntapes;
1525 }
1526 SEXP TransformADFunObjectTemplate(TMBad::ADFun<TMBad::ad_aug>* pf, SEXP control)
1527 {
1528  if (pf == NULL)
1529  Rf_error("Cannot transform '<pointer: (nil)>' (unloaded/reloaded DLL?)");
1530  typedef TMBad::ad_aug ad;
1531  typedef TMBad::ADFun<ad> adfun;
1532  // FIXME: Must require non parallel object !!!
1533  std::string method =
1534  CHAR(STRING_ELT(getListElement(control, "method"), 0));
1535  // Test adfun copy
1536  if (method == "copy") {
1537  *pf = adfun(*pf);
1538  return R_NilValue;
1539  }
1540  if (method == "set_compiled") {
1541  int i = 0;
1542 #ifdef _OPENMP
1543  i = omp_get_thread_num();
1544 #endif
1545  typedef void(*fct_ptr1)(double*);
1546  typedef void(*fct_ptr2)(double*,double*);
1547  pf->glob.forward_compiled =
1548  (fct_ptr1) R_ExternalPtrAddr(VECTOR_ELT(getListElement(control, "forward_compiled"), i));
1549  pf->glob.reverse_compiled =
1550  (fct_ptr2) R_ExternalPtrAddr(VECTOR_ELT(getListElement(control, "reverse_compiled"), i));
1551  return R_NilValue;
1552  }
1553  SEXP random_order = getListElement(control, "random_order");
1554  int nr = (Rf_isNull(random_order) ? 0 : LENGTH(random_order));
1555  std::vector<TMBad::Index> random;
1556  if (nr != 0) {
1557  random = std::vector<TMBad::Index>(INTEGER(random_order),
1558  INTEGER(random_order) + nr);
1559  for (size_t i=0; i<random.size(); i++)
1560  random[i] -= 1 ; // R index -> C index
1561  }
1562  TMB_TRY {
1563  if (method == "remove_random_parameters") {
1564  std::vector<bool> mask(pf->Domain(), true);
1565  for (size_t i = 0; i<random.size(); i++)
1566  mask[random[i]] = false;
1567  pf->glob.inv_index = TMBad::subset(pf->glob.inv_index, mask);
1568  }
1569  else if (method == "laplace") {
1570  SEXP config = getListElement(control, "config");
1571  newton::newton_config cfg(config);
1572  *pf = newton::Laplace_(*pf, random, cfg);
1573  }
1574  else if (method == "marginal_gk") {
1575  TMBad::gk_config cfg;
1576  SEXP config = getListElement(control, "config");
1577  if (!Rf_isNull(config)) {
1578  cfg.adaptive = getListInteger(config, "adaptive", 0);
1579  cfg.debug = getListInteger(config, "debug", 0);
1580  }
1581  *pf = pf -> marginal_gk(random, cfg);
1582  }
1583  else if (method == "marginal_sr") {
1584  SEXP config = getListElement(control, "config");
1585  std::vector<TMBad::sr_grid> grids;
1586  SEXP grid = getListElement(config, "grid");
1587  SEXP random2grid = getListElement(config, "random2grid");
1588  for (int i=0; i<LENGTH(grid); i++) {
1589  SEXP grid_i = VECTOR_ELT(grid, i);
1590  SEXP x = getListElement(grid_i, "x");
1591  SEXP w = getListElement(grid_i, "w");
1592  if (LENGTH(x) != LENGTH(w))
1593  Rf_error("Length of grid$x and grid$w must be equal");
1594  TMBad::sr_grid grid_sr;
1595  grid_sr.x = std::vector<double>(REAL(x), REAL(x) + LENGTH(x));
1596  grid_sr.w = std::vector<double>(REAL(w), REAL(w) + LENGTH(w));
1597  grids.push_back(grid_sr);
1598  }
1599  std::vector<TMBad::Index> r2g(INTEGER(random2grid),
1600  INTEGER(random2grid) + LENGTH(random2grid));
1601  for (size_t i=0; i<r2g.size(); i++)
1602  r2g[i] -= 1 ; // R index -> C index
1603  *pf = pf -> marginal_sr(random, grids, r2g, true);
1604  }
1605  else if (method == "parallelize")
1606  *pf = pf -> parallelize(2);
1607  else if (method == "compress") {
1608  int max_period_size = getListInteger(control, "max_period_size", 1024);
1609  TMBad::compress(pf->glob, max_period_size);
1610  }
1611  else if (method == "compress_and_compile") {
1612 #ifdef HAVE_COMPILE_HPP
1613  int max_period_size = getListInteger(control, "max_period_size", 1024);
1614  TMBad::compress(pf->glob, max_period_size);
1615  // if (config.optimize.instantly) pf->glob.eliminate();
1616  TMBad::compile(pf->glob);
1617 #else
1618  Rf_error("TMBad::compile() is unavailable");
1619 #endif
1620  }
1621  else if (method == "accumulation_tree_split")
1622  pf->glob = accumulation_tree_split(pf->glob, true);
1623  else if (method == "fuse_and_replay") {
1624  pf->glob.set_fuse(true);
1625  pf->replay();
1626  pf->glob.set_fuse(false);
1627  }
1628  else if (method == "reorder_random") {
1629  pf->reorder(random);
1630  }
1631  else if (method == "reorder_sub_expressions") {
1632  TMBad::reorder_sub_expressions(pf->glob);
1633  }
1634  else if (method == "reorder_depth_first") {
1635  TMBad::reorder_depth_first(pf->glob);
1636  }
1637  else if (method == "reorder_temporaries") {
1638  TMBad::reorder_temporaries(pf->glob);
1639  }
1640  else if (method == "parallel_accumulate") {
1641  // Known method - done elsewhere
1642  }
1643  else if (method == "optimize") {
1644  pf->optimize();
1645  } else {
1646  Rf_error("Method unknown: '%s'", method.c_str());
1647  }
1648  }
1649  TMB_CATCH {
1650  TMB_ERROR_BAD_ALLOC;
1651  }
1652  // for (size_t i=0; i<random.size(); i++) random[i] += 1 ; // C index -> R index
1653  // Rf_setAttrib(f, Rf_install("random_order"), asSEXP(random));
1654  return R_NilValue;
1655 }
1657 SEXP TransformADFunObject(SEXP f, SEXP control)
1658 {
1659  if (Rf_isNull(f))
1660  Rf_error("Expected external pointer - got NULL");
1661  SEXP tag = R_ExternalPtrTag(f);
1662  if (tag != Rf_install("ADFun"))
1663  if (tag != Rf_install("parallelADFun"))
1664  Rf_error("Expected ADFun or parallelADFun pointer");
1665  typedef TMBad::ad_aug ad;
1666  typedef TMBad::ADFun<ad> adfun;
1667  if(tag == Rf_install("ADFun")) {
1668  adfun* pf = (adfun*) R_ExternalPtrAddr(f);
1669  TransformADFunObjectTemplate(pf, control);
1670  } else if (tag == Rf_install("parallelADFun")) {
1671  // Warning: Most no meaningful for parallel models!:
1672  // OK : reorder_random etc
1673  // NOT OK : copy, set_compiled, marginal_sr etc
1674  parallelADFun<double>* ppf = (parallelADFun<double>*) R_ExternalPtrAddr(f);
1675  // Apply method for each component except for one special case:
1676  // 'Parallel accumulate'
1677  std::string method =
1678  CHAR(STRING_ELT(getListElement(control, "method"), 0));
1679  if (method == "parallel_accumulate") {
1680  int num_threads = getListInteger(control, "num_threads", 2);
1681  if (num_threads == 1) {
1682  // No need to parallelize
1683  return R_NilValue;
1684  }
1685  if (get_num_tapes(f) > 1) {
1686  // Already parallel (via parallel_accumulator or similar)
1687  return R_NilValue;
1688  }
1689  adfun* pf = (ppf->vecpf)[0]; // One tape - get it
1690  std::vector<adfun> vf = pf->parallel_accumulate(num_threads);
1691  if (config.trace.parallel) {
1692  Rcout << "Autopar work split\n";
1693  for (size_t i=0; i < vf.size(); i++) {
1694  Rcout << "Chunk " << i << ": ";
1695  Rcout << (double) vf[i].glob.opstack.size() / pf->glob.opstack.size() << "\n";
1696  }
1697  }
1698  parallelADFun<double>* new_ppf = new parallelADFun<double>(vf);
1699  delete ppf;
1700  R_SetExternalPtrAddr(f, new_ppf);
1701  return R_NilValue;
1702  }
1703 #ifdef _OPENMP
1704 #pragma omp parallel for num_threads(config.nthreads)
1705 #endif
1706  for (int i=0; i<ppf->ntapes; i++) {
1707  adfun* pf = (ppf->vecpf)[i];
1708  TransformADFunObjectTemplate(pf, control);
1709  }
1710  // Some methods change Domain or Range of individual tapes. This
1711  // is allowed when there is only one tape.
1712  if (ppf->ntapes == 1) {
1713  ppf->domain = (ppf->vecpf)[0]->Domain();
1714  ppf->range = (ppf->vecpf)[0]->Range();
1715  }
1716  // Now, check that it's ok. FIXME: Range() is not checked
1717  for (int i=0; i<ppf->ntapes; i++) {
1718  if (ppf->domain != (ppf->vecpf)[i]->Domain())
1719  Rf_warning("Domain has changed in an invalid way");
1720  }
1721  } else {
1722  Rf_error("Unknown function pointer");
1723  }
1724  return R_NilValue;
1725 }
1726 #endif
1727 
1728 #ifdef CPPAD_FRAMEWORK
1729 
1730 SEXP TransformADFunObject(SEXP f, SEXP control)
1731 {
1732  int mustWork = getListInteger(control, "mustWork", 1);
1733  if (mustWork)
1734  Rf_error("Not supported for CPPAD_FRAMEWORK");
1735  return R_NilValue;
1736 }
1737 #endif
1738 
1739  /* --- InfoADFunObject ---------------------------------------------------- */
1740 
1741 #ifdef TMBAD_FRAMEWORK
1742  SEXP InfoADFunObject(SEXP f) {
1743  typedef TMBad::ad_aug ad;
1744  typedef TMBad::ADFun<ad> adfun;
1745  if (Rf_isNull(f)) Rf_error("Expected external pointer - got NULL");
1746  int num_tapes = get_num_tapes(f);
1747  if (num_tapes >= 2)
1748  Rf_error("'InfoADFunObject' is only available for tapes with one thread");
1749  adfun* pf;
1750  if (num_tapes == 0)
1751  pf = (adfun*) R_ExternalPtrAddr(f);
1752  else {
1753  pf = ( (parallelADFun<double>*) R_ExternalPtrAddr(f) ) -> vecpf[0];
1754  }
1755  SEXP ans, names;
1756  PROTECT(ans = Rf_allocVector(VECSXP, 6));
1757  PROTECT(names = Rf_allocVector(STRSXP, 6));
1758  int i = 0;
1759 #define GET_INFO(EXPR) \
1760  SET_VECTOR_ELT(ans, i, asSEXP(EXPR)); \
1761  SET_STRING_ELT(names, i, Rf_mkChar(#EXPR)); \
1762  i++;
1763  // begin
1764  std::vector<bool> a = pf -> activeDomain();
1765  std::vector<int> ai(a.begin(), a.end());
1766  vector<int> activeDomain(ai);
1767  GET_INFO(activeDomain);
1768  int opstack_size = pf->glob.opstack.size();
1769  GET_INFO(opstack_size);
1770  int values_size = pf->glob.values.size();
1771  GET_INFO(values_size);
1772  int inputs_size = pf->glob.inputs.size();
1773  GET_INFO(inputs_size);
1774  int Domain = pf->Domain();
1775  GET_INFO(Domain);
1776  int Range = pf->Range();
1777  GET_INFO(Range);
1778  // end
1779 #undef GET_INFO
1780  Rf_setAttrib(ans,R_NamesSymbol,names);
1781  UNPROTECT(2);
1782  return ans;
1783  }
1784 #endif
1785 
1786 #ifdef CPPAD_FRAMEWORK
1787  SEXP InfoADFunObject(SEXP f)
1788  {
1789  ADFun<double>* pf;
1790  pf = (ADFun<double>*) R_ExternalPtrAddr(f);
1791  SEXP ans, names;
1792  PROTECT(ans = Rf_allocVector(VECSXP, 12));
1793  PROTECT(names = Rf_allocVector(STRSXP, 12));
1794  int i = 0;
1795 #define GET_MORE_INFO(MEMBER) \
1796  SET_VECTOR_ELT(ans, i, asSEXP(int(pf->MEMBER()))); \
1797  SET_STRING_ELT(names, i, Rf_mkChar(#MEMBER)); \
1798  i++;
1799  GET_MORE_INFO(Domain);
1800  GET_MORE_INFO(Range);
1801  GET_MORE_INFO(size_op);
1802  GET_MORE_INFO(size_op_arg);
1803  GET_MORE_INFO(size_op_seq);
1804  GET_MORE_INFO(size_par);
1805  GET_MORE_INFO(size_order);
1806  GET_MORE_INFO(size_direction);
1807  GET_MORE_INFO(size_text);
1808  GET_MORE_INFO(size_var);
1809  GET_MORE_INFO(size_VecAD);
1810  GET_MORE_INFO(Memory);
1811 #undef GET_MORE_INFO
1812  Rf_setAttrib(ans,R_NamesSymbol,names);
1813  UNPROTECT(2);
1814  return ans;
1815  }
1816 #endif
1817 
1818 #ifdef CPPAD_FRAMEWORK
1819 
1820  SEXP optimizeADFunObject(SEXP f)
1821  {
1822  SEXP tag=R_ExternalPtrTag(f);
1823  if(tag == Rf_install("ADFun")){
1824  ADFun<double>* pf;
1825  pf=(ADFun<double>*)R_ExternalPtrAddr(f);
1826  pf->optimize();
1827  }
1828  if(tag == Rf_install("parallelADFun")){
1829  parallelADFun<double>* pf;
1830  pf=(parallelADFun<double>*)R_ExternalPtrAddr(f);
1831  pf->optimize();
1832  }
1833  return R_NilValue;
1834  }
1835 #endif
1836 
1838  SEXP getTag(SEXP f){
1839  return R_ExternalPtrTag(f);
1840  }
1841 
1842 #ifdef TMBAD_FRAMEWORK
1843  SEXP EvalADFunObject(SEXP f, SEXP theta, SEXP control)
1844  {
1845  typedef TMBad::ad_aug ad;
1846  typedef TMBad::ADFun<ad> adfun;
1847  TMB_TRY {
1848  if(Rf_isNull(f))Rf_error("Expected external pointer - got NULL");
1849  SEXP tag=R_ExternalPtrTag(f);
1850  if(tag == Rf_install("ADFun"))
1851  return EvalADFunObjectTemplate< adfun >(f,theta,control);
1852  if(tag == Rf_install("parallelADFun"))
1853  return EvalADFunObjectTemplate<parallelADFun<double> >(f,theta,control);
1854  Rf_error("NOT A KNOWN FUNCTION POINTER");
1855  }
1856  TMB_CATCH {
1857  TMB_ERROR_BAD_ALLOC;
1858  }
1859  }
1860 #endif
1861 
1862 #ifdef CPPAD_FRAMEWORK
1863  SEXP EvalADFunObject(SEXP f, SEXP theta, SEXP control)
1864  {
1865  TMB_TRY {
1866  if(Rf_isNull(f))Rf_error("Expected external pointer - got NULL");
1867  SEXP tag=R_ExternalPtrTag(f);
1868  if(tag == Rf_install("ADFun"))
1869  return EvalADFunObjectTemplate<ADFun<double> >(f,theta,control);
1870  if(tag == Rf_install("parallelADFun"))
1871  return EvalADFunObjectTemplate<parallelADFun<double> >(f,theta,control);
1872  Rf_error("NOT A KNOWN FUNCTION POINTER");
1873  }
1874  TMB_CATCH {
1875  TMB_ERROR_BAD_ALLOC;
1876  }
1877  }
1878 #endif
1879 
1880 SEXP getSetGlobalPtr(SEXP ptr) {
1881 #ifdef TMBAD_FRAMEWORK
1882  SEXP global_ptr_tag = Rf_install("global_ptr");
1883  if (!Rf_isNull(ptr)) {
1884  SEXP tag = R_ExternalPtrTag(ptr);
1885  if (tag != global_ptr_tag) Rf_error("Invalid pointer type");
1886  TMBad::global_ptr = (TMBad::global**) R_ExternalPtrAddr(ptr);
1887  }
1888  SEXP res = R_MakeExternalPtr( (void*) TMBad::global_ptr, global_ptr_tag, R_NilValue);
1889  return res;
1890 #else
1891  return R_NilValue;
1892 #endif
1893 }
1894 
1895  SEXP tmbad_print(SEXP f, SEXP control) {
1896 #ifdef TMBAD_FRAMEWORK
1897  typedef TMBad::ad_aug ad;
1898  typedef TMBad::ADFun<ad> adfun;
1899  int num_tapes = get_num_tapes(f);
1900  adfun* pf;
1901  if (num_tapes == 0)
1902  pf = (adfun*) R_ExternalPtrAddr(f);
1903  else {
1904  int i = getListInteger(control, "i", 0);
1905  pf = ( (parallelADFun<double>*) R_ExternalPtrAddr(f) ) -> vecpf[i];
1906  }
1907  std::string method =
1908  CHAR(STRING_ELT(getListElement(control, "method"), 0));
1909  if (method == "num_tapes") { // Get number of tapes
1910  return Rf_ScalarInteger(num_tapes);
1911  }
1912  else if (method == "tape") { // Print tape
1913  int depth = getListInteger(control, "depth", 1);
1914  TMBad::global::print_config cfg;
1915  cfg.depth = depth;
1916  pf->glob.print(cfg);
1917  }
1918  else if (method == "dot") { // Print dot format
1919  graph2dot(pf->glob, true, Rcout);
1920  }
1921  else if (method == "inv_index") { // Print member
1922  using TMBad::operator<<;
1923  Rcout << pf->glob.inv_index << "\n";
1924  }
1925  else if (method == "dep_index") { // Print member
1926  using TMBad::operator<<;
1927  Rcout << pf->glob.dep_index << "\n";
1928  }
1929  else if (method == "src") { // Print C src code
1930  TMBad::code_config cfg;
1931  cfg.gpu = false;
1932  cfg.asm_comments = false;
1933  cfg.cout = &Rcout;
1934  *cfg.cout << "#include <cmath>" << std::endl;
1935  *cfg.cout
1936  << "template<class T>T sign(const T &x) { return (x > 0) - (x < 0); }"
1937  << std::endl;
1938  TMBad::global glob = pf->glob; // Invoke deep copy
1939  TMBad::compress(glob);
1940  write_forward(glob, cfg);
1941  write_reverse(glob, cfg);
1942  }
1943  else if (method == "op") {
1944  int name = getListInteger(control, "name", 0);
1945  int address = getListInteger(control, "address", 0);
1946  int input_size = getListInteger(control, "input_size", 0);
1947  int output_size = getListInteger(control, "output_size", 0);
1948  size_t n = pf->glob.opstack.size();
1949  SEXP ans = PROTECT(Rf_allocVector(STRSXP, n));
1950  for (size_t i=0; i<n; i++) {
1951  std::stringstream strm;
1952  if (address) strm << (void*) pf->glob.opstack[i] << " ";
1953  if (name) strm << pf->glob.opstack[i]->op_name() << " ";
1954  if (input_size) strm << pf->glob.opstack[i]->input_size();
1955  if (output_size) strm << pf->glob.opstack[i]->output_size();
1956  const std::string& tmp = strm.str();
1957  SET_STRING_ELT(ans, i, Rf_mkChar(tmp.c_str()));
1958  }
1959  UNPROTECT(1);
1960  return ans;
1961  }
1962  else {
1963  Rf_error("Unknown method: %s", method.c_str());
1964  }
1965 #endif
1966  return R_NilValue;
1967  }
1968 
1969 }
1970 
1971 /* Double interface */
1972 extern "C"
1973 {
1974 
1975  /* How to garbage collect a DoubleFun object pointer */
1976  void finalizeDoubleFun(SEXP x)
1977  {
1978  objective_function<double>* ptr=(objective_function<double>*)R_ExternalPtrAddr(x);
1979  if(ptr!=NULL)delete ptr;
1980  memory_manager.CallCFinalizer(x);
1981  }
1982 
1983  SEXP MakeDoubleFunObject(SEXP data, SEXP parameters, SEXP report, SEXP control)
1984  {
1985  /* Some type checking */
1986  if(!Rf_isNewList(data))Rf_error("'data' must be a list");
1987  if(!Rf_isNewList(parameters))Rf_error("'parameters' must be a list");
1988  if(!Rf_isEnvironment(report))Rf_error("'report' must be an environment");
1989 
1990  /* Create DoubleFun pointer */
1991  objective_function<double>* pF = NULL;
1992  TMB_TRY {
1993  pF = new objective_function<double>(data,parameters,report);
1994  }
1995  TMB_CATCH {
1996  if (pF != NULL) delete pF;
1997  TMB_ERROR_BAD_ALLOC;
1998  }
1999 
2000  /* Convert DoubleFun pointer to R_ExternalPtr */
2001  SEXP res,ans;
2002  PROTECT(res=R_MakeExternalPtr((void*) pF,Rf_install("DoubleFun"),R_NilValue));
2003  PROTECT(ans=ptrList(res));
2004  UNPROTECT(2);
2005  return ans;
2006  }
2007 
2008 
2009  SEXP EvalDoubleFunObject(SEXP f, SEXP theta, SEXP control)
2010  {
2011  TMB_TRY {
2012  int do_simulate = getListInteger(control, "do_simulate");
2013  int get_reportdims = getListInteger(control, "get_reportdims");
2014  objective_function<double>* pf;
2015  pf = (objective_function<double>*) R_ExternalPtrAddr(f);
2016  pf -> sync_data();
2017  PROTECT( theta=Rf_coerceVector(theta,REALSXP) );
2018  int n = pf->theta.size();
2019  if (LENGTH(theta)!=n) Rf_error("Wrong parameter length.");
2020  vector<double> x(n);
2021  for(int i=0;i<n;i++) x[i] = REAL(theta)[i];
2022  pf->theta=x;
2023  /* Since we are actually evaluating objective_function::operator() (not
2024  an ADFun object) we should remember to initialize parameter-index. */
2025  pf->index=0;
2026  pf->parnames.resize(0); // To avoid mem leak.
2027  pf->reportvector.clear();
2028  SEXP res;
2029  GetRNGstate(); /* Get seed from R */
2030  if(do_simulate) pf->set_simulate( true );
2031  PROTECT( res = asSEXP( pf->operator()() ) );
2032  if(do_simulate) {
2033  pf->set_simulate( false );
2034  PutRNGstate(); /* Write seed back to R */
2035  }
2036  if(get_reportdims) {
2037  SEXP reportdims;
2038  PROTECT( reportdims = pf -> reportvector.reportdims() );
2039  Rf_setAttrib( res, Rf_install("reportdims"), reportdims);
2040  UNPROTECT(1);
2041  }
2042  UNPROTECT(2);
2043  return res;
2044  }
2045  TMB_CATCH {
2046  TMB_ERROR_BAD_ALLOC;
2047  }
2048  }
2049 
2053  SEXP getParameterOrder(SEXP data, SEXP parameters, SEXP report, SEXP control)
2054  {
2055  TMB_TRY {
2056  /* Some type checking */
2057  if(!Rf_isNewList(data))Rf_error("'data' must be a list");
2058  if(!Rf_isNewList(parameters))Rf_error("'parameters' must be a list");
2059  if(!Rf_isEnvironment(report))Rf_error("'report' must be an environment");
2060  objective_function<double> F(data,parameters,report);
2061  F(); // Run through user template
2062  return F.parNames();
2063  }
2064  TMB_CATCH {
2065  TMB_ERROR_BAD_ALLOC;
2066  }
2067  }
2068 
2069 } /* Double interface */
2070 
2071 
2072 #ifdef TMBAD_FRAMEWORK
2073 TMBad::ADFun< TMBad::ad_aug >* MakeADGradObject_(SEXP data, SEXP parameters, SEXP report, SEXP control, int parallel_region=-1)
2074 {
2075  typedef TMBad::ad_aug ad;
2076  typedef TMBad::ADFun<ad> adfun;
2077  SEXP f = getListElement(control, "f");
2078  adfun* pf;
2079  bool allocate_new_pf = ( f == R_NilValue );
2080  if ( ! allocate_new_pf ) {
2081  if (parallel_region == -1)
2082  pf = (adfun*) R_ExternalPtrAddr(f);
2083  else
2084  pf = ((parallelADFun<double>*) R_ExternalPtrAddr(f))->vecpf[parallel_region];
2085  } else {
2086  SEXP control_adfun = R_NilValue;
2087  pf = MakeADFunObject_(data, parameters, report, control_adfun, parallel_region);
2088  }
2089  // Optionally skip gradient components (only need 'random' part of gradient)
2090  SEXP random = getListElement(control, "random");
2091  if (random != R_NilValue) {
2092  int set_tail = INTEGER(random)[0] - 1;
2093  std::vector<TMBad::Index> r(1, set_tail);
2094  pf -> set_tail(r);
2095  }
2096  adfun* pgf = new adfun (pf->JacFun());
2097  pf -> unset_tail(); // Not really needed
2098  if (allocate_new_pf) delete pf;
2099  return pgf;
2100 }
2101 #endif
2102 
2103 #ifdef CPPAD_FRAMEWORK
2104 ADFun< double >* MakeADGradObject_(SEXP data, SEXP parameters, SEXP report, SEXP control, int parallel_region=-1)
2105 {
2106  /* Create ADFun pointer */
2107  objective_function< AD<AD<double> > > F(data,parameters,report);
2108  F.set_parallel_region(parallel_region);
2109  int n=F.theta.size();
2110  Independent(F.theta);
2111  vector< AD<AD<double> > > y(1);
2112  y[0]=F.evalUserTemplate();
2113  ADFun<AD<double> > tmp(F.theta,y);
2114  tmp.optimize(); /* Remove 'dead' operations (could result in nan derivatives) */
2115  vector<AD<double> > x(n);
2116  for(int i=0;i<n;i++)x[i]=CppAD::Value(F.theta[i]);
2117  vector<AD<double> > yy(n);
2118  Independent(x);
2119  yy=tmp.Jacobian(x);
2120  ADFun< double >* pf = new ADFun< double >(x,yy);
2121  return pf;
2122 }
2123 #endif
2124 
2125 extern "C"
2126 {
2127 #ifdef TMBAD_FRAMEWORK
2128 
2129  SEXP MakeADGradObject(SEXP data, SEXP parameters, SEXP report, SEXP control)
2130  {
2131  typedef TMBad::ad_aug ad;
2132  typedef TMBad::ADFun<ad> adfun;
2133 
2134  adfun* pf = NULL;
2135  /* Some type checking */
2136  if(!Rf_isNewList(data))Rf_error("'data' must be a list");
2137  if(!Rf_isNewList(parameters))Rf_error("'parameters' must be a list");
2138  if(!Rf_isEnvironment(report))Rf_error("'report' must be an environment");
2139 
2140  /* Get the default parameter vector (tiny overhead) */
2141  SEXP par,res=NULL;
2142  objective_function< double > F(data,parameters,report);
2143 #ifdef _OPENMP
2144  SEXP f = getListElement(control, "f");
2145  int n = get_num_tapes(f);
2146  if (n==0) // No tapes? Count!
2147  n = F.count_parallel_regions(); // Evaluates user template
2148 #else
2149  F.count_parallel_regions(); // Evaluates user template
2150 #endif
2151  PROTECT(par=F.defaultpar());
2152 
2153  if(_openmp){ // Parallel mode
2154 #ifdef _OPENMP
2155  if(config.trace.parallel)
2156  std::cout << n << " regions found.\n";
2157  if (n==0) n++; // No explicit parallel accumulation
2158  start_parallel(); /* Start threads */
2159  vector< adfun* > pfvec(n);
2160  const char* bad_thread_alloc = NULL;
2161 #pragma omp parallel for num_threads(config.nthreads) if (config.tape.parallel && n>1)
2162  for(int i=0;i<n;i++){
2163  TMB_TRY {
2164  pfvec[i] = NULL;
2165  pfvec[i] = MakeADGradObject_(data, parameters, report, control, i);
2166  if (config.optimize.instantly) pfvec[i]->optimize();
2167  }
2168  TMB_CATCH {
2169  if (pfvec[i] != NULL) delete pfvec[i];
2170  bad_thread_alloc = excpt.what();
2171  }
2172  }
2173  if (bad_thread_alloc) {
2174  TMB_ERROR_BAD_THREAD_ALLOC;
2175  }
2176  parallelADFun<double>* ppf=new parallelADFun<double>(pfvec);
2177  /* Convert parallel ADFun pointer to R_ExternalPtr */
2178  PROTECT(res=R_MakeExternalPtr((void*) ppf,Rf_install("parallelADFun"),R_NilValue));
2179  // R_RegisterCFinalizer(res,finalizeparallelADFun);
2180 #endif
2181  } else { // Serial mode
2182  /* Actual work: tape creation */
2183  TMB_TRY {
2184  pf = NULL;
2185  pf = MakeADGradObject_(data, parameters, report, control, -1);
2186  if(config.optimize.instantly)pf->optimize();
2187  }
2188  TMB_CATCH {
2189  if (pf != NULL) delete pf;
2190  TMB_ERROR_BAD_ALLOC;
2191  }
2192  /* Convert ADFun pointer to R_ExternalPtr */
2193  PROTECT(res=R_MakeExternalPtr((void*) pf,Rf_install("ADFun"),R_NilValue));
2194  }
2195 
2196  /* Return ptrList */
2197  SEXP ans;
2198  Rf_setAttrib(res,Rf_install("par"),par);
2199  PROTECT(ans=ptrList(res));
2200  UNPROTECT(3);
2201  return ans;
2202  } // MakeADGradObject
2203 #endif
2204 
2205 #ifdef CPPAD_FRAMEWORK
2206 
2207  SEXP MakeADGradObject(SEXP data, SEXP parameters, SEXP report, SEXP control)
2208  {
2209  ADFun<double>* pf = NULL;
2210  /* Some type checking */
2211  if(!Rf_isNewList(data))Rf_error("'data' must be a list");
2212  if(!Rf_isNewList(parameters))Rf_error("'parameters' must be a list");
2213  if(!Rf_isEnvironment(report))Rf_error("'report' must be an environment");
2214 
2215  /* Get the default parameter vector (tiny overhead) */
2216  SEXP par,res=NULL;
2217  objective_function< double > F(data,parameters,report);
2218 #ifdef _OPENMP
2219  int n=F.count_parallel_regions(); // Evaluates user template
2220 #else
2221  F.count_parallel_regions(); // Evaluates user template
2222 #endif
2223  PROTECT(par=F.defaultpar());
2224 
2225  if(_openmp){ // Parallel mode
2226 #ifdef _OPENMP
2227  if(config.trace.parallel)
2228  std::cout << n << " regions found.\n";
2229  if (n==0) n++; // No explicit parallel accumulation
2230  start_parallel(); /* Start threads */
2231  vector< ADFun<double>* > pfvec(n);
2232  const char* bad_thread_alloc = NULL;
2233 #pragma omp parallel for num_threads(config.nthreads) if (config.tape.parallel && n>1)
2234  for(int i=0;i<n;i++){
2235  TMB_TRY {
2236  pfvec[i] = NULL;
2237  pfvec[i] = MakeADGradObject_(data, parameters, report, control, i);
2238  if (config.optimize.instantly) pfvec[i]->optimize();
2239  }
2240  TMB_CATCH {
2241  if (pfvec[i] != NULL) delete pfvec[i];
2242  bad_thread_alloc = excpt.what();
2243  }
2244  }
2245  if (bad_thread_alloc) {
2246  TMB_ERROR_BAD_THREAD_ALLOC;
2247  }
2248  parallelADFun<double>* ppf=new parallelADFun<double>(pfvec);
2249  /* Convert parallel ADFun pointer to R_ExternalPtr */
2250  PROTECT(res=R_MakeExternalPtr((void*) ppf,Rf_install("parallelADFun"),R_NilValue));
2251 #endif
2252  } else { // Serial mode
2253  /* Actual work: tape creation */
2254  TMB_TRY {
2255  pf = NULL;
2256  pf = MakeADGradObject_(data, parameters, report, control, -1);
2257  if(config.optimize.instantly)pf->optimize();
2258  }
2259  TMB_CATCH {
2260  if (pf != NULL) delete pf;
2261  TMB_ERROR_BAD_ALLOC;
2262  }
2263  /* Convert ADFun pointer to R_ExternalPtr */
2264  PROTECT(res=R_MakeExternalPtr((void*) pf,Rf_install("ADFun"),R_NilValue));
2265  }
2266 
2267  /* Return ptrList */
2268  SEXP ans;
2269  Rf_setAttrib(res,Rf_install("par"),par);
2270  PROTECT(ans=ptrList(res));
2271  UNPROTECT(3);
2272  return ans;
2273  } // MakeADGradObject
2274 #endif
2275 }
2276 
2277 
2284 #ifdef TMBAD_FRAMEWORK
2285 sphess_t< TMBad::ADFun< TMBad::ad_aug > > MakeADHessObject2_(SEXP data, SEXP parameters, SEXP report, SEXP control, int parallel_region=-1)
2286 {
2287  typedef TMBad::ad_aug ad;
2288  typedef TMBad::ADFun<ad> adfun;
2289  typedef sphess_t<adfun> sphess;
2290  SEXP gf = getListElement(control, "gf");
2291  adfun* pgf;
2292  bool allocate_new_pgf = ( gf == R_NilValue );
2293  if ( ! allocate_new_pgf ) {
2294  if (parallel_region == -1)
2295  pgf = (adfun*) R_ExternalPtrAddr(gf);
2296  else
2297  pgf = ((parallelADFun<double>*) R_ExternalPtrAddr(gf))->vecpf[parallel_region];
2298  } else {
2299  SEXP control_adgrad = R_NilValue;
2300  pgf = MakeADGradObject_(data, parameters, report, control_adgrad, parallel_region);
2301  }
2302  if (config.optimize.instantly) pgf->optimize();
2303  int n = pgf->Domain();
2304  std::vector<bool> keepcol(n, true);
2305  SEXP skip = getListElement(control, "skip");
2306  for(int i=0; i<LENGTH(skip); i++) {
2307  keepcol[ INTEGER(skip)[i] - 1 ] = false; // skip is R-index !
2308  }
2309  TMBad::SpJacFun_config spjacfun_cfg;
2310  spjacfun_cfg.index_remap = false;
2311  spjacfun_cfg.compress = config.tmbad.sparse_hessian_compress;
2312  TMBad::Sparse<adfun> h = pgf->SpJacFun(keepcol, keepcol, spjacfun_cfg);
2313  if (allocate_new_pgf) delete pgf;
2314  // NB: Lower triangle, column major =
2315  // Transpose of upper triangle, row major
2316  h.subset_inplace( h.row() <= h.col() ); // Upper triangle, row major
2317  h.transpose_inplace(); // Lower triangle, col major
2318  if (config.optimize.instantly) // Optimize now or later ?
2319  h.optimize();
2320  adfun* phf = new adfun( h );
2321  // Convert h.i and h.j to vector<int>
2322  vector<TMBad::Index> h_i(h.i);
2323  vector<TMBad::Index> h_j(h.j);
2324  sphess ans(phf, h_i.cast<int>(), h_j.cast<int>());
2325  return ans;
2326 } // MakeADHessObject2
2327 #endif
2328 
2335 #ifdef CPPAD_FRAMEWORK
2336 sphess MakeADHessObject2_(SEXP data, SEXP parameters, SEXP report, SEXP control, int parallel_region=-1)
2337 {
2338  /* Some type checking */
2339  if(!Rf_isNewList(data))Rf_error("'data' must be a list");
2340  if(!Rf_isNewList(parameters))Rf_error("'parameters' must be a list");
2341  if(!Rf_isEnvironment(report))Rf_error("'report' must be an environment");
2342 
2343  /* Prepare stuff */
2344  objective_function< AD<AD<AD<double> > > > F(data,parameters,report);
2345  F.set_parallel_region(parallel_region);
2346  int n = F.theta.size();
2347  SEXP skip = getListElement(control, "skip");
2348  vector<bool> keepcol(n); // Scatter for fast lookup
2349  for(int i=0; i<n; i++){
2350  keepcol[i]=true;
2351  }
2352  for(int i=0; i<LENGTH(skip); i++){
2353  keepcol[INTEGER(skip)[i]-1]=false; // skip is R-index !
2354  }
2355 #define KEEP_COL(col) (keepcol[col])
2356 #define KEEP_ROW(row,col) ( KEEP_COL(row) && (row>=col) )
2357 
2358  /* Tape 1: Function R^n -> R */
2359  Independent(F.theta);
2360  vector< AD<AD<AD<double> > > > y(1);
2361  y[0] = F.evalUserTemplate();
2362  ADFun<AD<AD<double> > > tape1(F.theta, y);
2363 
2364  /* Tape 2: Gradient R^n -> R^n (and optimize) */
2365  vector<AD<AD<double> > > xx(n);
2366  for(int i=0; i<n; i++) xx[i] = CppAD::Value(F.theta[i]);
2367  vector<AD<AD<double> > > yy(n);
2368  Independent(xx);
2369  yy = tape1.Jacobian(xx);
2370  ADFun<AD<double > > tape2(xx,yy);
2371  if (config.optimize.instantly) tape2.optimize();
2372 
2373  /* Tape 3: Hessian R^n -> R^m (optimize later) */
2374  tape2.my_init(keepcol);
2375  int colisize;
2376  int m=0; // Count number of non-zeros (m)
2377  for(int i=0; i<int(tape2.colpattern.size()); i++){
2378  colisize = tape2.colpattern[i].size();
2379  if(KEEP_COL(i)){
2380  for(int j=0; j<colisize; j++){
2381  m += KEEP_ROW( tape2.colpattern[i][j] , i);
2382  }
2383  }
2384  }
2385  // Allocate index vectors of non-zero pairs
2386  vector<int> rowindex(m);
2387  vector<int> colindex(m);
2388  // Prepare reverse sweep for Hessian columns
2389  vector<AD<double> > u(n);
2390  vector<AD<double> > v(n);
2391  for(int i = 0; i < n; i++) v[i] = 0.0;
2392  vector<AD<double> > xxx(n);
2393  for(int i=0; i<n; i++) xxx[i]=CppAD::Value(CppAD::Value(F.theta[i]));
2394  vector<AD<double> > yyy(m);
2395  CppAD::vector<int>* icol;
2396  // Do sweeps and fill in non-zero index pairs
2397  Independent(xxx);
2398  tape2.Forward(0, xxx);
2399  int k=0;
2400  for(int i = 0; i < n; i++){
2401  if (KEEP_COL(i)) {
2402  tape2.myReverse(1, v, i /*range comp*/, u /*domain*/);
2403  icol = &tape2.colpattern[i];
2404  for(int j=0; j<int(icol->size()); j++){
2405  if(KEEP_ROW( icol->operator[](j), i )){
2406  rowindex[k] = icol->operator[](j);
2407  colindex[k] = i;
2408  yyy[k] = u[icol->operator[](j)];
2409  k++;
2410  }
2411  }
2412  }
2413  }
2414  ADFun< double >* ptape3 = new ADFun< double >;
2415  ptape3->Dependent(xxx,yyy);
2416  sphess ans(ptape3, rowindex, colindex);
2417  return ans;
2418 } // MakeADHessObject2
2419 #endif
2420 
2421 // kasper: Move to new file e.g. "convert.hpp"
2422 template <class ADFunType>
2424 SEXP asSEXP(const sphess_t<ADFunType> &H, const char* tag)
2425 {
2426  SEXP par;
2427  par=R_NilValue;
2428  /* Convert ADFun pointer to R_ExternalPtr */
2429  SEXP res;
2430  PROTECT( res = R_MakeExternalPtr((void*) H.pf, Rf_install(tag), R_NilValue) );
2431  /* Return list */
2432  SEXP ans;
2433  /* Implicitly protected temporaries */
2434  SEXP par_symbol = Rf_install("par");
2435  SEXP i_symbol = Rf_install("i");
2436  SEXP j_symbol = Rf_install("j");
2437  Rf_setAttrib(res, par_symbol, par);
2438  Rf_setAttrib(res, i_symbol, asSEXP(H.i));
2439  Rf_setAttrib(res, j_symbol, asSEXP(H.j));
2440  PROTECT(ans=ptrList(res));
2441  UNPROTECT(2);
2442  return ans;
2443 }
2444 
2445 
2446 extern "C"
2447 {
2448 
2449 #ifdef TMBAD_FRAMEWORK
2450 #ifdef _OPENMP
2451  SEXP MakeADHessObject2(SEXP data, SEXP parameters, SEXP report, SEXP control){
2452  typedef TMBad::ad_aug ad;
2453  typedef TMBad::ADFun<ad> adfun;
2454  typedef sphess_t<adfun> sphess;
2455  if(config.trace.parallel)
2456  std::cout << "Count num parallel regions\n";
2457  objective_function< double > F(data,parameters,report);
2458  SEXP gf = getListElement(control, "gf");
2459  int n = get_num_tapes(gf);
2460  if (n==0) // No tapes? Count!
2461  n = F.count_parallel_regions(); // Evaluates user template
2462  if(config.trace.parallel)
2463  std::cout << n << " regions found.\n";
2464  if (n==0) n++; // No explicit parallel accumulation
2465  start_parallel(); /* FIXME: not needed */
2466  /* parallel test */
2467  const char* bad_thread_alloc = NULL;
2468  vector<sphess*> Hvec(n);
2469 #pragma omp parallel for num_threads(config.nthreads) if (config.tape.parallel && n>1)
2470  for (int i=0; i<n; i++) {
2471  TMB_TRY {
2472  Hvec[i] = NULL;
2473  Hvec[i] = new sphess( MakeADHessObject2_(data, parameters, report, control, i) );
2474  //optimizeTape( Hvec[i]->pf );
2475  }
2476  TMB_CATCH {
2477  if (Hvec[i] != NULL) {
2478  delete Hvec[i]->pf;
2479  delete Hvec[i];
2480  }
2481  bad_thread_alloc = excpt.what();
2482  }
2483  }
2484  if (bad_thread_alloc) {
2485  TMB_ERROR_BAD_THREAD_ALLOC;
2486  }
2487  parallelADFun<double>* tmp=new parallelADFun<double>(Hvec);
2488  return asSEXP(tmp->convert(),"parallelADFun");
2489  } // MakeADHessObject2
2490 #else
2491  SEXP MakeADHessObject2(SEXP data, SEXP parameters, SEXP report, SEXP control){
2492  typedef TMBad::ad_aug ad;
2493  typedef TMBad::ADFun<ad> adfun;
2494  typedef sphess_t<adfun> sphess;
2495  sphess* pH = NULL;
2496  SEXP ans;
2497  TMB_TRY {
2498  pH = new sphess( MakeADHessObject2_(data, parameters, report, control, -1) );
2499  //optimizeTape( pH->pf );
2500  ans = asSEXP(*pH, "ADFun");
2501  }
2502  TMB_CATCH {
2503  if (pH != NULL) {
2504  delete pH->pf;
2505  delete pH;
2506  }
2507  TMB_ERROR_BAD_ALLOC;
2508  }
2509  delete pH;
2510  return ans;
2511  } // MakeADHessObject2
2512 #endif
2513 #endif
2514 
2515 #ifdef CPPAD_FRAMEWORK
2516 #ifdef _OPENMP
2517  SEXP MakeADHessObject2(SEXP data, SEXP parameters, SEXP report, SEXP control){
2518  if(config.trace.parallel)
2519  std::cout << "Count num parallel regions\n";
2520  objective_function< double > F(data,parameters,report);
2521  int n=F.count_parallel_regions();
2522  if(config.trace.parallel)
2523  std::cout << n << " regions found.\n";
2524  if (n==0) n++; // No explicit parallel accumulation
2525 
2526  start_parallel(); /* Start threads */
2527 
2528  /* parallel test */
2529  const char* bad_thread_alloc = NULL;
2530  vector<sphess*> Hvec(n);
2531 #pragma omp parallel for num_threads(config.nthreads) if (config.tape.parallel && n>1)
2532  for (int i=0; i<n; i++) {
2533  TMB_TRY {
2534  Hvec[i] = NULL;
2535  Hvec[i] = new sphess( MakeADHessObject2_(data, parameters, report, control, i) );
2536  optimizeTape( Hvec[i]->pf );
2537  }
2538  TMB_CATCH {
2539  if (Hvec[i] != NULL) {
2540  delete Hvec[i]->pf;
2541  delete Hvec[i];
2542  }
2543  bad_thread_alloc = excpt.what();
2544  }
2545  }
2546  if (bad_thread_alloc) {
2547  TMB_ERROR_BAD_THREAD_ALLOC;
2548  }
2549  parallelADFun<double>* tmp=new parallelADFun<double>(Hvec);
2550  for(int i=0; i<n; i++) {
2551  delete Hvec[i];
2552  }
2553  // Adds finalizer for 'tmp' !!! (so, don't delete tmp...)
2554  SEXP ans = asSEXP(tmp->convert(),"parallelADFun");
2555  return ans;
2556  } // MakeADHessObject2
2557 #else
2558  SEXP MakeADHessObject2(SEXP data, SEXP parameters, SEXP report, SEXP control){
2559  sphess* pH = NULL;
2560  SEXP ans;
2561  TMB_TRY {
2562  pH = new sphess( MakeADHessObject2_(data, parameters, report, control, -1) );
2563  optimizeTape( pH->pf );
2564  ans = asSEXP(*pH, "ADFun");
2565  }
2566  TMB_CATCH {
2567  if (pH != NULL) {
2568  delete pH->pf;
2569  delete pH;
2570  }
2571  TMB_ERROR_BAD_ALLOC;
2572  }
2573  delete pH;
2574  return ans;
2575  } // MakeADHessObject2
2576 #endif
2577 #endif
2578 }
2579 
2580 extern "C"
2581 {
2582 
2583 #ifdef TMBAD_FRAMEWORK
2584  SEXP usingAtomics(){
2585  SEXP ans;
2586  PROTECT(ans = Rf_allocVector(INTSXP,1));
2587  INTEGER(ans)[0] = 1; // TMBAD doesn't benefit from knowing if 'false'
2588  UNPROTECT(1);
2589  return ans;
2590  }
2591 #endif
2592 
2593 #ifdef CPPAD_FRAMEWORK
2594  SEXP usingAtomics(){
2595  SEXP ans;
2596  PROTECT(ans = Rf_allocVector(INTSXP,1));
2597  INTEGER(ans)[0] = atomic::atomicFunctionGenerated;
2598  UNPROTECT(1);
2599  return ans;
2600  }
2601 #endif
2602 
2603  SEXP getFramework() {
2604  // ans
2605  SEXP ans;
2606 #ifdef TMBAD_FRAMEWORK
2607  ans = Rf_mkString("TMBad");
2608 #elif defined(CPPAD_FRAMEWORK)
2609  ans = Rf_mkString("CppAD");
2610 #else
2611  ans = Rf_mkString("Unknown");
2612 #endif
2613  PROTECT(ans);
2614  // openmp_sym (Not strictly necessary to PROTECT)
2615  SEXP openmp_sym = Rf_install("openmp");
2616  PROTECT(openmp_sym);
2617  // openmp_res
2618  SEXP openmp_res;
2619 #ifdef _OPENMP
2620  openmp_res = Rf_ScalarLogical(1);
2621 #else
2622  openmp_res = Rf_ScalarLogical(0);
2623 #endif
2624  PROTECT(openmp_res);
2625  // Assemble
2626  Rf_setAttrib(ans, openmp_sym, openmp_res);
2627  UNPROTECT(2);
2628  // Add more stuff
2629 #ifdef TMBAD_FRAMEWORK
2630  SEXP index_size_sym = Rf_install("sizeof(Index)");
2631  PROTECT(index_size_sym);
2632  SEXP index_size = Rf_ScalarInteger(sizeof(TMBad::Index));
2633  PROTECT(index_size);
2634  Rf_setAttrib(ans, index_size_sym, index_size);
2635  UNPROTECT(2);
2636 #endif
2637  UNPROTECT(1); // ans
2638  return ans;
2639  }
2640 }
2641 
2642 extern "C"
2643 {
2644  void tmb_forward(SEXP f, const Eigen::VectorXd &x, Eigen::VectorXd &y) {
2645 #ifdef CPPAD_FRAMEWORK
2646  SEXP tag=R_ExternalPtrTag(f);
2647  if(tag == Rf_install("ADFun")) {
2648  ADFun<double>* pf;
2649  pf = (ADFun<double>*) R_ExternalPtrAddr(f);
2650  y = pf->Forward(0, x);
2651  } else
2652  if(tag == Rf_install("parallelADFun")) {
2653  parallelADFun<double>* pf;
2654  pf = (parallelADFun<double>*) R_ExternalPtrAddr(f);
2655  y = pf->Forward(0, x);
2656  } else
2657  Rf_error("Unknown function pointer");
2658 #endif
2659 #ifdef TMBAD_FRAMEWORK
2660  typedef TMBad::ad_aug ad;
2661  typedef TMBad::ADFun<ad> adfun;
2662  SEXP tag=R_ExternalPtrTag(f);
2663  if(tag == Rf_install("ADFun")) {
2664  adfun* pf = (adfun*) R_ExternalPtrAddr(f);
2665  y = pf->forward(x);
2666  } else
2667  if(tag == Rf_install("parallelADFun")) {
2668  parallelADFun<double>* pf;
2669  pf = (parallelADFun<double>*) R_ExternalPtrAddr(f);
2670  y = pf->forward(x);
2671  } else
2672  Rf_error("Unknown function pointer");
2673 #endif
2674  }
2675  void tmb_reverse(SEXP f, const Eigen::VectorXd &v, Eigen::VectorXd &y) {
2676 #ifdef CPPAD_FRAMEWORK
2677  SEXP tag=R_ExternalPtrTag(f);
2678  if(tag == Rf_install("ADFun")) {
2679  ADFun<double>* pf;
2680  pf = (ADFun<double>*) R_ExternalPtrAddr(f);
2681  y = pf->Reverse(1, v);
2682  } else
2683  if(tag == Rf_install("parallelADFun")) {
2684  parallelADFun<double>* pf;
2685  pf = (parallelADFun<double>*) R_ExternalPtrAddr(f);
2686  y = pf->Reverse(1, v);
2687  } else
2688  Rf_error("Unknown function pointer");
2689 #endif
2690 #ifdef TMBAD_FRAMEWORK
2691  typedef TMBad::ad_aug ad;
2692  typedef TMBad::ADFun<ad> adfun;
2693  SEXP tag=R_ExternalPtrTag(f);
2694  if(tag == Rf_install("ADFun")) {
2695  adfun* pf = (adfun*) R_ExternalPtrAddr(f);
2696  y = pf->reverse(v);
2697  } else
2698  if(tag == Rf_install("parallelADFun")) {
2699  parallelADFun<double>* pf;
2700  pf = (parallelADFun<double>*) R_ExternalPtrAddr(f);
2701  y = pf->reverse(v);
2702  } else
2703  Rf_error("Unknown function pointer");
2704 #endif
2705  }
2706 }
2707 
2708 #endif /* #ifndef WITH_LIBTMB */
2709 
2710 
2711 
2712 
2713 
2714 #ifdef WITH_LIBTMB
2715 
2716 template class objective_function<double>;
2717 #ifdef CPPAD_FRAMEWORK
2718 template class objective_function<AD<double> >;
2719 template class objective_function<AD<AD<double> > >;
2720 template class objective_function<AD<AD<AD<double> > > >;
2721 #endif
2722 #ifdef TMBAD_FRAMEWORK
2723 template class objective_function<TMBad::ad_aug>;
2724 #endif
2725 
2726 extern "C"
2727 {
2728  SEXP MakeADFunObject(SEXP data, SEXP parameters, SEXP report, SEXP control);
2729  SEXP InfoADFunObject(SEXP f);
2730  SEXP tmbad_print(SEXP f, SEXP control);
2731  SEXP optimizeADFunObject(SEXP f);
2732  SEXP EvalADFunObject(SEXP f, SEXP theta, SEXP control);
2733  SEXP MakeDoubleFunObject(SEXP data, SEXP parameters, SEXP report, SEXP control);
2734  SEXP EvalDoubleFunObject(SEXP f, SEXP theta, SEXP control);
2735  SEXP getParameterOrder(SEXP data, SEXP parameters, SEXP report, SEXP control);
2736  SEXP MakeADGradObject(SEXP data, SEXP parameters, SEXP report, SEXP control);
2737  SEXP MakeADHessObject2(SEXP data, SEXP parameters, SEXP report, SEXP control);
2738  SEXP usingAtomics();
2739  SEXP getFramework();
2740  SEXP getSetGlobalPtr(SEXP ptr);
2741  SEXP TransformADFunObject(SEXP f, SEXP control);
2742  void tmb_forward(SEXP f, const Eigen::VectorXd &x, Eigen::VectorXd &y);
2743  void tmb_reverse(SEXP f, const Eigen::VectorXd &v, Eigen::VectorXd &y);
2744 }
2745 
2746 #endif /* #ifdef WITH_LIBTMB */
2747 
2748 /* Register native routines (see 'Writing R extensions'). Especially
2749  relevant to avoid symbol lookup overhead for those routines that
2750  are called many times e.g. EvalADFunObject. */
2751 extern "C"{
2752  /* Some string utilities */
2753 #define xstringify(s) stringify(s)
2754 #define stringify(s) #s
2755  /* May be used as part of custom calldef tables */
2756 #define TMB_CALLDEFS \
2757  {"MakeADFunObject", (DL_FUNC) &MakeADFunObject, 4}, \
2758  {"FreeADFunObject", (DL_FUNC) &FreeADFunObject, 1}, \
2759  {"InfoADFunObject", (DL_FUNC) &InfoADFunObject, 1}, \
2760  {"tmbad_print", (DL_FUNC) &tmbad_print, 2}, \
2761  {"EvalADFunObject", (DL_FUNC) &EvalADFunObject, 3}, \
2762  {"TransformADFunObject",(DL_FUNC) &TransformADFunObject,2}, \
2763  {"MakeDoubleFunObject", (DL_FUNC) &MakeDoubleFunObject, 4}, \
2764  {"EvalDoubleFunObject", (DL_FUNC) &EvalDoubleFunObject, 3}, \
2765  {"getParameterOrder", (DL_FUNC) &getParameterOrder, 4}, \
2766  {"MakeADGradObject", (DL_FUNC) &MakeADGradObject, 4}, \
2767  {"MakeADHessObject2", (DL_FUNC) &MakeADHessObject2, 4}, \
2768  {"usingAtomics", (DL_FUNC) &usingAtomics, 0}, \
2769  {"getFramework", (DL_FUNC) &getFramework, 0}, \
2770  {"getSetGlobalPtr", (DL_FUNC) &getSetGlobalPtr, 1}, \
2771  {"TMBconfig", (DL_FUNC) &TMBconfig, 2}
2772  /* May be used as part of custom R_init function
2773  C-callable routines (PACKAGE is 'const char*') */
2774 #define TMB_CCALLABLES(PACKAGE) \
2775  R_RegisterCCallable(PACKAGE, "tmb_forward", (DL_FUNC) &tmb_forward); \
2776  R_RegisterCCallable(PACKAGE, "tmb_reverse", (DL_FUNC) &tmb_reverse);
2777  /* Default (optional) calldef table. */
2778 #ifdef TMB_LIB_INIT
2779 #include <R_ext/Rdynload.h>
2780 static R_CallMethodDef CallEntries[] = {
2781  TMB_CALLDEFS
2782  ,
2783  /* User's R_unload_lib function must also be registered (because we
2784  disable dynamic lookup - see below). The unload function is
2785  mainly useful while developing models in order to clean up
2786  external pointers without restarting R. Should not be used by TMB
2787  dependent packages. */
2788 #ifdef LIB_UNLOAD
2789  {xstringify(LIB_UNLOAD), (DL_FUNC) &LIB_UNLOAD, 1},
2790 #endif
2791  /* End of table */
2792  {NULL, NULL, 0}
2793 };
2794 void TMB_LIB_INIT(DllInfo *dll){
2795  R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
2796  R_useDynamicSymbols(dll, (Rboolean)FALSE);
2797  // Example: TMB_LIB_INIT = R_init_mypkg
2798  // ^
2799  // +-------+
2800  // ^
2801  TMB_CCALLABLES(&(xstringify(TMB_LIB_INIT)[7]));
2802 }
2803 #endif /* #ifdef TMB_LIB_INIT */
2804 #undef xstringify
2805 #undef stringify
2806 }
data_indicator::cdf_upper
VT cdf_upper
Logarithm of upper CDF
Definition: tmb_core.hpp:445
TMBad::subset
std::vector< T > subset(const std::vector< T > &x, const std::vector< bool > &y)
Vector subset by boolean mask.
Definition: graph_transform.hpp:40
vector
Vector class used by TMB.
Definition: vector.hpp:17
TMBad::reorder_temporaries
void reorder_temporaries(global &glob)
Re-order computational graph to make it more compressible.
Definition: TMBad.cpp:567
tmbutils::array
Array class used by TMB.
Definition: tmbutils.hpp:23
data_indicator::segment
data_indicator segment(int pos, int n)
Extract segment of indicator vector or array.
Definition: tmb_core.hpp:478
data_indicator::data_indicator
data_indicator()
Default CTOR.
Definition: tmb_core.hpp:451
TMBad::reorder_sub_expressions
void reorder_sub_expressions(global &glob)
Re-order computational graph to make it more compressible.
Definition: TMBad.cpp:537
TMBad::reorder_depth_first
void reorder_depth_first(global &glob)
Depth-first reordering of computational graph.
Definition: TMBad.cpp:595
newton::newton_config
Newton configuration parameters.
Definition: newton.hpp:637
TMBad::ADFun
Automatic differentiation function object.
Definition: TMBad.hpp:117
TMBad::global::print_config
Configuration of print method.
Definition: global.hpp:1479
TMBad::global::ad_aug
Augmented AD type.
Definition: global.hpp:2831
matrix
Matrix class used by TMB.
Definition: vector.hpp:101
data_indicator::osa_active
bool osa_active()
Are we inside a oneStepPredict calculation?
Definition: tmb_core.hpp:510
TMBad::ADFun::reorder
void reorder(std::vector< Index > last)
Reorder computational graph to allow quick updates of selected inputs.
Definition: TMBad.hpp:237
config_struct::optimize
bool optimize
Trace tape optimization.
Definition: config.hpp:30
data_indicator
Utilities for OSA residuals.
Definition: tmb_core.hpp:441
TMBad::global
Struct defining the main AD context.
Definition: global.hpp:797
TMBad::global::ad_aug::Value
Scalar Value() const
Return the underlying scalar value of this ad_aug.
Definition: TMBad.cpp:2188
data_indicator::data_indicator
data_indicator(VT obs, bool init_one=false)
Construct from observation vector.
Definition: tmb_core.hpp:456
config_struct::autopar
bool autopar
Enable automatic parallization (if OpenMP is enabled) ?
Definition: config.hpp:49
data_indicator::order
vector< int > order()
Get order in which the one step conditionals will be requested by oneStepPredict. ...
Definition: tmb_core.hpp:489
PARAMETER_VECTOR
#define PARAMETER_VECTOR(name)
Get parameter vector from R and declare it as vector<Type>
Definition: tmb_core.hpp:220
TMBad::SpJacFun_config
Configuration parameters for SpJacFun()
Definition: TMBad.hpp:77
asSEXP
SEXP asSEXP(const matrix< Type > &a)
Convert TMB matrix, vector, scalar or int to R style.
Definition: convert.hpp:30
config_struct::sparse_hessian_compress
bool sparse_hessian_compress
Reduce memory of sparse hessian even if reducing speed ?
Definition: config.hpp:45
TMBad::ADFun::replay
void replay()
Replay this operation sequence to a new sequence.
Definition: TMBad.hpp:750
data_indicator::cdf_lower
VT cdf_lower
Logarithm of lower CDF
Definition: tmb_core.hpp:443
sum
Type sum(Vector< Type > x)
Definition: convenience.hpp:33
TMBad::sr_grid
Utilility class for sequential_reduction.
Definition: graph_transform.hpp:520
asMatrix
matrix< Type > asMatrix(const vector< Type > &x, int nr, int nc)
Vector <-> Matrix conversion (for row-major matrices)
Definition: convert.hpp:120
TMBad::ADFun::parallel_accumulate
std::vector< ADFun > parallel_accumulate(size_t num_threads)
Parallel split this operation sequence Split function f:R^n->R by its accumulation tree...
Definition: TMBad.hpp:717
config_struct::nthreads
int nthreads
Number of OpenMP threads to use (if OpenMP is enabled)
Definition: config.hpp:50
parallel_accumulator
Helper to manage parallel accumulation.
Definition: tmb_core.hpp:957
TMBad::ADFun::optimize
void optimize()
Tape optimizer.
Definition: TMBad.hpp:195
data_indicator::ord
vector< int > ord
Subset argument (zero-based) passed from oneStepPredict. See data_indicator::order().
Definition: tmb_core.hpp:447
config_struct::parallel
bool parallel
Trace info from parallel for loops.
Definition: config.hpp:29
data_indicator::fill
void fill(vector< Type > p, SEXP ord_)
Fill with parameter vector.
Definition: tmb_core.hpp:464
data_indicator::osa_flag
bool osa_flag
Flag to tell if OSA calculation is active. See data_indicator::osa_active().
Definition: tmb_core.hpp:449
License: GPL v2