Generic Species.Rmd

---
  title: "Simple MixFishSim Example"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Simple MixFishSim Example}
%\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
---
  
  This is a simple example of how to use \textbf{'MixFishSim'} to generate
simulations of the dynamics in a mixed fishery. We describe how to calibrate
the habitat fields, the population models, the fishery model and implement a
simple fixed spatial closure. \\

First, load the packages and set a seed for reproducibility.

# Load MixFishSim

```{r packages}
#install.packages('devtools')
#install.packages('githubinstall')

# library(devtools)
# library(githubinstall)
# install_github("pdolder/MixFishSim")
# 
# 
# githubinstall("MixFishSim")


#install.packages('Rcpp')
#install.packages('Rtools')
#install.packages('MixFishSim')
library(MixFishSim)
#library(knitr)
#opts_chunk$set(tidy = TRUE)

#library(reshape2)

#set.seed(123)

```

# Initialise the simulation

This vignette is a paired down example of how to construct a simulation using
MixFishSim. We include only a basic example and encourage users to explore the
other features of the package. \\


## Base parameters

First we specify the basic parameters of the simulation. This includes the
dimensions of the spatial domain, the number of years to simulate, the number
of fleets and vessels per fleet and the number of species and how often (in
                                                                         weeks) the fish move.

The object returned is used internally by MixFishSim a list with two levels: 
  
  * sim$idx : The different units of different processes 
* sim$brk.idx: breaks for each of the key processes in units of a timestep

```{r basic}

#NEW VERSION that has week breaks for entire simulation
# source("R/init_sim_Bens.R")
# sim <- init_sim_Bens(nrows = 100, ncols = 100, n_years = 20, n_tows_day = 1,
#                 n_days_wk_fished = 1, n_fleets = 1, n_vessels = 0, n_species = 2,
#                 move_freq = 1)


#NEW VERSION that has week breaks for entire simulation and allows fishing on just 1 day per week
source("R/init_sim_Bens_nofish.R")
sim <- init_sim_Bens_nofish(nrows = 100, ncols = 100, n_years = 1, n_tows_day = 1,
                            n_days_wk_fished = 1, n_fleets = 1, n_vessels = 1, n_species = 2,
                            move_freq = 1)

class(sim)
sim$idx
names(sim$brk.idx)

```

## Habitat setup

This function creates the spatial fields which support the fish populations and
determine their spatial distributions. You define the parameters for the
matern covariance function for each population and optionally the location of
any spawning closure areas. 

It returns a list of suitable habitat for each species (hab), the habitat as
adjusted during the spawning period (spwn_hab) and the binary location of
spawning areas (spwn_loc). It also returns the locations as x1,x2,y1,y2 and the
multiplier of attractiveness to the spawning area during spawning periods
(spwn_mult).

If plot.dist = TRUE, it returns the plots to a file.

