---
title: "Introducing the Fisheries Integrated Modeling System (FIMS)"
output: github_document
vignette: >
  %\VignetteIndexEntry{Introducing the Fisheries Integrated Modeling System (FIMS)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
# library(dplyr)
```

## Fisheries Integrated Modeling System
The NOAA Fisheries Integrated Modeling System (FIMS) is a new modeling framework for fisheries modeling. FIMS is a software system designed and architected to support next-generation fisheries stock assessment, ecosystem, and socioeconomic modeling. It's important to note that FIMS itself is not a model, but rather a framework for creating models. The framework is made up of many modules that come together to create a "the best model" that suites the needs of the end-user. What follows is a demo of creating a catch-at-age assessment model using FIMS. 

## Creating Models in FIMS
To begin, we import the FIMS and TMB libraries. Calling `library(FIMS)` automatically loads the Rcpp functions and modules into the R environment. The function call, `clear()`, ensures C++ memory from any previous fims model run is cleared out. 
```{r fims1, warning=FALSE, message=FALSE}
# automatically loads fims Rcpp module
library(FIMS)
library(TMB)

# clear memory
clear()
```

## Setting up Data
Data and variable values are taken from the [Li et. al.](https://spo.nmfs.noaa.gov/content/fishery-bulletin/comparison-4-primary-age-structured-stock-assessment-models-used-united) Model Comparison project ([github site](https://github.com/Bai-Li-NOAA/Age_Structured_Stock_Assessment_Model_Comparison)). See [R/data_mile1.R](https://github.com/NOAA-FIMS/FIMS/blob/main/R/data_mile1.R) and [tests/testthat/test-fims-estimation.R](https://github.com/NOAA-FIMS/FIMS/blob/integration-test-estimation-R-clear-FIMS/tests/testthat/test-fims-estimation.R) for details on how data and variable values are read into FIMS from the Model Comparison project.

First let's set up the dimensions of the model based on the Model Comparison project:
```{r fims-dims}
nyears <- 30 # the number of years which we have data for.
nseasons <- 1 # the number of seasons in each year. FIMS currently defaults to 1
ages <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12) # age vector.
nages <- 12 # the number of age groups.
```

### Preparing Data using FIMSFrame

 We will be reading data into the model using the FIMSFrame S4 R class set up in [R/fimsframe.R](https://github.com/NOAA-FIMS/FIMS/blob/main/R/fimsframe.R)

```{r fimsframe}
# use FIMS data frame
data(package = "FIMS")
fims_frame <- FIMSFrame(data_mile1)
```

The `fims_frame` object contains a `@data` slot that holds a long data frame with catch data for the fishery and index data for the survey:

```{r ageframe}
str(fims_frame)
fims_frame@data |>
  dplyr::filter(type == "landings") |>
  utils::head()
fims_frame@data |>
  dplyr::filter(type == "index") |>
  utils::head()
```


Using this data frame, we will start setting up the FIMS data objects. This example from the Model Comparison project sets up a single fishery fleet with age composition and catch data and a single survey with age composition data and an index. Data are read into FIMS as long vectors, regardless of their original dimension, hence the motivation behind the long data frames created with the fimsframe S4 classes. 

```{r data}
# fishery data
fishery_catch <- FIMS::m_landings(fims_frame)
fishery_agecomp <- FIMS::m_agecomp(fims_frame, "fleet1")

# survey data
survey_index <- FIMS::m_index(fims_frame, "survey1")

