Simple_MixFishSim_Example_RecreatePublication.Rmd

---
title: "Simple MixFishSim Example"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Simple MixFishSim Example}
    %\VignetteEngine{knitr::rmarkdown}
      \usepackage[utf8]{inputenc}
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}

sim <- init_sim(nrows = 100, ncols = 100, n_years = 20, n_tows_day = 4,
		n_days_wk_fished = 5, n_fleets = 2, n_vessels = 20, n_species = 4,
	  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_create_hab.R")

source("R/BENS_plot_habitat.R")

#setup controls outside to use in muiltiple places


# #original from example script
# 		  spp.ctrl = list(
# 				  "spp.1" = list('nu' = 1/0.015,
# 						'var' = 1,
# 						'scale' = 1,
# 						'Aniso' = matrix(nc = 2, c(1.5, 3, -3, 4))),
# 				  "spp.2" = list('nu' = 1/0.05,
# 						'var'  = 2,
# 						'scale' = 12,
# 						'Aniso' = matrix(nc = 2, c(1, 2, -1, 2))),
# 				 plot.dist = TRUE,
# 				 plot.file = "C:/Users/benjamin.levy/Desktop/Github/READ-PDB-blevy2-toy/testfolder"
# 				 )


#values from the MixFishSim publication (two are those given above in the example)
		  spp.ctrl = list(
				  "spp.1" = list('nu' = 1/0.015, #inconsistent btwn paper and moveCovTest
						'var' = 1,
						'scale' = 10, #given as 10 in FullTest.R
						'Aniso' = matrix(nc = 2, c(1.5, 3, -3, 4))),
				  "spp.2" = list('nu' = 1/0.05, #inconsistent btwn paper and moveCovTest
						'var'  = 2,
						'scale' = 20, #given as 20 in FullTest.R
						'Aniso' = matrix(nc = 2, c(1, 2, -1, 2))),
				 "spp.3" = list('nu' = 1/0.01, #inconsistent btwn paper and moveCovTest
						'var' = 1,
						'scale' = 5, #given as 5 in moveCovTest.R
						'Aniso' = matrix(nc = 2, c(2.5, 1, -1, 2))),
				 "spp.4" = list('nu' = 1/.005, #inconsistent btwn paper and moveCovTest
						'var' = 1,
						'scale' = 30, #given as 30 in moveCovTest.R
						'Aniso' = matrix(nc = 2, c(0.1, 2, -1, 0.2))),
				
				 		 plot.dist = TRUE, 
				 plot.file = "testfolder"
				 )

		  
hab <- BENS_create_hab(sim_init = sim, 
		  spp.ctrl = spp.ctrl,
		  
		  spawn_areas = list(
				     "spp1" = list(
						   'area1' = c(40,50,40,50),   #of the form  c(x1, x2, y1, y2) THESE ARE BOUNDARIES OF MATRIX VALUES
						   'area2' = c(80,90,60,70)    #of the form  c(x1, x2, y1, y2) IE, X FROM 2 TO 5 AND Y FROM 2 TO 5
						   ),
				     "spp2" = list(
						   'area1' = c(50,60,30,40),
						   'area2' = c(80,90,90,90)
						   ),
				     	 "spp3" = list(
						   'area1' = c(30,34,10,20),   #of the form  c(x1, x2, y1, y2) THESE ARE BOUNDARIES OF MATRIX VALUES
						   'area2' = c(60,70,20,30)    #of the form  c(x1, x2, y1, y2) IE, X FROM 2 TO 5 AND Y FROM 2 TO 5
						   ),
				     "spp4" = list(
						   'area1' = c(50,55,80,85),
						   'area2' = c(30,40,30,40)
						   ),
				     spwn_mult = 10, 
				    plot.dist = TRUE, 
				    plot.file = "testfolder" ),
		  
		  
		  #created this new part defining strata
		  
		  strata = list(                              #Form: upper left, upper right, lower left, lower right
				     "strata1" = c(1,sim$idx[["nrows"]]/2,1,sim$idx[["ncols"]]/2),   #of the form  c(x1, x2, y1, y2) THESE ARE BOUNDARIES OF STRATA VALUES
				      
				      "strata2" = c(1,sim$idx[["nrows"]]/2,sim$idx[["ncols"]]/2+1,sim$idx[["ncols"]]),
				     
				     	"strata3" = c(sim$idx[["nrows"]]/2 + 1,sim$idx[["nrows"]], 1 ,sim$idx[["ncols"]]/2),
				     
              "strata4" = c(sim$idx[["nrows"]]/2 + 1,sim$idx[["nrows"]],sim$idx[["ncols"]]/2 +1 , sim$idx[["ncols"]])
				       )
		  )
		  

#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}

Pop <- init_pop(sim_init = sim, Bio = c("spp1" = 1e5, "spp2" = 2e5,"spp3" = 1e5, "spp4" = 1e4),
		hab = hab[["hab"]], start_cell = c(5,5), #single start cell for all species?
		lambda = c("spp1" = 0.1, "spp2" = 0.1,"spp3" = 0.1, "spp4" = 0.1), #same lambda for all?
		init_move_steps = 40,
		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),
				  "spp3" = c("model" = "BH", "a" = 18, "b" = 11,"cv" = 0.7),
				  "spp4" = c("model" = "BH", "a" = 0.3, "b" = 0.5,"cv" = 0.6)),
		rec_wk = list("spp1" = 13:16, "spp2" = 12:16,"spp3" = 14:16, "spp4" = 16:20),
		spwn_wk = list("spp1" = 16:18, "spp2" = 16:19,"spp3" = 16:18, "spp4" = 18:20),
		M = c("spp1" = 0.2, "spp2" = 0.1,"spp3" = 0.2, "spp4" = 0.1),
		K = c("spp1" = 0.3, "spp2" = 0.3,"spp3" = 0.3, "spp4" = 0.3) #all the same for now
		)