```{r habitat}


source("R/BENS_plot_habitat.R")

#values settled on from anisotropy and habtest scripts

spp.ctrl = list(
  "spp.1" = list('nu' = 1/0.05, 
                 'var' = 1,
                 'scale' = 5, 
                 'Aniso' = matrix(nc = 2,  c(1.5, -3, 3, 4) )),
  "spp.2" = list('nu' = 1/0.015, 
                 'var'  = 1,
                 'scale' = 25, 
                 'Aniso' = matrix(nc = 2,c(1, -2, 1, 2))),
  plot.dist = TRUE, 
  plot.file = "testfolder"
)


#DEFINING UNIFORM NXN GRID
#defining strata coordinates
#number of rows and columns to define (assume want n equal sized strata)
nrows <- 4
ncols <- 4

#of the form  c(x1, x2, y1, y2) THESE ARE BOUNDARIES OF STRATA VALUES
# xs are rows ys columns
strata_num <- 1
stratas <- list()

 for(r in seq(nrows)){
   
   for(c in seq(ncols)){
  
   stratas[[paste("strata",strata_num,sep="")]] = c( (r-1)*(sim$idx[["nrows"]]/nrows)+1, (r)*(sim$idx[["nrows"]]/nrows) , (c-1)*(sim$idx[["ncols"]]/ncols)+1  , (c)*(sim$idx[["ncols"]]/ncols) ) 
  
    strata_num <- strata_num +1  
    
  }
  
}


#DEFINING RANDOM NXN GRID
#size of our domain
totalrows <- 100
totalcols <- 100

#desired dimentions of strata
nrows <- 5
ncols <- 6

#of the form  c(x1, x2, y1, y2) THESE ARE BOUNDARIES OF STRATA VALUES
# xs are rows ys columns
strata_num <- 1
stratas <- list()

#generate randow sequence of rows
rowind <- sort(sample(seq(totalrows),nrows-1,replace=FALSE))
rowind <- c(1,rowind,totalrows)

row_idx <- 0
 for(r in seq(length(rowind)-1)){
row_idx <- row_idx+1

   #generate random sequence of columns
    colind <- sort(sample(seq(totalcols),ncols-1,replace=FALSE))
    colind <- c(1,colind,totalcols)
    
    col_idx <- 0
    
   for(c in seq((length(colind)-1))){
     
     col_idx <- col_idx+1
  
   stratas[[paste("strata",strata_num,sep="")]] = c( rowind[row_idx] , rowind[row_idx+1] , colind[col_idx] , colind[col_idx+1]) 
  
    strata_num <- strata_num +1  
    
  }
  
}


source("R/BENS_create_hab.R")

hab <- BENS_create_hab(sim_init = sim, 
                       spp.ctrl = spp.ctrl,
                       
                       spawn_areas = list(
                         "spp1" = list(
                           'area1' = c(30,45,55,65),   #of the form  c(x1, x2, y1, y2) THESE ARE BOUNDARIES OF MATRIX VALUES
                           'area2' = c(70,90,50,60)   #need to revisit closure areas
                         ),
                         "spp2" = list(
                           'area1' = c(30,45,55,65),
                           'area2' = c(70,90,50,60)
                         )),
                       spwn_mult = 10,#THIS DOES NOT CHANGE ANYTHING. MUST CHANGE VALUE AT TOP OF BENS_CREATE_HAB 
                       plot.dist = F, 
                       plot.file = "testfolder",
                       
                       
                       #created this new part defining strata
                      strata = stratas

                       )

#plot strata
par(mar=c(1,1,1,1))
fields::image.plot(hab$stratas)


#print(hab)
source("R/BENS_plot_habitat.R")
## Plot the unadjusted habitat fields
BENS_plot_habitat(hab = hab$hab, 	spp.ctrl = spp.ctrl)



#old version
plot_habitat(hab$hab)

## Plot the adjusted habitat fields
plot_habitat(hab$spwn_hab)

```

## Population models

Now we need to set up the population models for the simulations. We do this
with the init_pop function. We set the initial population biomasses, movement
rates, recruitment parameter and growth and natural mortality rates.

The object created stores all the starting conditions and containers for
recording the changes in the populations during the simulations.

We can plot the starting distributions for each population as a check.

