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){
162 /* Drop out */
163 return;
164 }
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
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>
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 }
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 }
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 }
486 return ans;
487 }
490 int n = this->size();
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 }
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();
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);
561 vector<int> tmp = (pf->RevSparseHes(n,Py)).
template cast<int>();
563 }
564
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
582 dim << x.rows(), x.cols();
583 return dim;
584 }
586 return x.dim;
587 }
588 template<class Other> // i.e. vector or expression
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){
605 xvec[0] = x;
606 push(xvec, name);
607 }
608 // Eval: cast to vector<Type>
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;
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;
663 report_stack<Type> reportvector;
664 bool reversefill; // used to find the parameter order in user template (not anymore - use pushParname instead)
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 */
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
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. */
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 }
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()){
920 ans += ( this->reportvector() * TMB_epsilon_ ).
sum();
921 }
922 return ans;
923 }
924
925 }; // objective_function
926
956 template<class Type>
958 Type result;
959 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");
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));
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);
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");
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 }
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){
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 */
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);
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 */
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
1254 SEXP report, SEXP control, int parallel_region=-1,
1255 SEXP &info=R_NilValue)
1256 {
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
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
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 {
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
1379 std::cout << n << " regions found.\n";
1380 if (n==0) n++; // No explicit parallel accumulation
1381 start_parallel(); /* FIXME: NOT NEEDED */
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
1464 std::cout << n << " regions found.\n";
1465 if (n==0) n++; // No explicit parallel accumulation
1466 start_parallel(); /* Start threads */
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 }
1527 {
1528 if (pf == NULL)
1529 Rf_error("Cannot transform '<pointer: (nil)>' (unloaded/reloaded DLL?)");
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");
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");
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);
1626 pf->glob.set_fuse(false);
1627 }
1628 else if (method == "reorder_random") {
1630 }
1631 else if (method == "reorder_sub_expressions") {
1633 }
1634 else if (method == "reorder_depth_first") {
1636 }
1637 else if (method == "reorder_temporaries") {
1639 }
1640 else if (method == "parallel_accumulate") {
1641 // Known method - done elsewhere
1642 }
1643 else if (method == "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");
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
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) {
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());
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);
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 {
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
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);
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;
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.");
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
2074 {
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);
2112 y[0]=F.evalUserTemplate();
2113 ADFun<AD<double> > tmp(F.theta,y);
2114 tmp.optimize(); /* Remove 'dead' operations (could result in nan derivatives) */
2116 for(int i=0;i<n;i++)x[i]=CppAD::Value(F.theta[i]);
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 {
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
2156 std::cout << n << " regions found.\n";
2157 if (n==0) n++; // No explicit parallel accumulation
2158 start_parallel(); /* Start threads */
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
2228 std::cout << n << " regions found.\n";
2229 if (n==0) n++; // No explicit parallel accumulation
2230 start_parallel(); /* Start threads */
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 {
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 }
2310 spjacfun_cfg.index_remap = false;
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>
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");
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);
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) */
2366 for(int i=0; i<n; i++) xx[i] = CppAD::Value(F.theta[i]);
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
2388 // Prepare reverse sweep for Hessian columns
2391 for(int i = 0; i < n; i++) v[i] = 0.0;
2393 for(int i=0; i<n; i++) xxx[i]=CppAD::Value(CppAD::Value(F.theta[i]));
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){
2454 typedef sphess_t<adfun> sphess;
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
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;
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){
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){
2519 std::cout << "Count num parallel regions\n";
2520 objective_function< double > F(data,parameters,report);
2521 int n=F.count_parallel_regions();
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;
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
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
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 }
VT cdf_upper
Logarithm of upper CDF
std::vector< T > subset(const std::vector< T > &x, const std::vector< bool > &y)
Vector subset by boolean mask.
Vector class used by TMB.
void reorder_temporaries(global &glob)
Re-order computational graph to make it more compressible.
data_indicator segment(int pos, int n)
Extract segment of indicator vector or array.
data_indicator()
Default CTOR.
void reorder_sub_expressions(global &glob)
Re-order computational graph to make it more compressible.
void reorder_depth_first(global &glob)
Depth-first reordering of computational graph.
Newton configuration parameters.
Automatic differentiation function object.
Configuration of print method.
Matrix class used by TMB.
bool osa_active()
Are we inside a oneStepPredict calculation?
void reorder(std::vector< Index > last)
Reorder computational graph to allow quick updates of selected inputs.
bool optimize
Trace tape optimization.
Utilities for OSA residuals.
Struct defining the main AD context.
Scalar Value() const
Return the underlying scalar value of this ad_aug.
data_indicator(VT obs, bool init_one=false)
Construct from observation vector.
bool autopar
Enable automatic parallization (if OpenMP is enabled) ?
vector< int > order()
Get order in which the one step conditionals will be requested by oneStepPredict. ...
#define PARAMETER_VECTOR(name)
Get parameter vector from R and declare it as vector<Type>
Configuration parameters for SpJacFun()
SEXP asSEXP(const matrix< Type > &a)
Convert TMB matrix, vector, scalar or int to R style.
bool sparse_hessian_compress
Reduce memory of sparse hessian even if reducing speed ?
void replay()
Replay this operation sequence to a new sequence.
VT cdf_lower
Logarithm of lower CDF
Type sum(Vector< Type > x)
Utilility class for sequential_reduction.
matrix< Type > asMatrix(const vector< Type > &x, int nr, int nc)
Vector <-> Matrix conversion (for row-major matrices)
std::vector< ADFun > parallel_accumulate(size_t num_threads)
Parallel split this operation sequence Split function f:R^n->R by its accumulation tree...
int nthreads
Number of OpenMP threads to use (if OpenMP is enabled)
Helper to manage parallel accumulation.
void optimize()
Tape optimizer.
vector< int > ord
Subset argument (zero-based) passed from oneStepPredict. See data_indicator::order().
bool parallel
Trace info from parallel for loops.
void fill(vector< Type > p, SEXP ord_)
Fill with parameter vector.
bool osa_flag
Flag to tell if OSA calculation is active. See data_indicator::osa_active().