add random effects example - not working

inst/include/common/information.hpp

@@ -8,6 +8,7 @@
 
 #include "../distributions/density_components_base.hpp"
 #include "../distributions/normal_lpdf.hpp"
+#include "../distributions/ar1_lpdf.hpp"
 #include "../pop_dy/population.hpp"
 #include "../pop_dy/von_bertalanffy.hpp"
 
@@ -91,12 +92,21 @@ class Information {
           variable_map_iterator vmit;
           vmit = this->variable_map.find(n->key[0]); 
           n->observed_value = *(*vmit).second;
-
           for(size_t i=1; i<n->key.size(); i++){
             vmit = this->variable_map.find(n->key[i]); 
-            n->observed_value.insert(std::end(n->expected_value), 
+            n->observed_value.insert(std::end(n->observed_value), 
                                      std::begin(*(*vmit).second), std::end(*(*vmit).second));
           } 
+          std::shared_ptr<DensityComponentBase<Type> > density_components_base = density_components[n->id];
+          AR1LPDF<Type>* ar1 = (AR1LPDF<Type>*) density_components_base.get();
+          Rcout << "inside setup_re, obs.val size is: " << n->observed_value.size() << std::endl;
+          Rcout << "inside setup_re, rho is: " << ar1->rho[0] << std::endl;
+          for (pop_iterator it = this->pop_models.begin();
+               it != this->pop_models.end(); ++it) {
+               std::shared_ptr<Population<Type> > pop = (*it).second;
+            Rcout << "pop u size is: " << pop->u.size() << std::endl;
+            Rcout << "first pop u is: " << pop->u[0] << std::endl;
+          }
         }
       }
     }
@@ -116,6 +126,21 @@ class Information {
       }
     }
   }
+
+    bool CreateModel(){
+      bool valid_model = true;
+      for (pop_iterator it = this->pop_models.begin(); it != this->pop_models.end();
+      ++it) {
+        setup_population();
+        setup_priors();
+        setup_random_effects();
+        Rcout << "setup_random_effects successful!" << std::endl;
+        setup_data();
+      }
+      return valid_model;
+    }
+
+
 };
 
 template<typename Type>

inst/include/common/model.hpp