```{r pop_init}

#load rcpp exports
Rcpp::sourceCpp(file= "src/Movement.cpp")
Rcpp::sourceCpp(file= "src/RcppExports.cpp")
Rcpp::compileAttributes() #this updates RcppExports.R file, which contains function definitions



#CALULATE INDICES OF NONZERO VALUES IN HAB TO PASS TO MOVE_POPULAITON DURING MOVEMENT
nonzero_idx <- lapply(paste0("spp", seq_len(sim$idx[["n.spp"]])), function(s) {
  
  which(hab[["hab"]][[s]] >0 , arr.ind=T)
  
})

names(nonzero_idx) <- paste("spp",seq_len(sim$idx[["n.spp"]]), sep ="")

#Week 13 is first week of april
#half way through april to half way through may would be week 15-18

source("R/init_pop_Bens.R")

#original settings
# Pop <- init_pop_Bens(sim_init = sim, Bio = c("spp1" = 2e5, "spp2" = 4e5), #these values from paper : 1e5 and 2e5
# 		hab = hab[["hab"]], start_cell = c(50,50), 
# 		lambda = c("spp1" = 0.1, "spp2" = 0.1), #same lambda for all?
# 		init_move_steps = 20,
# 		rec_params = list("spp1" = c("model" = "BH", "a" = 6, "b" = 4, "cv" = 0.7),
# 				  "spp2" = c("model" = "BH", "a" = 27, "b" = 4,"cv" = 0.6)), #these values from paper 
# 		rec_wk = list("spp1" = 12:15, "spp2" = 12:15),
# 		spwn_wk = list("spp1" = 15:18, "spp2" = 15:18),
# 		M = c("spp1" = 0.2, "spp2" = 0.1), #these values from paper: c("spp1" = 0.2, "spp2" = 0.1)
# 		K = c("spp1" = 0.3, "spp2" = 0.3) #all the same for now
# 		)

#decreasing population settings
Pop <- init_pop_Bens(sim_init = sim, Bio = c("spp1" = 4e5, "spp2" = 10e5), #these values from paper : 1e5 and 2e5
                     hab = hab[["hab"]], start_cell = c(50,50),
                     lambda = c("spp1" = 0.1, "spp2" = 0.1), #same lambda for all?
                     init_move_steps = 20,
                     rec_params = list("spp1" = c("model" = "BH", "a" = 2, "b" = 4, "cv" = 0),
                                       "spp2" = c("model" = "BH", "a" = 7, "b" = 4,"cv" = 0)), #these values from paper
                     rec_wk = list("spp1" = 15:18, "spp2" = 15:18),
                     spwn_wk = list("spp1" = 15:18, "spp2" = 15:18),
                     M = c("spp1" = 0.275, "spp2" = 0.225), #these values from paper: c("spp1" = 0.2, "spp2" = 0.1)
                     K = c("spp1" = 0.3, "spp2" = 0.3), #all the same for now
                     nz = nonzero_idx
                     )



#I PUT THIS BACK INSIDE RUN_SIM
# #Calculate movement probabilities (used to be in run_sim)
# Move_Prob  <- lapply(paste0("spp", seq_len(sim[["idx"]][["n.spp"]])), function(s) { move_prob_Lst(lambda = Pop[["dem_params"]][[s]][["lambda"]], hab = hab[["hab"]][[s]])})
# 
# 
# Move_Prob_spwn <- lapply(paste0("spp", seq_len(sim[["idx"]][["n.spp"]])), function(s) { move_prob_Lst(lambda = Pop[["dem_params"]][[s]][["lambda"]], hab = hab[["spwn_hab"]][[s]])})
# 
# 
# names(Move_Prob)      <- paste0("spp", seq_len(sim[["idx"]][["n.spp"]]))
# names(Move_Prob_spwn) <- paste0("spp", seq_len(sim[["idx"]][["n.spp"]]))





names(Pop)

Pop$dem_params

par(mfrow = c(2,1))

image(Pop$Start_pop[[1]], main = "spp1 starting biomass")
image(Pop$Start_pop[[2]], main = "spp2 starting biomass")




```

## Population movement

Now we set up the population tolerance to different temperatures which
determines how the populations move during the course of a year. We can then
plot the combined spatiotemporal suitable habitat to examine how these
interact.

```{r temp}

#set temperature preferences manually. 
#The following assumes moveCov has been created and already has an empty spp_tol sublist
moveCov[["spp_tol"]] <- list() #just in case
moveCov[["spp_tol"]] <- list("spp1" = list("mu" = 7.98, "va" = 3),  #8.13 IF TEMP INCREASES   7.98 if temp constant
                             "spp2" = list("mu" = 7.98, "va" = 3) )



plot(norm_fun(x = 0:25, mu = 8.3, va = 3)/max(norm_fun(0:25, 8.3, 3)),
     type = "l", xlab = "Temperature", ylab = "Tolerance", lwd = 2)
lines(norm_fun(x = 0:25, mu = 15, va = 9)/ max(norm_fun(0:25, 15, 9)),
      type = "l", col = "blue", lwd = 2)
lines(norm_fun(x = 0:25, mu = 17, va = 7)/ max(norm_fun(0:25, 17, 7)),
      type = "l", col = "green", lwd = 2)

legend(x = 2, y = 0.9, legend = c("spp1", "spp2"), lwd = 2, col = c("black", "blue"))


# plot_spatiotemp_hab(hab = hab, moveCov = moveCov, spwn_wk = list("spp1" = 15:18, "spp2" = 15:18), plot.file =  "testfolder")
# 
# 
# 
# #to plot just the temp preferences over time
# source("R/BENS_plot_spatiotemp_hab_justtemp.R")
# 
# BENS_plot_spatiotemp_hab_justtemp(hab = hab, moveCov = moveCov, spwn_wk = list("spp1" = 16:18, "spp2" = 16:19,"spp3" = 16:18, "spp4" = 18:20), plot.file =  "testfolder")
# 

```