# survey agecomp not set up in fimsframe yet
survey_agecomp <- FIMS::m_agecomp(fims_frame, "survey1")
```

## Creating Modules in FIMS
Now that we've prepared the data, let's pass it into FIMS. Each module in the FIMS-R interface is made of S4 objects. These S4 objects serve as a interface between R and the underlining C++ code that defines FIMS. Modules are instantiated using the `methods::new()` function. We can use `methods::show()` to view all the fields (i.e. variables) and methods (i.e. functions) available in a given module.

### The Fleet Module

#### Fleet Data

Each fleet is required to have data in order to evaluate the objective function. Currently FIMS only has a fleet module that is used to set up both fleets and surveys. FIMS contains an Index module and AgeComp module to pass data objects into the fleet module. Each of these data modules require a dimension be added to indicate the dimensions of the raw data (e.g. nyears x nages matrix). Any years with missing data should be specified with a value set to -999. Given this information, FIMS is able to correctly apply dimension folding for model output.

Using the `methods::show()` function, we can see that the Index module has a vector field named *index_data* and the AgeComp module has a vector field names *age_comp_data*. 

```{r fleet-show}
show(Index)
show(AgeComp)
```

We'll create both index and age composition modules for the fleet using the `methods::new()` function and pass in the data defined above from the Model Comparison project.
```{r fleet-set-data}
# fleet index data
fishing_fleet_index <- methods::new(Index, nyears)
# fleet age composition data
fishing_fleet_age_comp <- methods::new(AgeComp, nyears, nages)
fishing_fleet_index$index_data <- fishery_catch # unit: mt
# Effective sampling size is 200
fishing_fleet_age_comp$age_comp_data <- fishery_agecomp * 200 # unit: number at age; proportion at age also works
```

#### Fleet Selectivity
Now that we've passed in data for the fishing fleet, we need to set up its selectivity module. We will set this to be selectivity function using the LogisticSelectivity module. The`methods::show()` function indicates this module has two parameter fields: *inflection_point* and *slope*, and an `evaluate()` and `get_id()` function. 

Each variable of [Parameter class](https://github.com/NOAA-FIMS/FIMS/blob/main/inst/include/interface/rcpp/rcpp_objects/rcpp_interface_base.hpp) has three additional fields: *value*, *is_random_effect*, and *estimated*. Currently, FIMS is not set up to run random effects. The default value for this field and the *estimate* field is currently set to `FALSE`. We can use the *value* field to input variables defined in the Model Comparison project.
 
```{r fleet_selectivity}
methods::show(LogisticSelectivity)
fishing_fleet_selectivity <- methods::new(LogisticSelectivity)
fishing_fleet_selectivity$inflection_point$value <- 2.0
fishing_fleet_selectivity$inflection_point$is_random_effect <- FALSE
fishing_fleet_selectivity$inflection_point$estimated <- TRUE
fishing_fleet_selectivity$slope$value <- 1.0
fishing_fleet_selectivity$slope$is_random_effect <- FALSE
fishing_fleet_selectivity$slope$estimated <- TRUE
```

#### Creating the Fleet Object
Now that we've created everything that a fleet needs, lets create the actual fleet object. First let's run `methods::show(Fleet)` to see all the fields and methods available from R. 

```{r show-Fleet}
show(Fleet)
```
We can see that there are five boolean flags: estimate_F, estimate_q, and is_survey, random_F, and random_q. There are two vectors, log_Fmort and log_obs_error, and a double, log_q. There are two integer fields for the number of ages and years. Additionally, there are five Methods: SetAgeCompLikelihood, SetIndexLikelihood, SetObservedAgeCompData, SetObservedIndexData, and setSelectivity. The last three of these will be used to link up the AgeComp, Index, and Selectivity modules defined above with the fleet module defined below. 


```{r fleet}
# Create fleet module
fishing_fleet <- methods::new(Fleet)
# Set nyears and nages
fishing_fleet$nages <- nages
fishing_fleet$nyears <- nyears
# Set values for log_Fmort
fishing_fleet$log_Fmort <- log(c(
  0.009459165, 0.02728886, 0.04506364,
  0.06101782, 0.04860075, 0.08742055,
  0.0884472, 0.1866079, 0.109009, 0.1327043,
  0.1506155, 0.161243, 0.1166402, 0.1693461,
  0.1801919, 0.1612405, 0.3145732, 0.2572476,
  0.2548873, 0.2514621, 0.3491014, 0.2541077,
  0.4184781, 0.3457212, 0.3436855, 0.3141712,
  0.3080268, 0.4317453, 0.3280309, 0.4996754
))
# Turn on estimation for F
fishing_fleet$estimate_F <- TRUE
fishing_fleet$random_F <- FALSE
# Set value for log_q
fishing_fleet$log_q <- log(1.0)
fishing_fleet$estimate_q <- FALSE
fishing_fleet$random_q <- FALSE
fishing_fleet$log_obs_error <- rep(log(sqrt(log(0.01^2 + 1))), nyears)
fishing_fleet$estimate_obs_error <- FALSE
# Next two lines not currently used by FIMS
fishing_fleet$SetAgeCompLikelihood(1)
fishing_fleet$SetIndexLikelihood(1)
# Set Index, AgeComp, and Selectivity using the IDs from the modules defined above
fishing_fleet$SetObservedIndexData(fishing_fleet_index$get_id())
fishing_fleet$SetObservedAgeCompData(fishing_fleet_age_comp$get_id())
fishing_fleet$SetSelectivity(fishing_fleet_selectivity$get_id())
```

### The Survey Module
We will now repeat the steps from Fleet to set up the Survey. A survey object is essentially the same as a fleet object with a catchability (q) variable.

#### Survey Data

```{r survey-set-data}
# fleet index data
survey_fleet_index <- methods::new(Index, nyears)
# survey age composition data
survey_fleet_age_comp <- methods::new(AgeComp, nyears, nages)
survey_fleet_index$index_data <- survey_index # unit: mt; it's possible to use other units as long as the index is assumed to be proportional to biomass
# Effective sampling size is 200
survey_fleet_age_comp$age_comp_data <- survey_agecomp * 200 # unit: number at age; proportion at age also works
```

#### Survey Selectivity

```{r survey-selectivity}
survey_fleet_selectivity <- new(LogisticSelectivity)
survey_fleet_selectivity$inflection_point$value <- 1.5
survey_fleet_selectivity$inflection_point$is_random_effect <- FALSE
survey_fleet_selectivity$inflection_point$estimated <- TRUE
survey_fleet_selectivity$slope$value <- 2.0
survey_fleet_selectivity$slope$is_random_effect <- FALSE
survey_fleet_selectivity$slope$estimated <- TRUE
```

#### Creating the Survey Object

```{r survey}
survey_fleet <- methods::new(Fleet)
survey_fleet$is_survey <- TRUE
survey_fleet$nages <- nages
survey_fleet$nyears <- nyears
survey_fleet$estimate_F <- FALSE
survey_fleet$random_F <- FALSE
survey_fleet$log_q <- log(3.315143e-07)
survey_fleet$estimate_q <- TRUE
survey_fleet$random_q <- FALSE
# sd = sqrt(log(cv^2 + 1)), sd is log transformed
survey_fleet$log_obs_error <- rep(log(sqrt(log(0.2^2 + 1))), nyears)
survey_fleet$estimate_obs_error <- FALSE
survey_fleet$SetAgeCompLikelihood(1)
survey_fleet$SetIndexLikelihood(1)
survey_fleet$SetSelectivity(survey_fleet_selectivity$get_id())
survey_fleet$SetObservedIndexData(survey_fleet_index$get_id())
survey_fleet$SetObservedAgeCompData(survey_fleet_age_comp$get_id())
```  


### Creating a Population
The final step is to set up the population module. Before doing so, we first need to set up each component of the population (e.g. recruitment, growth, etc.). 

#### Recruitment
We'll use a Beverton Holt recruitment module. We first instantiate a module using the `methods::new()` function. We can use `methods::show()` to view all the fields (i.e. variables) and methods (i.e. functions) available in `BevertonHoltRecruitment` module.

```{r recruitment}
# Recruitment
recruitment <- methods::new(BevertonHoltRecruitment)
methods::show(BevertonHoltRecruitment)
```
There are three parameters we need to set-up: *log_sigma_recruit*, *log_rzero*, and *logit_steep*.

```{r set-up-recruitment}
recruitment$log_sigma_recruit$value <- log(0.4)
recruitment$log_rzero$value <- log(1e+06) # unit: log(number)
recruitment$log_rzero$is_random_effect <- FALSE
recruitment$log_rzero$estimated <- TRUE
recruitment$logit_steep$value <- -log(1.0 - 0.75) + log(0.75 - 0.2)
recruitment$logit_steep$is_random_effect <- FALSE
recruitment$logit_steep$estimated <- FALSE
```

We also need to set up log recruitment deviations. FIMS recruitment modules have a boolean, *estimate_log_devs* to specify whether or not log deviations are estimated; and a vector, *log_devs* to set the log deviation values. 

```{r rec-devs}
recruitment$estimate_log_devs <- FALSE
recruitment$log_devs <- c(
  0.08904850, 0.43787763, -0.13299042, -0.43251973,
  0.64861200, 0.50640852, -0.06958319, 0.30246260,
  -0.08257384, 0.20740372, 0.15289604, -0.21709207,
  -0.13320626, 0.11225374, -0.10650836, 0.26877132,
  0.24094126, -0.54480751, -0.23680557, -0.58483386,
  0.30122785, 0.21930545, -0.22281699, -0.51358369,
  0.15740234, -0.53988240, -0.19556523, 0.20094360,
  0.37248740, -0.07163145
)
```
#### Growth
Now, we'll define the growth module for our population using an empirical weight at age model.
```{r growth}
# Growth
ewaa_growth <- methods::new(EWAAgrowth)
ewaa_growth$ages <- ages
ewaa_growth$weights <- c(
  0.0005306555, 0.0011963283, 0.0020582654,
  0.0030349873, 0.0040552124, 0.0050646975,
  0.0060262262, 0.0069169206, 0.0077248909,
  0.0084461128, 0.0090818532, 0.0096366950
) # unit: mt
```   
#### Maturity
Each population will also need a maturity model. Here we define a logistic maturity model. 
```{r maturity}
# Maturity
maturity <- new(LogisticMaturity)
maturity$inflection_point$value <- 2.25
maturity$inflection_point$is_random_effect <- FALSE
maturity$inflection_point$estimated <- FALSE
maturity$slope$value <- 3
maturity$slope$is_random_effect <- FALSE
maturity$slope$estimated <- FALSE
```

Now that our life history sub-models are defined, lets define the actual population.

```{r population}
# Population
population <- new(Population)
population$log_M <- rep(log(0.2), nyears * nages)
population$estimate_M <- FALSE
population$log_init_naa <- log(c(
  993947.5, 811707.8, 661434.4,
  537804.8, 436664.0, 354303.4,
  287397.0, 233100.2, 189054.0,
  153328.4, 124353.2, 533681.3
)) # unit: in number
population$estimate_init_naa <- TRUE
population$nages <- nages
population$ages <- ages
population$nfleets <- 2 # 1 fleet and 1 survey
population$nseasons <- nseasons
population$nyears <- nyears
```

Now we need to link up the recruitment, growth, and maturity modules we set above with this new population module. We do this by calling `get_id()` from each respective module and passing that unique ID into each respective `Set` function from population. 

```{r set-pop-modules}
population$SetMaturity(maturity$get_id())
population$SetGrowth(ewaa_growth$get_id())
population$SetRecruitment(recruitment$get_id())
```


## Putting It All Together 

### Creating the FIMS Model and Making the TMB Function

```{r model}
sucess <- CreateTMBModel()
parameters <- list(p = get_fixed())
obj <- MakeADFun(data = list(), parameters, DLL = "FIMS", silent = TRUE)
```
## Fitting the Model
```{r fit_model}
opt <- nlminb(obj$par, obj$fn, obj$gr,
  control = list(eval.max = 800, iter.max = 800)
) # , method = "BFGS",
#   control = list(maxit=1000000, reltol = 1e-15))

