inst/markdown/04.snowcrab_fishery_model_turing.md

In jae0/snowcrab: Snow crab assessment suite for aegis* and associated projects.

title: "Dynamical Fishery Model" author: "Jae S. Choi" toc: true number-sections: true highlight-style: pygments editor: render-on-save: false format: html: code-fold: true html-math-method: katex self-contained: true pdf: pdf-engine: lualatex docx: default

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 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.

#| 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( "1999_present", snowcrab_filter_class, sep="_" )
  fishery_model_label = "turing1"
  carstm_results_directory = file.path( modeldir, carstm_model_label, "fishery_model_results", fishery_model_label ) # default

  fishery_model_data_inputs( year.assessment=year.assessment,  type="biomass_dynamics", for_julia=TRUE, 
    save_location=carstm_results_directory, fishery_model_label=fishery_model_label  ) 

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

  # for development, save at alternate locations
  # carstm_results_directory = file.path( homedir, "projects", "dynamical_model", "snowcrab", "data" )
  # fishery_model_data_inputs( year.assessment=year.assessment,  type="biomass_dynamics", for_julia=TRUE, save_location=carstm_results_directory )

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, where more options are also shown. But it will require learning the language Julia and Turing library, but is very simple.

  year_assessment =2023
  yrs = 1999:year_assessment

  rootdir = homedir()
  project_directory = joinpath( rootdir, "bio", "bio.snowcrab", "inst", "julia" ) 
  bio_data_directory = joinpath( rootdir, "bio.data", "bio.snowcrab", "modelled", "1999_present_fb", "fishery_model_results", "turing1"  )
  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 )

  # choose:
  aulab ="cfanorth"     # 41 hrs (size struct), 20 min (logistic_discrete_hist)
  aulab ="cfasouth"     # 65 hrs (size struct), min (logistic_discrete_hist)  
  aulab ="cfa4x"        # 24 hrs (size struct), min (logistic_discrete_hist)

  include( joinpath(project_directory, "snowcrab_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
  include( joinpath(project_directory, "logistic_discrete.jl" ) )   

OR,

1. 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.

#| 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 with report generation.

For example, using Rmarkdown reports.

Here some of the outputs can be seen below:

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', '1999_present_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)

jae0/snowcrab documentation built on Feb. 27, 2024, 2:42 p.m.