C++ Review

class: center, middle, inverse, title-slide

# C++ Review
## TMB Training Session II
### Andrea Havron<br>NOAA Fisheries, OST

---

layout: true

.footnote[U.S. Department of Commerce | National Oceanic and Atmospheric Administration | National Marine Fisheries Service]

code.cpp{
  font-size: 14px;
}
code.r{
  font-size: 14px;
}

</style>

---
# C++ and Memory

.column-60[
C++ code

```cpp
#include <iostream>

int main() {
  //initiate a variable
  float a = 3.1459;
  //initiate a new variable with the same value as a
  float b = a;
  //initiate a variable c that points to the same address as a
  float* c = &a;
  std::cout << "c is equal to the address of a; *c = a" << std::endl;
  std::cout << "c = " << c << std::endl;
  std::cout << "a = " << a << "; *c = " << *c << std::endl;
  return 0;
}
```
]

.column-40-left[
Console output

```r
g++ memoryexample.cpp -o a.exe
a.exe

c is equal to the address of a; *c = a
c = 0x78fe14
a = 3.1459; *c = 3.1459
```
]

---
# C++ and Memory

.column-60[
C++ code

```cpp
#include <iostream>

int main() {
  //initiate a variable
  float a = 3.1459;
  //initiate a new variable with the same value as a
  float b = a;
  //initiate a variable c that points to the same address as a
  float* c = &a;
  *c = 100;
  std::cout << "a and *c have been updated; b has not"  << std::endl;
  std::cout << "a = " << a << "; *c = " << *c << std::endl;
  std::cout << "b = " << b << std::endl;
  return 0;
}
```
]

.column-40-left[
Console output

```r
g++ memoryexample.cpp -o a.exe
a.exe

a and *c have been updated; b has not
a = 100; *c = 100
b = 3.1459
```
]
---
# C++ and Memory

.column-60[
C++ code

```cpp
#include <iostream>

int main() {
  //initiate a variable
  float a = 3.1459;
  //initiate a new variable with the same value as a
  float b = a;
  //initiate a variable c that points to the same address as a
  float* c = &a;
  *c = 100;
  c = &b;
  b = 10;
  a = 3;
  std::cout << "c now equals the address of b" << std::endl;
  std::cout << "c = " << c << "; &b = " << &b << std::endl;
  std::cout << "a = " << a << std::endl;
  std::cout << "b = " << b << std::endl; 
  std::cout << "*c = " << *c ;
  return 0;
}
```
]

.column-40-left[
Console output

```r
g++ memoryexample.cpp -o a.exe
a.exe

c now equals the address of b
c = 0x78fe10; &b = 0x78fe10
a = 3
b = 10
 *c = 10
```
]
---
# Update values using pointers

.column-60[
C++ code

```cpp
#include <iostream>
#include <cmath>

template <class T>
  T my_add(T *x, T y){
    *x += y;
    return 0;
}
template <class T>
  T my_exp(T *x){
    *x = std::exp(*x);
    return 0;
}
int main() {
  //initiate a variable
  float par = 0;
  //initiate a new constant used to update par
  float b = 0.6931472;
  //intiate a pointer to par
  float *c = &par;
  //update par using functions
  my_add(c, b);
  std::cout << "par = " << par << std::endl;
  my_exp(c);
  std::cout << "par = " << par << std::endl;
  return 0;
}
```
]

.column-40-left[
Console output

```r
g++ pointerexample.cpp -o a.exe
a.exe

par = 0.6931472
par = 2
```
]

---
# C++ Anatomy of a Class

```cpp
class ClassName{
  Access specifier: //can be private, public, or protected
  Data members; //variables to be used
  
  ClassName(){} //Constructor - automatically called when new object created
  
  Member functions(){} //Methods to access data members
  
}; //class name ends with semicolon

//create instance of class
ClassName<T>* ClassName<T>::instance = new ClassName<T>();

```
- Memory is not allocated when a class is defined
- An instance of the class needs to be created for memory allocation (new object created)
- **new** is the C++ initialize function that allocates memory (creates pointer)
- Once the class is instantiated, members can be accessed using
 - the dot(**.**) when class is declared
 - the pointer(**->**) when an instance of the class is created using **new**