@@ -63,11 +63,8 @@ class Model{
   Type evaluate(){
     Type jnll = 0.0;
 
-    this->info->setup_population();
-
     //maybe here, setup functions can take a simulate flag and simulation can be controlled from model
     //setup pointers for priors
-    this->info->setup_priors();
      for(density_components_iterator it = this->density_components.begin(); it!= this->density_components.end(); ++it){
       std::shared_ptr<DensityComponentBase<Type> > n = (*it).second;
       #ifdef TMB_MODEL
@@ -80,24 +77,24 @@ class Model{
 
 
     //setup pointers for random effects
-    info->setup_random_effects();
     //evaluate nlls for priors and random effects
      for(density_components_iterator it = this->density_components.begin(); it!= this->density_components.end(); ++it){
       std::shared_ptr<DensityComponentBase<Type> > n = (*it).second;
       #ifdef TMB_MODEL
         n->of = this->of;
       #endif
-      if(n->input_type == "random_effects"){
+        Rcout << "inside model, just before re evaluate" << std::endl;
+      if(n->input_type == "re"){
         jnll -= n->evaluate();
       }
+      Rcout << "inside model, just after re evaluate" << std::endl;
     }
 
     for (pop_iterator it = this->pop_models.begin();
          it != this->pop_models.end(); ++it) {
       std::shared_ptr<Population<Type> > pop = (*it).second;
       pop->evaluate(); 
     }
-    this->info->setup_data();
     for(density_components_iterator it = this->density_components.begin(); it!= this->density_components.end(); ++it){
       std::shared_ptr<DensityComponentBase<Type> > n = (*it).second;
       #ifdef TMB_MODEL

inst/include/distributions/ar1_lpdf.hpp

@@ -25,13 +25,16 @@ struct AR1LPDF : public DensityComponentBase<Type> {
     virtual ~AR1LPDF() {}
 
     virtual const Type evaluate(){
-        //stationarity assumption: -1>rho<1
-        //randomwalk assumption: rho = 1
-        if(rho[0] != 1){
+      Rcout << "AR1 observed size is: " << this->observed_value.size() << std::endl;
+      Rcout << "AR1 sd is: " << sd[0] << std::endl;
+      Rcout << "AR1 rho is: " << rho[0] << std::endl;
+      if(rho[0] != 1){
             rho[0] = 1 / (1 + exp(-logit_rho[0])) * 2 - 1;
         }
 
         sd[0] = exp(log_sd[0]);
+        Rcout << "AR1 sd is: " << sd[0] << std::endl;
+        Rcout << "AR1 rho is:" << rho[0] << std::endl;
         log_likelihood = -density::SCALE(density::AR1(rho[0]), sd[0])(this->observed_value);
         this->log_likelihood_vec[0] = log_likelihood;
 

inst/include/interface/rcpp/rcpp_interface.hpp

@@ -13,6 +13,27 @@ bool CreateModel(){
        i++) {
     RcppInterfaceBase::interface_objects[i]->prepare();
   }
+
+  // base model
+  std::shared_ptr<Information<TMB_FIMS_REAL_TYPE>> d0 =
+    Information<TMB_FIMS_REAL_TYPE>::getInstance();
+  d0->CreateModel();
+
+  // first-order derivative
+  std::shared_ptr<Information<TMB_FIMS_FIRST_ORDER>> d1 =
+    Information<TMB_FIMS_FIRST_ORDER>::getInstance();
+  d1->CreateModel();
+
+  // second-order derivative
+  std::shared_ptr<Information<TMB_FIMS_SECOND_ORDER>> d2 =
+    Information<TMB_FIMS_SECOND_ORDER>::getInstance();
+  d2->CreateModel();
+
+  // third-order derivative
+  std::shared_ptr<Information<TMB_FIMS_THIRD_ORDER>> d3 =
+    Information<TMB_FIMS_THIRD_ORDER>::getInstance();
+  d3->CreateModel();
+
   return true;
 }
 
@@ -147,6 +168,7 @@ RCPP_MODULE(growth) {
     Rcpp::class_<PopulationInterface>("Population")
     .constructor()
     .field("ages", &PopulationInterface::ages)
+    .field("u", &PopulationInterface::u)
     .method("get_id", &PopulationInterface::get_id)
     .method("set_growth", &PopulationInterface::SetGrowth)
     .method("get_module_name", &PopulationInterface::get_module_name);
@@ -179,6 +201,7 @@ RCPP_MODULE(growth) {
       .field("observed_value", &AR1LPDFInterface::observed_value)
       .field("log_sd", &AR1LPDFInterface::log_sd)
       .field("logit_rho", &AR1LPDFInterface::logit_rho)
+      .field("rho", &AR1LPDFInterface::rho)
       .field("input_type", &AR1LPDFInterface::input_type)
       .field("log_likelihood_vec", &AR1LPDFInterface::log_likelihood_vec)
       .field("simulate_flag", &AR1LPDFInterface::simulate_flag)

inst/include/interface/rcpp/rcpp_objects/rcpp_distributions.hpp

@@ -387,13 +387,15 @@ class AR1LPDFInterface : public DensityComponentsInterface{
         ar1->key.resize(this->key.size());
         for(int i=0; i<key.size(); i++){
             ar1->key[i] = this-> key[i];
+          Rcout << "ar1 key is: " << ar1->key[i] << std::endl;
         }
         ar1->simulate_flag = this->simulate_flag;
         ar1->osa_flag = false;
 
         ar1->observed_value.resize(this->observed_value.size());
         for(size_t i=0; i<this->observed_value.size(); i++){
             ar1->observed_value[i] = this->observed_value[i].value;
+          Rcout << "ar1 observed value is: " << ar1->observed_value[i] << std::endl;
             if(this->observed_value[i].is_random_effect){
                 model->random_effects.push_back(&(ar1)->observed_value[i]);
             }
@@ -403,6 +405,7 @@ class AR1LPDFInterface : public DensityComponentsInterface{
         ss << this->get_module_name() << "_" << this->id << "_log_sd";
         for(size_t i=0; i<this->log_sd.size(); i++){
             ar1->log_sd[i] = this->log_sd[i].value;
+          Rcout << "ar1 log_sd is: " << ar1->log_sd[i] << std::endl;
             if(this->log_sd[i].estimable){
                 model->parameters.push_back(&(ar1)->log_sd[i]);
                 model->pnames.push_back(ss.str());
@@ -414,6 +417,8 @@ class AR1LPDFInterface : public DensityComponentsInterface{
         ss << this->get_module_name() << "_" << this->id << "_logit_rho";
         for(size_t i=0; i<this->logit_rho.size(); i++){
             ar1->logit_rho[i] = this->logit_rho[i].value;
+
+          Rcout << "ar1 logit_rho is: " << ar1->logit_rho[i] << std::endl;
             if(this->logit_rho[i].estimable){
                 model->parameters.push_back(&(ar1)->logit_rho[i]);
                 model->pnames.push_back(ss.str());
@@ -425,6 +430,7 @@ class AR1LPDFInterface : public DensityComponentsInterface{
         ss << this->get_module_name() << "_" << this->id << "_rho";
         for(size_t i=0; i<this->rho.size(); i++){
             ar1->rho[i] = this->rho[i].value;
+          Rcout << "ar1 rho is: " << ar1->rho[i] << std::endl;
             if(this->rho[i].estimable){
                 model->parameters.push_back(&(ar1)->rho[i]);
                 model->pnames.push_back(ss.str());
@@ -434,6 +440,9 @@ class AR1LPDFInterface : public DensityComponentsInterface{
 
 
         model->density_components[ar1->id] = ar1;
+       // Rcout << "inside rcpp_distribution, model rho is: " << model->density_components[ar1->id]->rho[0] << std::endl;
+        Rcout << "inside rcpp_distribution, ar1 rho is: " << ar1 -> rho[0] << std::endl;
+
         info->density_components[ar1->id] = ar1;
         return true;
     }
@@ -458,14 +467,15 @@ class AR1LPDFInterface : public DensityComponentsInterface{
      * Update the model parameter values and finalize. Sets the parameter values and evaluates the
      * portable model once and transfers values back to the Rcpp interface.
      */
-    void finalize(Rcpp::NumericVector v) {
+    void finalize(Rcpp::NumericVector v, Rcpp::NumericVector re) {
 
         std::shared_ptr< Model<double> > model = Model<double>::getInstance();
         std::shared_ptr<DensityComponentBase<double> > density_components_base = model->density_components[this->id];
         AR1LPDF<double>* ar1 = (AR1LPDF<double>*) density_components_base.get();  
 
         for (int i = 0; i < v.size(); i++) {
             (*model->parameters[i]) = v[i];
+            (*model->random_effects[i]) = re[i];
         }
 
         double f = model->evaluate();

inst/include/interface/rcpp/rcpp_objects/rcpp_population.hpp

@@ -40,6 +40,7 @@ uint32_t PopulationInterfaceBase::id_g = 1;
 class PopulationInterface  : public PopulationInterfaceBase {
 public:
     Rcpp::NumericVector ages;
+    VariableVector u;
     uint32_t growth_id;        /**< id of the growth function*/
 
 
@@ -82,6 +83,29 @@ class PopulationInterface  : public PopulationInterfaceBase {
         ss << key << "_length";
         info->variable_map[ss.str()] = &(pop)->length; 
         ss.str("");
+
+
+        if(this->u.size() == 0){
+          pop->u.resize(ages.size());
+          for(int i=0; i<ages.size(); i++){
+            pop->u[i] = 0.0;
+          }
+        } else {
+          pop->u.resize(this->u.size());
+          for(int i =0; i < this->u.size(); i++){
+            pop->u[i] = this->u[i];
+            if(this->u[i].is_random_effect){
+              model->random_effects.push_back(&(pop)->u[i]);
+            }
+          }
+        }
+
+        ss << this->get_module_name() << "_" << this->id;
+        ss.str("");
+        ss << key << "_u";
+        info->variable_map[ss.str()] = &(pop)->u; 
+        ss.str("");
+
 
         info->pop_models[pop->id] = pop;
         model->pop_models[pop->id] = pop;

inst/include/pop_dy/population.hpp

@@ -13,17 +13,20 @@ struct Population : public ModelObject<Type> {
 
   fims::Vector<Type> length;
   fims::Vector<Type> ages;
+  fims::Vector<Type> u;
 
   std::shared_ptr< VonBertalanffy<Type> > vb;
 
   Population() { 
     this->id = Population::id_g++;
   }
-
+  
 
   inline void CalculateLength(){ 
+
     for(int i =0; i < ages.size(); i++){
-        length[i] = vb -> evaluate(ages[i]);
+      Rcout << "u, " << i << " is: " << u[i] << std::endl;
+        length[i] = vb -> evaluate(ages[i]) + u[i];
     }
   }
 

src/ModularTMBExample.cpp

@@ -33,7 +33,15 @@ Type objective_function<Type>::operator()(){
     DATA_VECTOR(y);
     DATA_VECTOR_INDICATOR(keep,y);
     //get the parameter values
-    PARAMETER_VECTOR(p)
+    PARAMETER_VECTOR(p);
+    PARAMETER_VECTOR(re);
+
+    Rcout << "inside cpp, model data size is: " << model->data.size() << std::endl;
+    Rcout << "inside cpp, y size is: " << y.size() << std::endl;
+    Rcout << "inside cpp, model parameter size is: " << model->parameters.size() << std::endl;
+    Rcout << "inside cpp, p size is: " << p.size() << std::endl;
+    Rcout << "inside cpp, model random effects size is: " << model->random_effects.size() << std::endl;
+    Rcout << "inside cpp, re size is: " << re.size() << std::endl;;
 
     //update the data values for type Type
     for(int i =0; i < model->data.size(); i++){
@@ -44,6 +52,9 @@ Type objective_function<Type>::operator()(){
     for(int i =0; i < model->parameters.size(); i++){
         *model->parameters[i] = p[i];
     }
+    for(int i =0; i < model->random_effects.size(); i++){
+      *model->random_effects[i] = re[i];
+    }
 
     model -> of = this;
     model -> keep = keep;