print(opt)
```

### TMB Reporting
```{r tmb_report}
sdr <- TMB::sdreport(obj)
sdr_fixed <- summary(sdr, "fixed")
report <- obj$report(obj$env$last.par.best)

print(sdr_fixed)

# report out nll components
report$rec_nll # recruitment
report$index_nll # fishery catch and survey index
report$age_comp_nll # fishery and survey age composition
```


## Plotting Results
```{r plots}
library(ggplot2)
index_results <- data.frame(
  observed = survey_fleet_index$index_data,
  expected = report$exp_index[[2]]
)
print(index_results)

ggplot(index_results, aes(x = 1:nyears, y = observed)) +
  geom_point() +
  xlab("Year") +
  ylab("Index (mt)") +
  geom_line(aes(x = 1:nyears, y = expected), color = "blue") +
  theme_bw()

catch_results <- data.frame(
  observed = fishing_fleet_index$index_data,
  expected = report$exp_index[[1]]
)
print(catch_results)

ggplot(catch_results, aes(x = 1:nyears, y = observed)) +
  geom_point() +
  xlab("Year") +
  ylab("Index (mt)") +
  geom_line(aes(x = 1:nyears, y = expected), color = "blue") +
  theme_bw()
```

## Clear C++ objects from memory
```{r clear}
clear()
```