inst/scripts/04.snowcrab_fishery_model_turing.r

In jae0/snowcrab: Snow crab assessment suite using aegis* and associated projects and data streams.

---
title: "Logistic Fishery Model"
author: "Jae S. Choi"
toc: true
number-sections: true
highlight-style: pygments
editor:
  render-on-save: false
  markdown:
    references:
      location: document
format:
  html: 
    code-fold: true
    html-math-method: katex
    self-contained: true
    embed-resources: true
  pdf:
    pdf-engine: lualatex
  docx: default 
---
 
   
<!-- Preamble


This is a Markdown document ... To create HTML or PDF, etc, run: 

As it runs R and Julia, Quarto is probably the better tool to render. Alternatively, you can run everything directly and then render just the figures by turning off the R and Julia code chunks.

  make quarto FN=04.snowcrab_fishery_model_turing YR=2023 SOURCE=~/projects/bio.snowcrab/inst/markdown WK=~/bio.data/bio.snowcrab/assessments DOCEXTENSION=html # {via Quarto}

  # use Rmarkdown only if already run code chunks manually and all figures have been generated.
  # make rmarkdown FN=04.snowcrab_fishery_model_turing YR=2023 SOURCE=~/projects/bio.snowcrab/inst/markdown WK=~/bio.data/bio.snowcrab/assessments DOCTYPE=html_document DOCEXTENSION=html # {via Rmarkdown}

  # or a direct pdf construction through pandoc ... not tested yet
  # make pdf FN=04.snowcrab_fishery_model_turing  # {via pandoc}


Alter year and directories to reflect setup or copy Makefile and alter defaults to your needs.
  
 
-->
 

## Purpose

Prepare aggregated time-series data for inference using a Bayesian fishery model. Currently, the default model is the discrete logistic formulation. It is labelled, "logistic_historical" in the following. There are other variations possible. See the calling code. 

---

## Prepare data for discrete logistic models