## Fleet models

Here we initialise the fleet with fish landings price per tonne, catchability
coefficients per population, fuel cost, the coefficients for the step function
and fleet behaviour.

We can plot the behaviour of the step function to check its suitable for our
simulations. This determines the relationship between the monetary value gained
from a fishing tow and the next move by the vessel when using the correlated
random walk function.

```{r fleets}
#initial settings
# fleets <- init_fleet(sim_init = sim, VPT = list("spp1" = 4, "spp2" = 3),
# 		     Qs = list("fleet 1" = c("spp1" = 1e-5, "spp2" = 3e-5),
# 			       "fleet 2" = c("spp1" = 5e-5, "spp2" = 1e-5)
# 			       ),
# 		     fuelC = list("fleet1" = 3, "fleet 2" = 8),
# 		     step_params = list("fleet 1" = c("rate" = 3, "B1" = 1, "B2" = 2, "B3" = 3),
# 					"fleet 2" = c("rate" = 3, "B1" = 2, "B2" = 4, "B3" = 4)
# 					),				
# 		     past_knowledge = TRUE,
# 		     past_year_month = TRUE,
# 		     past_trip = TRUE,
# 		     threshold = 0.7
# 		     )


#no fishing
fleets <- init_fleet(sim_init = sim, VPT = list("spp1" = 0, "spp2" = 0), #VPT = value per ton
                     Qs = list("fleet 1" = c("spp1" = 0, "spp2" = 0)   #Q = catchability
                     ),
                     fuelC = list("fleet1" = 3),
                     step_params = list("fleet 1" = c("rate" = 3, "B1" = 1, "B2" = 2, "B3" = 3)
                     ),				
                     past_knowledge = FALSE,  #dont use past knowledge
                     past_year_month = TRUE,
                     past_trip = TRUE,
                     threshold = 0.7
)



test_step(step_params = fleets$fleet_params[[1]]$step_params, rev.max = 1e2)
test_step(step_params = fleets$fleet_params[[2]]$step_params, rev.max = 1e2)

```

## Spatial closure

We set up a spatial closure. There are multiple options in defining
this, but we simply define a static fixed site closure for demonstration
purposes. 

```{r close}
# #practice creating a larger data.frame than just single points
# library(tidyr)
# 
# #set x and y min/max which are coordinates on the grid
# xmin <- 6
# xmax <- 10
# ymin <- 26
# ymax <- 40
# 
# 
# x <- xmin:xmax
# y<- ymin:ymax
# 
# #View(crossing(x,y))
# 
# closure <- init_closure(input_coords = data.frame(x = x, y = y),
# 			spp1 = "spp1", year_start = 2) 

```

## Survey

Its also possible to define a survey design using the init_survey function, but
we do not do so for this demonstration. Please refer to the function help file
if this is required.

```{r survey}

source("R/BENS_init_survey.R")


#CURRENTLY NEED TO MAKE SURE THAT N_STATIONS*#YEARS / #STRATA IS A WHOLE NUMBER OTHERWISE DAY, TOW, YEAR WONT LINEUP WITH NUMBER OF STATIONS
#ALSO NEED N_STATION TO BE DIVISIBLE BY STATIONS_PER_DAY
#ALSO NEED N_STATIONS / STATIONS_PER_DAY <= 52 otherwise wont get to all of them in a year results in NA in the matrix

#setup catch log
surv_random <- BENS_init_survey(sim_init = sim,design = 'random_station', n_stations = 80, #this is total per year (20 in each of 4 strata)
                                start_day = 1, stations_per_day = 1, Qs = c("spp1" = 0.1, "spp2"= 0.2),
                                strata_coords = hab$strata, strata_num = hab$stratas )

```

## Run simulation

Finally we run the simulation. The output is a list of objects containing all
the information on fisheries catches, the population dynamics and population
distributions. These can be examined with some inbuilt plotting functions. 