names(Pop)

Pop$dem_params

par(mfrow = c(2,2))

image(Pop$Start_pop[[1]], main = "spp1 starting biomass")
image(Pop$Start_pop[[2]], main = "spp2 starting biomass")
image(Pop$Start_pop[[3]], main = "spp3 starting biomass")
image(Pop$Start_pop[[4]], main = "spp4 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}




moveCov <- init_moveCov(sim_init = sim, steps = 52,
			spp_tol = list("spp1" = list("mu" = 12, "va" = 8),
				       "spp2" = list("mu" = 15, "va" = 9), #given as 7 in code but 9 in paper
				       "spp3" = list("mu" = 17, "va" = 7),
				       "spp4" = list("mu" = 14, "va" = 10)
				       )
			)

plot(norm_fun(x = 0:25, mu = 12, va = 8)/max(norm_fun(0:25, 12, 8)),
     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)
lines(norm_fun(x = 0:25, mu = 14, va = 10)/ max(norm_fun(0:25, 14, 10)),
      type = "l", col = "red", lwd = 2)
legend(x = 2, y = 0.9, legend = c("spp1", "spp2"), lwd = 2, col = c("black", "blue","green","red"))


plot_spatiotemp_hab(hab = hab, moveCov = moveCov, spwn_wk = list("spp1" = 16:18, "spp2" = 16:19,"spp3" = 16:18, "spp4" = 18:20), 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}

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
		     )

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")


surv_fixed <- BENS_init_survey(sim_init = sim,design = 'fixed_station', n_stations = 48, start_day = 1, stations_per_day = 7, Qs = c("spp1" = 0.1, "spp2"= 0.2))

#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


surv_random <- BENS_init_survey(sim_init = sim,design = 'random_station', n_stations = 48, start_day = 1, stations_per_day = 7, 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}

source("R/run_sim.R")

res <- run_sim(sim_init = sim,
	       pop_init = Pop,
	       move_cov = moveCov,
	       fleets_init = fleets,
	       hab_init = hab,
	       save_pop_bio = TRUE,
	       survey = surv_fixed,
	       closure = closure,
	       InParallel = TRUE)

```

## 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
p1 <- plot_pop_summary(results = res,
		 timestep = "annual",
		 save = TRUE, 
		 save.location = "C:/Users/benjamin.levy/Desktop/Github/READ-PDB-blevy2-toy/testfolder"
		 )

plot_pop_summary(results = res, timestep = "daily", save = TRUE, 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.