This step requires [bio.snowcrab/inst/scripts/03.biomass_index_carstm.r](https://github.com/jae0/bio.snowcrab/inst/scripts/03.biomass_index_carstm.r)   to be completed. That step statistically models the survey samples in space and time and then aggregates across space for an estimate/index of abundance. 

Here we read in these results and prepare (format) them for timeseries analysis.
```{r}
#| eval: true
#| output: false
  require(aegis)
  year.assessment = 2023
  p = bio.snowcrab::load.environment( year.assessment=year.assessment )
  modeldir = p$modeldir
  snowcrab_filter_class = "fb"     # fishable biomass (including soft-shelled )  "m.mat" "f.mat" "imm"

  carstm_model_label = paste( "default", snowcrab_filter_class, sep="_" )
  fishery_model_label = "turing1"
 
  # save to : carstm_directory = file.path(modeldir, carstm_model_label)

  fishery_model_data_inputs( 
    year.assessment=year.assessment,   
    type="biomass_dynamics", 
    snowcrab_filter_class="fb",
    modeldir=modeldir, 
    carstm_model_label=carstm_model_label,
    for_julia=TRUE,
    fishery_model_label="turing1"  
  ) 

  # the name/location of the above creates a data file (which will be passed to Julia, later) 
  fndat  = file.path( modeldir, carstm_model_label, "biodyn_biomass.RData" )  # 

  # for development using dynamical_model/snowcrab, save/copy at alternate locations
  # modeldir = file.path( homedir, "projects", "dynamical_model", "snowcrab", "data" ) 

```


---

## Model fitting and parameter inference

Now that data inputs are ready, we can continue on to model fitting and inference. This is done in (Julia)[https://julialang.org/]. 

First install or download Julia and install appropriate libraries from withing Julia. 

There are two ways you can use this. 

  1. Run within Julia directly.

    This method is most flexible and [documented here](https://github.com/jae0/dynamical_model/blob/master/snowcrab/04.snowcrab_fishery_model.md), where more options are also shown.  But it will require learning the language Julia and Turing library, but is very simple.

```{julia}
#| eval: false 
#| output: false

  year_assessment =2023
  yrs = 1999:year_assessment

  rootdir = homedir()
 
  project_directory  = joinpath( rootdir, "projects", "bio.snowcrab", "inst", "julia" ) 
  bio_data_directory = joinpath( rootdir, "bio.data", "bio.snowcrab", "modelled", "default_fb" )
  outputs_directory  = joinpath( rootdir, "bio.data", "bio.snowcrab", "fishery_model" ) 
 
  model_variation = "logistic_discrete_historical"   
  model_outdir = joinpath( outputs_directory, string(year_assessment), model_variation )
  
  include( joinpath(project_directory, "startup.jl" ) )  # load libs
  include( joinpath(project_directory, "logistic_discrete_functions.jl" ) )  # load core modelling functions
 
  o = load( joinpath( bio_data_directory, "biodyn_biomass.RData" ), convert=true)   # load data; alter file path as required   
 
  # run the whole script or in parts inside "bio.snowcrab/inst/julia/logistic_discrete.jl" for more control
  aulab ="cfanorth"     
  include( joinpath(project_directory, "logistic_discrete.jl" ) )   

  aulab ="cfasouth"    
  include( joinpath(project_directory, "logistic_discrete.jl" ) )   

  aulab ="cfa4x"      
  include( joinpath(project_directory, "logistic_discrete.jl" ) )   


```



<!--  **OR,**

  2. Run within R and call Julia indirectly, using the JuliaCall R-library (install that too). 

As this and the other Julia files exist in a Git repository and we do not want to add more files there, we copy the necessary files to a temporary work location such that new output files (figures, html documents, Julia logs, etc.) can be created there. 


```{r}
#| eval: true
#| output: false

  # Define location of data and outputs
  outputs_directory = file.path(data_root, "bio.snowcrab", "fishery_model" )
  model_outdir = file.path( outputs_directory, year.assessment, "logistic_discrete_historical" )
  work_dir = file.path( work_root, "fishery_model" )
  
  dir.create( outputs_directory, recursive=TRUE )
  dir.create( model_outdir, recursive=TRUE )
  dir.create( work_dir, recursive=TRUE)

  # copy julia scripts to a temp location where julia can write to (R-library is usually write protected)
  julia_scripts_location = system.file( "julia", package="bio.snowcrab" )
  file.copy( file.path(julia_scripts_location, "snowcrab_startup.jl"), work_dir, overwrite=TRUE  )
  file.copy( file.path(julia_scripts_location, "logistic_discrete_functions.jl"), work_dir, overwrite=TRUE )
  file.copy( file.path(julia_scripts_location, "logistic_discrete.jl"), work_dir, overwrite=TRUE )
 
      # to set up:
      # julia_setup(JULIA_HOME = "the folder that contains julia binary")
      # options(JULIA_HOME = "the folder that contains julia binary")
      # Set JULIA_HOME in command line environment.

 
  if ( !any( grepl("JuliaCall", o[,"Package"] ) ) ) install.packages("JuliaCall")

  # load JuliaCall interface
  library(JuliaCall)

  julia = try( julia_setup( install=FALSE, installJulia=FALSE ) )

  if ( inherits(julia, "try-error") ) {
    install_julia()
    julia = try( julia_setup( install=FALSE, installJulia=FALSE ) ) # make sure it works
    if ( inherits(julia, "try-error") )  stop( "Julia install failed, install manually?") 
  }

  ## transfer params to julia environment
  julia_assign( "year_assessment", year.assessment )  # copy data into julia session
  julia_assign( "bio_data_directory", data_root) 
  julia_assign( "outputs_directory", outputs_directory )  # copy data into julia session
  julia_assign( "model_outdir", model_outdir )  # copy data into julia session


  # julia_source cannot traverse directories .. temporarily switch directory

  currentwd = getwd() 
   
  setwd(work_dir) 

    # load/install libraries and setup directories 
    # .. if not all libraries are installed, you might need to re-run this a few times 
    # .. until there are no longer any pre-compilation messages
    julia_source( "snowcrab_startup.jl" )  
    
    # load core modelling functions
    julia_source( "logistic_discrete_functions.jl" )  

    # load data; alter file path as required
    julia_assign( "fndat", fndat )  # load data into julia session
    julia_command("o = load( fndat, convert=true) ")
 
    # Model and compute predictions for each area
    julia_command( "Random.seed!(1234)" )

    # modelling and outputs for each region 
    julia_assign( "yrs", p$yrs )   

    # cfanorth:    
    julia_assign( "aulab", "cfanorth" )   
    julia_source( "logistic_discrete.jl" )   

    # cfasouth:    
    julia_assign( "aulab", "cfasouth")  # copy data into julia session
    julia_source( "logistic_discrete.jl" )   # modelling and outputs

    # cfa4x:    
    julia_assign( "aulab", "cfa4x" )  # copy data into julia session
    julia_source( "logistic_discrete.jl" )   # modelling and outputs

  setwd(currentwd) # revert work directory

  # Example:
  # to move data into R
  # objects that might be useful:  m, num, bio, trace, Fkt, FR, FM 
  # bio = julia_eval("bio")  
```
 


-->


And that is it. Now we continue to look at the results: 



## Stock status: Biomass Index (aggregate) ... {.c}
    
\begin{tiny}
```{r fbindex-timeseries, out.width='60%', echo=FALSE, fig.align='center', fig.cap = 'The fishable biomass index (t) predicted by CARSTM of Snow Crab survey densities. Error bars represent Bayesian 95\\% Credible Intervals. Note large errors in 2020 when there was no survey.' }
include_graphics( file.path( SCD, 'modelled', 'default_fb', 'aggregated_biomass_timeseries' , 'biomass_M0.png') )
# \@ref(fig:fbindex-timeseries)
```
\end{tiny}




## Stock status: Modelled Biomass (pre-fishery) ... {.c}

\begin{tiny}
```{r logisticPredictions, out.width='32%', echo=FALSE, fig.show='hold', fig.align='center', fig.cap = 'Model 1 fishable, posterior mean modelled biomass (pre-fishery; kt) are shown in dark orange for N-ENS, S-ENS and 4X (left, middle and right). Light orange are posterior samples of modelled biomass (pre-fishery; kt) to illustrate the variability of the predictions. The biomass index (post-fishery, except prior to 2004) after model adjustment by the model catchability coefficient is in gray.' } 
loc = file.path( SCD, 'fishery_model', year.assessment, 'logistic_discrete_historical' )
fn1 = file.path( loc, 'plot_predictions_cfanorth.pdf' ) 
fn2 = file.path( loc, 'plot_predictions_cfasouth.pdf' ) 
fn3 = file.path( loc, 'plot_predictions_cfa4x.pdf' ) 
include_graphics(c(fn1, fn2, fn3) )
# \@ref(fig:logisticPredictions)
```
\end{tiny}



## Stock status: Fishing Mortality ... {.c}

```{r logisticFishingMortality, out.width='32%', echo=FALSE,  fig.show='hold', fig.align='center', fig.cap = 'Time-series of modelled instantaneous fishing mortality from Model 1, for N-ENS (left), S-ENS (middle), and 4X (right). Samples of the posterior densities are presented, with the darkest line being the mean.' }
  odir = file.path( fishery_model_results, year.assessment, "logistic_discrete_historical" )
  fn1 = file.path( odir, "plot_fishing_mortality_cfanorth.pdf" ) 
  fn2 = file.path( odir, "plot_fishing_mortality_cfasouth.pdf" ) 
  fn3 = file.path( odir, "plot_fishing_mortality_cfa4x.pdf" ) 
include_graphics(c(fn1, fn2, fn3) )
# \@ref(fig:logisticFishingMortality)
```



## Stock status: Reference Points ... {.c}

```{r logistic-hcr, out.width='32%', echo=FALSE, fig.show='hold', fig.align='center', fig.cap = 'Reference Points (fishing mortality and modelled biomass) from Model 1, for N-ENS (left), S-ENS (middle), and 4X (right). The large yellow dot indicates most recent year and the 95\\% CI.' }
  odir = file.path( fishery_model_results, year.assessment, "logistic_discrete_historical" )
  fn1 = file.path( odir, 'plot_hcr_cfanorth.pdf' ) 
  fn2 = file.path( odir, 'plot_hcr_cfasouth.pdf' ) 
  fn3 = file.path( odir, 'plot_hcr_cfa4x.pdf' ) 
  include_graphics(c(fn1, fn2, fn3) )
#  \@ref(fig:logistic-hcr)
```