```{r sim}
run_simulation <- function(x){
  
  source("R/run_sim.R")
  
  
  #to source a new go_fish where I edited to skips most things:
  #1: load file
  source("R/go_fish_Bens.R") #my edited version that skips most things
  #2: allow the function to call other hidden functions from mixfishsim 
  environment(go_fish_Bens) <- asNamespace('MixFishSim')
  #3: replace go_fish with go_fish_Bens in the MixFishSim package
  assignInNamespace("go_fish", go_fish_Bens, ns = "MixFishSim")
  
  
    source("R/RcppExports_Bens.R")
  
    #to source a new move_population where I edited to skips most things:
  #1: load file
  Rcpp::sourceCpp("src/Movement_Bens.cpp") #my edited version that skips most things
  #2: allow the function to call other hidden functions from mixfishsim 
  environment(move_population_Bens) <- asNamespace('MixFishSim')
  #3: replace move_population with go_fish_Bens in the MixFishSim package
  assignInNamespace("move_population", move_population_Bens, ns = "MixFishSim")
  
  
  
  library(MixFishSim) #each core needs to load library
  
  #load data from CPU
  #  load("C:/Users/benjamin.levy/Desktop/Github/READ-PDB-blevy2-toy/20 year moveCov matrices/Final/Final_moveCov_12_9_Bensmethod.RData")
  
  #load constant temp data from Mars
  #load("/net/home5/blevy/Bens_R_Projects/READ-PDB-blevy2-toy/20 year moveCov matrices/Final/Final_constanttemp_20yr.RData")
  
  #load increase temp data from Mars
  load("/net/home5/blevy/Bens_R_Projects/READ-PDB-blevy2-toy/20 year moveCov matrices/Final/Increase_temp_20yr_LargerTempVA7.RData")
  
  
  res<- run_sim(sim_init = sim,
                pop_init = Pop,
                move_cov = moveCov,
                fleets_init = fleets,
                hab_init = hab,
                save_pop_bio = TRUE,
                survey = NULL, #surv_random, will try to survey after simulation
                closure = NULL,
                InParallel = TRUE,
                nz = nonzero_idx)  #does it runin parallel? Doesnt seem like it
  
  
}



#run in parallel

library(parallel)

nCoresToUse <- detectCores() - 1   #this show 16 cores but I think I have 6??

nCoresToUse <- 6


cl <- parallel::makeCluster(nCoresToUse,revtunnel = TRUE, outfile = "", verbose = TRUE, master=nsl(Sys.info()['nodename'])) #options from https://stackoverflow.com/questions/41639929/r-cannot-makecluster-multinode-due-to-cannot-open-the-connection-error


result <- list()


result <- parallel::parLapply(cl,1:6,run_simulation)
parallel::stopCluster(cl)


# 
# #below throws an error
# #result <- foreach(i=1:3) %dopar% run_simulation(i)
# 
# #below throws a similar error: 3 nodes produced errors; first error: missing value where TRUE/FALSE needed
#  parLapply(cl,1:3,run_simulation)
# stopCluster(cl)


```

## Summary plots



There are a series of input plotting functions to visualise the results of the
simulation. For example, we can explore:
  
  * the population dynamics for each species
* Seasonal patterns in exploitation
* the location choice of a vessel
* the realised step function for a vessel

Users will wish to define their own plots, depending on the issues of interest
and all the results are saved in the output from the run_sim function.

```{r plots}

## Biological
source("R/plot_pop_summary.R")
p1 <- plot_pop_summary(results = res,
                       timestep = "annual",
                       save = FALSE, 
                       save.location = NULL
)

plot_pop_summary(results = res, timestep = "daily", save = FALSE, save.location = "C:/Users/benjamin.levy/Desktop/Github/READ-PDB-blevy2-toy/testfolder")
p1

p2 <- plot_daily_fdyn(res)
p2

## Fishery

logs <- combine_logs(res[["fleets_catches"]])

p3 <- plot_vessel_move(sim_init = sim, logs = logs, fleet_no = 1, vessel_no = 5,
                       year_trip = 5, trip_no = 10)
p3

p4 <- plot_realised_stepF(logs = logs, fleet_no = 1, vessel_no = 1)
p4      



#try some new ones

#plot survey 2 ways
p5 <- plot_survey(survey = res$survey, type = "spatial")
p5

p6 <- plot_survey(survey = res$survey, type = "index")
p6


#spatio temporal population plots
source("R/Bens_plot_pop_spatiotemp.R")
p7 <- Bens_plot_pop_spatiotemp(results = res,
                               save = TRUE, 
                               save.location = "testfolder/spatialplots"
)
p7


```

Note in our example how the fishing mortality rate for species 2 changes
following the spatial closure, which was set to cover some of the core
distribution of the population.