resources:<br>
[**geeksforgeeks**](https://www.geeksforgeeks.org/c-classes-and-objects/)<br>
[**cplusplus**](https://cplusplus.com/doc/tutorial/classes/)
---
# Class Example

C++ code
.pull-left[

```cpp
#include <iostream>

template <class T>
class MyClass{
    public:
    T a;
    T b;
    MyClass(){}
    T evaluate(){
        return a + b;
    }
};
int main() {
  //initiate class with pointer
  MyClass<double> *m = new MyClass<double>();
  m -> a = 1.2;
  m -> b = 3.4;
  std::cout << m -> evaluate() << std::endl;  
  //intiate class with declaration
  MyClass<double> myclass;
  myclass.a = 2;
  myclass.b = 3;
  std::cout << myclass.evaluate() << std::endl;
  return 0;
}
```
]
.pull-right[
Console output

```r
g++ classexample.cpp -o a.exe
a.exe

4.6
5
```
]
---
# Modular TMB structure: .hpp files
.pull-left[
[interface.hpp](https://github.com/NOAA-FIMS/FIMS/blob/main/inst/include/interface/interface.hpp) - setup TMB dependencies and define data types

```cpp
#ifdef TMB_MODEL
#include <TMB.hpp>

template<typename T>
struct ModelTraits{
 typedef typename CppAD::vector<Type> DataVector;
 typedef typename CppAD::vector<T> ParameterVector;
 typedef typename tmbutils::vector<T> EigenVector;
};

#endif /* TMB_MODEL */
```
]
.pull-right[
[test_dnorm_distribution.hpp](https://github.com/NOAA-FIMS/FIMS/blob/main/inst/extdata/TMB_tests/distributions/test_dnorm_distribution.hpp) - model

```cpp
#include "interface.hpp"
#include "distributions.hpp"

template<typename T>
class Model{
  
  using DataVector = typename ModelTraits<T>::DataVector;
      
  public:
  DataVector y; /*!< observation */
  T mean;  /*!< mean of the normal distribution */
  T sd; /*!< standard deviation of the normal distribution, must be strictly positive.*/

// Initiate pointer to link .cpp to .hpp
  static Model<T>* instance;

/** @brief Constructor.*/
  Model(){}

// @brief Create new singleton class
  // @return Model<T>* 
  static Model<T>* getInstance(){
    return Model<T>::instance;
  }
  ...
```
  ]
---
# Modular TMB structure: .hpp files
.pull-left[
[test_dnorm_distribution.hpp](https://github.com/NOAA-FIMS/FIMS/blob/main/inst/extdata/TMB_tests/distributions/test_dnorm_distribution.hpp) - model

```cpp
#include "interface.hpp"
#include "distributions.hpp"
namespace fims {
template<typename T>
class Model{
  
  using DataVector = typename ModelTraits<T>::DataVector;
      
  public:
  DataVector y; /*!< observation */
  T mean;  /*!< mean of the normal distribution */
  T sd; /*!< standard deviation of the normal distribution, must be strictly positive.*/

// Initiate pointer to link .cpp to .hpp
  static Model<T>* instance;

/** @brief Constructor.*/
  Model(){}

// @brief Create new singleton class
  // @return Model<T>* 
  static Model<T>* getInstance(){
    return Model<T>::instance;
  }
  ...
```
]
.pull-right[
<br>

```cpp
  // @brief Function that calculates the negative log-likelihood given the data and parameters
  // @return negative log-likelihood (nll)
  T evaluate(){
    T nll = 0;
    int i;
    int n = y.size();
    fims::Dnorm<T> nll_dnorm;
    nll_dnorm.mean = mean;
    nll_dnorm.sd = sd;
    for(i =0; i < n; i++){
      nll_dnorm.x = y[i];
      nll -= nll_dnorm.evaluate(true);
    }
    return nll;
  }
};

// @brief Create new instance of Model
  // @tparam T 
  template<class T>
  Model<T>* Model<T>::instance = new Model<T>();
}}
```
]
---
# Modular TMB structure: .cpp file
[test_dnorm_distribution.cpp](https://github.com/NOAA-FIMS/FIMS/blob/main/inst/extdata/TMB_tests/distributions/test_dnorm_distribution.cpp) - model

```cpp
#include "test_dnorm_distribution.hpp"

template<class Type>
Type objective_function<Type>::operator()(){
  // create pointer, inst, that points to singleton class of Model in test_dnorm_distribution.hpp
  //getinstance is defined in test_dnorm_distribution.hpp
  fims::Model<Type>* inst = fims::Model<Type>::getInstance();

DATA_VECTOR(y);
  PARAMETER_VECTOR(p);
  Type mean = p[0];
  Type sd = exp(p[1]);

// access and assign members of Model class using inst pointer
  inst -> y = y;
  inst -> mean = mean;
  inst -> sd = sd;

// create Type nll and assign value to the return of the 
  // evaluate() function defined in test_dnorm_distribution.hpp
  Type nll = inst -> evaluate();
  return nll;
}
```

---
# Modular TMB structure: R unit test
From [test_distributions.R](https://github.com/NOAA-FIMS/FIMS/blob/main/inst/extdata/TMB_tests/distributions/test_distributions.R)

```r
library(testthat)
library(TMB)

# Compile code with TMB_MODEL defined
TMB::compile(paste0(path, "/test_dnorm_distribution.cpp"), flags = "-DTMB_MODEL")

# dnorm unit test
# load test
dyn.load(dynlib(paste0(path, "/test_dnorm_distribution")))

#Simulate new data with R
set.seed(123)
y = rnorm(10, 5, 3)

# Calculate negative log-likelihood with R dnorm
nll = -sum(stats::dnorm(y, 5,3, TRUE))

# Initialize TMB model object with true values
mod = MakeADFun(data = list(y =y),
                parameters = list(p = c(5, log(3))),
                DLL = "test_dnorm_distribution")

# Compare R nll to TMB nll
expect_equal(nll, mod$fn())
```
---
# Rcpp Interface
.pull-left[
[Matthew's example](https://github.com/NOAA-FIMS/ModularTMBExample) using a vonBertalanffy model

```cpp

class Variable{
public:
    static std::vector<Variable*> parameters;
    bool estimable = FALSE;
    double value = 0;
    
    Variable(){
      Variable::parameters.push_back(this);
    }
    
};

std::vector<Variable*> Variable::parameters;

// Returns the initial values for the parameter set
Rcpp::NumericVector get_parameter_vector(){
  Rcpp::NumericVector p;
    
  for(int i = 0; i < Variable::parameters.size(); i++){
    if(Variable::parameters[i]->estimable){
        p.push_back(Variable::parameters[i]->value);
    }
  }
  return p;
}
```
]

.pull-right[
TMB model, .cpp file

```cpp
template<typename Type>
Type objective_function<Type>::operator()(){
    
    //get the singleton instance for type Type
    vonBertalanffyInterface::instance->prepare_template<Type>();
    VonBertalanffyModel<Type>* model =
        VonBertalanffyModel<Type>::getInstance();
    
    //get the parameter values
    PARAMETER_VECTOR(p)
    
    //update the parameter values for type Type
    for(int i =0; i < model->parameters.size(); i++){
        *model->parameters[i] = p[i];
    }
    
    //evaluate the model objective function value
    return model->evaluate();
}
```
]
---
#Rcpp Interface, calling from R
.pull-left[

```r
#get the Rcpp module
g <- Rcpp::Module(module = "growth",
                  PACKAGE = "ModularTMBExample")
vonB<-new(g$vonBertalanffy)

#initialize k
vonB$k$value<-.05
vonB$k$estimable<-TRUE

#initialize a_min
vonB$a_min$value<-.01
vonB$a_min$estimable<-FALSE

#initialize l_inf
vonB$l_inf$value<-7
vonB$l_inf$estimable<-TRUE

#set data
vonB$data <-data

#set ages 
vonB$ages<-ages

#prepare for interfacing with TMB
vonB$prepare()
```
]

.pull-right[

```r
#create an empty data list (data set above)
data <- list()

#create a parameter list
parameters <- list(
  p = g$get_parameter_vector()
)

#Run TMB
obj <- MakeADFun(data, parameters, DLL="ModularTMBExample")
opt <- nlminb(obj$par, obj$fn, obj$gr)
```
]
---
# TMB and FIMS
[**fimsflow**](https://docs.google.com/presentation/d/1U8N6J6YQIdd0Wer-o5vZvGTHgw1GLA8_6a0r8sYfqXo/edit#slide=id.g13fa6bc6e99_0_0)
<img src="data:image/png;base64,#static/fims-tmb1.png" width="60%" style="display: block; margin: auto auto auto 0;" />
---
# TMB and FIMS
[**fimsflow**](https://docs.google.com/presentation/d/1U8N6J6YQIdd0Wer-o5vZvGTHgw1GLA8_6a0r8sYfqXo/edit#slide=id.g1405a9bb7f8_0_0)
<img src="data:image/png;base64,#static/fims-tmb2.png" width="80%" style="display: block; margin: auto auto auto 0;" />