TasmModel_RunningCode.Rmd

---
title: "Tasmanian model run"
author: "Asta"
date: "4 Feb 2019"
output: 
  html_document:
    code_folding: show
    toc: FALSE
    toc_float: FALSE
  pdf_document:
    toc:no
---

### Thoughts and tasks 

Mar 5: Chat with Scott Ling about urchin parameters 

Spines are ca 30% of body weight, so discount that (and water inside the urchin) from body mass, as it is really not available for consumption and does not really determine their feeding preferences. 
Centrstephanus - Wmax = 700g, Wmat = 65g with spines. For the model purpose therefore 
Wmax = 500g, Wmat = 45g
Heliocedaris in tasmania, Wmax = 300g, Wmat = 45g, average weight of an individual = 150g. Without spines Wmax = 200g, Wmat = 30g, average weight ca 100g (so an average individual is quite big!)
Heliocedaris in Port Phillip Bay, Wmax (with spines) = 210g, Wmat = 25, average weight = 85 

Since we pool the urchins and have an average set of parameters we use 
Wmax = 350, Wmat = 40

Maximum intake:

Maximum intake rates have been calcualted based on gut weight and evacuation rate of 5 days (determined in experiments). For Centrostephanus daily intake then is = 0.01-0.03 (mean of 0.02) g/g/day. For Heliocedaris = 0.026g/g/day (0.006-0.04 for min-max gut weight versus body weight) in Tasmania and 0.032g/g/day in Port Phillip bay. 

In mizer parameters are given per year, so average daily intake of 0.02g/g/day would become 7.3 g/g/year. This seems way too small. In my current setup h for urchins is 145, so a 50g urchin would eat (145*50^0.8)/365= 9g of algae per day. Seems too much. If h is reduced in half 
(75*50^0.8)/365  = 4.7 g per day or 8-9% of body weight. 

Olistophs wmat of 440g (say 500g) and h of 180, means their maximum possible intake per year would be 
180*500^0.8 = 26000g or 26 kg. Per day it is 70g per 500 g fish or 14% of body mass. If h is 90, then maximum intake is 
(90*500^0.8)/365 = 35g or 7% of body weight. This seems more reasonable. 

For D_lewini of 250g and h = 50, the max intake per day is (50*250^0.8)/365 = 11g of fish per individual. This seems totally reasonable, in fact possibly too small... 

**

Feb 25: The difference in long term projections is big if we use 200 or 1000 size groups. So the number of size groups is a very important parameter and should be explored as such. Same applies to the time step or dt

Feb 15: After talking to Scott I realised that urchins are much too small! They can grow to 900g in size!

Issues as of Feb 11

Added urchins and lobsters. They all seem to coexist ok, but the problem is with diets and the fact that herbivores don't eat macroalgae and exclusive benthivores don't eat benthos anymore. Lobsters eat all the benthos, it seems. 

Meeting on Feb 6, 2019

Lobsters eat small fish and urchins. They used to be the dominant predator of small fish (hypothesis)
Southern and eastern rock lobster 
Southern calamari and snapper thermal preference 
Use relative abundance and use species specific size spectra from the first 5 years 

1. Fishing mortality using data inside and outside MPA
2. Benthic resource slope and regeneration rate and how they change with temperature and what does it mean to the community 
3. How does temperature alter the community structure inside MPA using long term data

how are these systems functioning?
how does fishing and climate affect that?

Adding urchins and lobsters
For urchins h based on standard equation is 26, but given that they are herbivores it should be at least 3 time higher. Still it only is 75, while h for herbivorous fish is closer to 180. So if urchin h = 200, it means they can eat 200g per 1 g per year, or 2kg per 10g urchin per year, or 5.5g per day. Need references from herbivory studies 

Parameterisation thoughts

What will control species calibrated through erepro? If their recruitment is directly proportional to spawn, then with fishing mortality added the only way it would not lead to crashing biomass, is through density dependence processes. This means that when biomass is removed through fishing it has to be compensated by emergent increased recruitment. This will only work through competition processesses and starvation mortality. If under fishing smaller number of larvae compete less for food, they will grow faster and will have less background mortality before reproduction. To get this happening we want to make sure that background food can be limited and lead to starvation mortality. 

Generally we find that erepro for most species is around 0.0001 but is much higher (0.04) in Trachinus. I suspect this is because they have such small body size and not enough size classes between maturation and w_inf to produce lots of eggs. This is also partly because psi function is identical among species when I imagine such small species shift allocation much faster to reproduction. At the momnet I use exponent 5 in the psi function (mizer default is 10). Chaning the exponent for all species (in MizerParams-class.R) does not make much effect on Trachinus because it just increases allocation speed in all species. So it makes sense that their erepro is much higher, as it would encompass higher allocation to reproduction (sort off...)

Pooling of species 

The two Notolabrus (wrasse) species are now pool into one group, as they are too similar for the model to coexist 
Nev's thesis say that M. australis and A. vittiger both grow to about 320mm and both feed on invertebrates, such as molluscs, echinoids and poriferans. M. australis has some algae in its stomachs but Nev's view is that algae are not digested but are accidentally swallowed. For M. australis the MP(5-95) is 14.8 for A. vittiger it is 17.3. Growth trajectories for the two species in Nev's thesis also look similar, they reach about 15-18cm by the age of 2, and 30 cm by the age of 5-6 and older ages are not seen. I pool them into leather_cool category


### Install mizer/rewiring  

```{r warning=F, message=FALSE, echo=T}

rm(list=ls())

#install.packages("mizer")
#library(devtools)
#install_github("astaaudzi/mizer", ref = "rewire-temp")
library(mizer)

```

### Data files

```{r warning=T, message=FALSE, echo=T}
mariaParams <- read.csv(file = "Tasm_ParametersMar6.csv") #16 species
#mariaParams <- read.csv(file = "mariaParamsMar6.csv") #16 species
inter <- read.csv(file = "Tasm_Interaction.csv") #16 species
inter <- as.matrix(inter)

#number of size groups 
no_size_groups = 200
#timestep used in the integration 
dt = 0.2
## scalar to get initial abundance close to what is observed per survey. Initial abundances are calculated based on n, q, w_inf, and kappa, so they should not change during the optimisation. Basically the default initial abundances for small species are way too high and for large species too small. 
scalarForInitial <-c(38.5 ,  0.3, 403.4,  10.6 ,  5.0  , 3.4 ,  0.2 ,  0.01, 114.5 ,  4.9  , 1.5 , 15.9,  10.3 ,  2.5 ,  0.4  , 7.1  , 0.4)

```

### Background spectra

The slope for the plankton spectrum is assumed to be a standard 2. 
For the benthic spectrum I use AbNoUr <- 0.8 - 0.85*log10(InvDataBinned$wgtGroup) equation and calibrate the kappa to ensure initial abundances are within the range of observed abundances per m2

```{r warning=T, message=FALSE, echo=T}

## resource params
kappa = 20 # 20 # intercept assuming g/m2
lambda = 2.1 # 
w_pp_cutoff = 1 #g
r_pp = 2 #2 # rate of regeneration
min_w_pp = 1e-10 #g

kappa_ben = 80 #80 # intercept assuming g/m2  
lambda_ben = 1.85 #this slope does not include urchins and lobsters
w_bb_cutoff = 50 #
r_bb = 1 # 1.5 # something to be calibrated. Default mizer option is 10
min_w_bb = 0.001 # 0.01

kappa_alg = 100 # 100 #intercept assuming g/m2
lambda_alg = 1.6 #slope is much more shallow, to allow for lots of large kelp, but size structure does not really make much sense here
w_aa_cutoff = 100 
r_aa = 2 #1 #something to be calibrated
min_w_aa = 0.001

```

### Fixed params

```{r warning=T, message=FALSE, echo=T, eval = F}

#### egg size, rmax, exponents ####
#mariaParams$w_min <- (0.0015* mariaParams$w_inf^0.14)*10  # Barneche et al. 2018 equation for egg sizes, but let's start with at least 10 times bigger individuals, as they don't start as eggs

mariaParams$w_min <- 0.05 
mariaParams$w_min[which(mariaParams$species == "C_laticeps")] <- 1  #assume that shark eggs are 5 grams?

mariaParams$n <- 0.75 # maximum intake scaling, mizer default is 2/3, but the growth curves don't look so good
mariaParams$p <- 0.8 # metabolism scaling
mariaParams$q <- 0.75 # search rate scaling

#### h & ks  ####

# by default mizer will set  h values, calcualted as h = [3*k_vb/(alfa*f0)]*W_inf^(1/3)
# I use size dependent scaling, derived from earlier explorations 

mariaParams$h = 50 * (mariaParams$w_inf/1000)^0.15
# then for herbivores, who have low assimilation efficiency, we double the maximum intake rate (based on Graham Edgar daily intake calculations, ref in Soler et al. 2018)
mariaParams$h[which(mariaParams$funcgr == "herbi")] <- mariaParams$h[which(mariaParams$funcgr == "herbi")]*2
#See thoughts above on whether these values give reasonable estimates

#set size dependent ks based on my earlier explorations 
mariaParams$ks <- 20*mariaParams$w_inf^(-0.25)
#however this relationship gives Trachinops very high ks, so I decrease it. The relationship seems too steep for smallest species. Even though small species have high per unit body mass ks, this apppears much too high
mariaParams$ks[which(mariaParams$species == "T_caudimaculatus")] <- 4.5 #instead of 8

## Is ks for urchins too high? 
#which(mariaParams$species == "urchins")
#mariaParams$h[16]
# Let's assume they are are feeding at maximum possible rate
#energy <- (mariaParams$h[16]*mariaParams$w_mat[16]^mariaParams$n[16])*mariaParams$alpha[16]
#metab <- mariaParams$ks[16]*mariaParams$w_mat[16]^mariaParams$p[16]
#energy- metab # plenty of net energy left

### Background mortality 

# Version 1
#I am now using mizer default values for background mortality, as z0 = z0pre * w_inf ^ (n-1)
#mariaParams$z0 = mariaParams$z0pre * mariaParams$w_inf ^ (mariaParams$n -1 )
#min(mariaParams$z0)

# Version 2
#Actually I am not using mizer default, because I am not sure why background mortlaity should depend on n. I am just using exponent of -0.3, but scale down Trachinops mortality from 0.2 to 0.15
#mariaParams$z0 = (mariaParams$z0pre) * mariaParams$w_inf ^ (-0.3)
#min(mariaParams$z0)
#max(mariaParams$z0)
#plot(mariaParams$w_inf, mariaParams$z0)
#exponential scaling gives very high background mort for Trachinops, so I set it lower 
#mariaParams$z0[which(mariaParams$species =="T_caudimaculatus")] <- 0.15
#mariaParams$z0[which(mariaParams$species =="P_laticlavius")] <- 0.08

# Version 3 
# Mortality of smallest species is already very high from predation, so I am not sure why background mortality for large species should be lower (as in mizer default), especially when they are in larval stage. In fact there is evidence that fast growth results in higher mortality, and this would apply for larvae of large species. So I am setting background mortality to same value for all species. The initial value of background mortality of 0.1 I used is actually quite high, in Atlantis we tried to use values of 0.0001 or lower.  So I will set it at 0.05 for now
mariaParams$z0 <- 0.05

```

### Dynamic parameters: gamma and erepro

```{r warning=T, message=FALSE, echo=T, eval = F}

## Gamma 

## In the MizerParams we asked mizer to set gamma's based on the defaults: 
## gamma = f0*h*beta^(2-lambda)/((1-f0)*2*3.14*kappa*sigma)

mean_lam = mean(lambda, lambda_ben, lambda_alg)
mean_kappa = mean(kappa, kappa_ben, kappa_alg)

#as an initial setup we use default mizer gamma parameters, but note this is a very crude approximation because we have 3 size spectra and species have different access to these spectra
gammaDef = mariaParams$f0 *mariaParams$h*mariaParams$beta^(2-mean_lam)/((1-mariaParams$f0)*2*3.14*mean_kappa*mariaParams$sigma)

#for herbivores gamma values are doubled
#gammaDef[which(mariaParams$funcgr == "herbi")] <- gammaDef[which(mariaParams$funcgr == "herbi")] *2

# for predators gamma values should be much higher or else they don't get enough food, because they generally search for food actively 
gammaDef[which(mariaParams$funcgr == "predat")] <- gammaDef[which(mariaParams$funcgr == "predat")] *4
mariaParams$gamma <- gammaDef *20

areaSearchedPerDay <- (mariaParams$gamma*mariaParams$w_mat^mariaParams$q)/365
areaSearchedPerDay

## Erepro 

# Initial erepro is set using a general approximation now but will be optimised later. This way the steady state optimiser finds values easier
mariaParams$erepro <- 0.05 * mariaParams$w_inf^(-0.5)

#write.csv(mariaParams, file = "mariaParamsMar6.csv")

```

### Senescence and juvenile mort functions and params 

My initial explorations show that juvenile mortality function is quite important here, because increasing exponent (juv.e) from 0.3 to 1 (very steep and quick decrease in juv mortality) made some erepros crazy high (25K). Even values of 0.7 are too high and lead to Trachinus erepro of 6. 

```{r warning=T, message=FALSE, echo=T}

## Senescence

## function to set senescence mortality and add it on top of mu_b. The function needs mizer params object to be set up first to know the number of size classes. So I run MizerParams here, but can overwrite it later as long as no_w does not change

params <- MizerParams(mariaParams, interaction = inter, no_w = no_size_groups, kappa = kappa, lambda = lambda, w_pp_cutoff = w_pp_cutoff, r_pp = r_pp, min_w_pp = min_w_pp, kappa_ben = kappa_ben, lambda_ben = lambda_ben, w_bb_cutoff = w_bb_cutoff, r_bb = r_bb, min_w_bb = min_w_bb, kappa_alg = kappa_alg, lambda_alg = lambda_alg, w_aa_cutoff = w_aa_cutoff, r_aa = r_aa, min_w_aa = min_w_aa)

#parameters for senescence mortality as used in Law et al. 2009
k.sm <- 0.8 #mortality per year at the threshold size (should be 0.5 originally)
xsw <- 0.9 #proportion of w_inf at which mortality is at k.sm (should be 0.9)
sen.e <- 3 #exponent of the senescence mortality (larger value will give steeper increase in the last few sizes) (should be 3)

sen_mort = function(sppParams, params, k.sm, xsw, sen.e) {
  
  sen.mort.m = array(0, dim = c(length(sppParams$species),length(params@dw)))
  
  for (i in 1: length(sppParams$species)) {
    
    mu_Sen = k.sm * 10^(sen.e*(log(params@w)-log(xsw*sppParams$w_inf[i])))
    
    sen.mort.m[i,] <- mu_Sen    
  }

  # For really small species, like Trachinops predation mortality will be so high that senescence mortality is unlikely to be the case and perhaps should not even be applied
  sen.mort.m[which(params@species_params$w_inf < 100),] <- 0
    
  return(sen.mort.m)
  
}

## Juvenile mort 

## I also add extra juvenile mortality, because in reality background spectrum species would be imposing mortality on tiny fish, yet we don't model it here. As a result mortality curves look bell-shaped when instead they should be exponential 

juv.sm <- 10 #5 # mortality per year at threshold
juvsw <- 2 # Size of w_min at which mortality is at juv.sm. Note that at the moment w_min is setup at 10x the egg size. 
juv.e <- 0.3 #0.5 #exponent, larger values will give steeper decrease (orignal was 0.3, but it may be too slow decrease)

juv_mort = function(sppParams, params, juv.sm, juvsw, juv.e) {
  
  juv.mort.m = array(0, dim = c(length(sppParams$species),length(params@dw)))
  
  for (i in 1: length(sppParams$species)) {
    mu_Juv = juv.sm * 10^-(juv.e*(log(params@w) - log(juvsw*sppParams$w_min[i])))
    juv.mort.m[i,] <- mu_Juv
  }
  
  # It is most likely that for sharks or other species with large eggs (w_min of at least 1 gram), this early steep mortalit does not apply. So get rid of it
  juv.mort.m[which(params@species_params$w_min > 0.99),] <- 0
  
  return(juv.mort.m)
}

```

### Interaction 

```{r warning=T, message=FALSE, echo=T, eval = T}

availUr = 0.1 #urchins to fish
availUrLob = 0.7 #Urchins to lobsters
availSchooling = 0.35 #avaiability of schooling fish to predators

# and overwrite default interaction values with this (only for predators that feed on fish, i.e. have values > 0)
inter[c(which(inter[,which(mariaParams$species == "T_caudimaculatus")] >0)),which(mariaParams$species == "T_caudimaculatus")] <- availSchooling
inter[c(which(inter[,which(mariaParams$species == "P_laticlavius")] >0)),which(mariaParams$species == "P_laticlavius")] <- availSchooling
inter[c(which(inter[,which(mariaParams$species == "C_rasor")] >0)),which(mariaParams$species == "C_rasor")] <- availSchooling


inter[c(which(inter[,which(mariaParams$species == "urchins")] >0)),which(mariaParams$species == "urchins")] <- availUr
inter[which(mariaParams$species == "lobsters"), which(mariaParams$species == "urchins")] <- availUrLob


```

### Run   

```{r warning=FALSE, message=FALSE, echo=T}

### setup again with new erepro
params <- MizerParams(mariaParams, interaction = inter, no_w = no_size_groups, kappa = kappa, lambda = lambda, w_pp_cutoff = w_pp_cutoff, r_pp = r_pp, min_w_pp = min_w_pp, kappa_ben = kappa_ben, lambda_ben = lambda_ben, w_bb_cutoff = w_bb_cutoff, r_bb = r_bb, min_w_bb = min_w_bb, kappa_alg = kappa_alg, lambda_alg = lambda_alg, w_aa_cutoff = w_aa_cutoff, r_aa = r_aa, min_w_aa = min_w_aa)

#update initial abundances
params@initial_n <- params@initial_n/scalarForInitial

## add senescence and juvenile mortality to background mortality 
params@mu_b <- params@mu_b + sen_mort(mariaParams, params, k.sm, xsw, sen.e) + juv_mort(mariaParams, params, juv.sm, juvsw, juv.e)

## setup run time 
tmax = 100
tasm1 <- project(params, t_max = tmax, effort= 0, dt = dt, diet_steps = 0)
plot(tasm1)

```


### Get biomass and RDD/rmax

```{r warning=FALSE, message=FALSE, echo=T}

## total biomass at last step
getBiomass(tasm1)[tmax,]
#plotDietComp(tasm1)

#Check where they are at rmax level

model <- tasm1

a1 <- getRDI(params,model@n[tmax,,],model@n_pp[tmax,], model@n_bb[tmax,], model@n_aa[tmax,], model@intTempScalar[,,(tmax/dt)], model@metTempScalar[,,(tmax/dt)])

#get RDD 
a2 <- getRDD(params,model@n[tmax,,],model@n_pp[tmax,], model@n_bb[tmax,], model@n_aa[tmax,], sex_ratio = 0.5, model@intTempScalar[,,(tmax/dt)], model@metTempScalar[,,(tmax/dt)])

#get RDD to rmax ratio
rmaxratio <- a2/mariaParams$r_max
plot(rmaxratio)
#1-a2/a1


```

### Run with fishing   

```{r warning=FALSE, message=FALSE, echo=T}

#now run with some fishing mortality (equal to all)

tasm1ef <- project(params, t_max = tmax, effort=0.2, dt = 0.2, diet_steps = 0)  
plot(tasm1ef)

getBiomass(tasm1ef)[tmax,]/getBiomass(tasm1)[tmax,]
#plotDietComp(tasm1ef)

```

###  Calibration: functions

I will calibrate erepro, gamma, and 3 inter matrix terms, plus r_bb, r_aa, r_pp. 
Data to use for error functions: 

relative biomasses 
trends in biomasses
no extinction 
absolute biomasses 
diets 


```{r warning=FALSE, message=FALSE, echo=T}

#startpars is a list where each range of paramters get their value

startpars <- list()
startpars$gamma <- mariaParams$gamma
startpars$erepro <- mariaParams$erepro
startpars$inter <- c(availSchooling, availUr, availUrLob)
startpars$r <- c(r_pp, r_bb, r_aa)

length(startpars[['r']])

# this function will take the upper multiplier values and set upper bounds
upper_bounds <- function(startpars, inter.mult = inter.mult.par, erepro.mult = erepro.mult.par, gamma.mult = gamma.mult.par, r.mult = r.mult.par, no_sp = numb_of_spp) {
  
  #set an empty vector for upper values
  upper_value <- list()
  
  #if we are optimising erepro, we sample them on a log scale and make sure the upper bound is below 1 (or 0.5 or whatever upper value seems reasoable)
 if (length(startpars[['erepro']]) > 0) {
    startpars_erepro <- starpars$erepro
    upper_valueE <- startpars_erepro + log(erepro.mult) #we allow upper ereproboundary to be erepro.mult times bigger
    upper_valueE[c(which(exp(upper_valueE) > 0.2))] <- log(0.2) ## here we assume that erepro cannot be higher than 0.2
    upper_value$erepro <- upper_valueE
  }
  
 if(length(startpars[['gamma']]) > 0) {
    startpars_gamma <- startpars$gamma
    upper_valueG <- startpars_gamma * gamma.mult
    upper_valueG[c(which(upper_valueG > (max(startpars_gamma)*2)))] <- (max(startpars_gamma)*2) # set a upper limit of gamma at twice the maximum value of the total dataset (i.e. any species)
    upper_value$gamma <- upper_valueG
 }
  
   if (inter.opt) {
   startpars_inter <- startpars$inter 
   upper_valueI <- startpars_inter*inter.mult 
   upper_valueI[c(which(upper_valueI > 0.9))] <- 0.9 #here we assume that vulnerability cannot be higher than 0.9
   upper_value$inter <- upper_valueI
   }
  
     if (r.opt) {
   startpars_r <- startpars$r
   upper_valueR <- startpars_r*r.mult 
   upper_valueR[c(which(upper_valueR > 6))] <- 6 #I assume that regeneration rate cannot be higher than 6
   upper_value$r <- upper_valueR
  }
 
  return(upper_value)
}

lower_bounds <- function(startpars, inter.opt = T, erepro.opt = F, gamma.opt = F, inter.mult = inter.mult.par, erepro.mult = erepro.mult.par, gamma.mult = gamma.mult.par, no_sp = numb_of_spp) {
  
  lower_value <- as.numeric(vector())
  
  #if we are optimising erepro, we sample them on a log scale and make sure the upper bound is below 1 (or 0.5 or whatever upper value seems reasoable)
  if (erepro.opt) {
    startpars_erepro <- startpars[c(1:no_sp)]
    lower_valueE <- startpars_erepro - log(erepro.mult) #we allow upper ereproboundary to be erepro.mult times lower
    lower_valueE[c(which(exp(lower_valueE) < 1e-7))] <- log(1e-7) ## here we assume that erepro cannot be lower than 1e-7 (not sure this is neded...)
    lower_value <- c(lower_value, lower_valueE)
 }

    if(gamma.opt) {
    startpars_gamma <- startpars[c((no_sp+1):(no_sp*2))]
    lower_valueG <- startpars_gamma / gamma.mult
    lower_valueG[c(which(lower_valueG < (min(startpars_gamma)/2)))] <- (min(startpars_gamma)/2) # set a lower limit of gamma at half the minimum value of the total dataset 
    lower_value <- c(lower_value, lower_valueG)
    }
  
    if (inter.opt) {
    startpars_inter <- startpars[c((no_sp*2+1):(no_sp*2+4))]
    lower_valueI <- startpars_inter/inter.mult 
    lower_valueI[c(which(lower_valueI < 0.1))] <- 0.1 #here we assume that vulnerability cannot be lower than 0.1
    lower_value <- c(lower_value, lower_valueI)
  }

  
  return(lower_value)
}

#The optim function uses fn function to calculate errors
#This function will take parameter values from the optim function, will take my mizer object params to run the model (run_model function), and a few other thing

calibratePar_Tasm <- function(startpars, mariaParams, inter, meansteps= meansteps.par, effort = effort, trend.mult = trend.error.mult) {
  
# Function arguments: 
# startpars       - vectors of initial parameters for optimising (erepro, gamma  and interaction matrix values)
# mariaParams     - my species parameter file
# meansteps       - number of years to average biomass, at the end of the simulation run 

  optimizer_count <- optimizer_count + 1        #Keep runing count of function evaluations by optim()  
  assign("optimizer_count", optimizer_count, pos = .GlobalEnv)
  
  no_sp <- length(mariaParams$species)
  
### If calibrating erepro use this    
 #put optimised parameters into the param matrix
 startpars_erepro <- startpars[c(1:no_sp)]
 mariaParams$erepro <- exp(startpars_erepro) #for erepro we optimise the parameter on the log scale 
 
 #Then put the gamma values
 startpars_gamma <- startpars[c((no_sp+1):(no_sp*2))]
 mariaParams$gamma <- startpars_gamma
 
 ###now update the interaction matrix
 
 startpars_inter <- startpars[c((no_sp*2+1):(no_sp*2+4))]
  
 inter[c(which(inter[,which(mariaParams$species == "T_caudimaculatus")] >0)),which(mariaParams$species == "T_caudimaculatus")] <- startpars_inter[1]
inter[c(which(inter[,which(mariaParams$species == "P_laticlavius")] >0)),which(mariaParams$species == "P_laticlavius")] <- startpars_inter[2]
inter[c(which(inter[,which(mariaParams$species == "urchins")] >0)),which(mariaParams$species == "urchins")] <- startpars_inter[3]
inter[which(mariaParams$species == "lobsters"), which(mariaParams$species == "urchins")] <- startpars_inter[4]

  #setup mizerParams object with them
   params <- MizerParams(mariaParams, interaction = inter, no_w = no_size_groups, kappa = kappa, lambda = lambda, w_pp_cutoff = w_pp_cutoff, r_pp = r_pp, min_w_pp = min_w_pp, kappa_ben = kappa_ben, lambda_ben = lambda_ben, w_bb_cutoff = w_bb_cutoff, r_bb = r_bb, min_w_bb = min_w_bb, kappa_alg = kappa_alg, lambda_alg = lambda_alg, w_aa_cutoff = w_aa_cutoff, r_aa = r_aa, min_w_aa = min_w_aa)

## update initial abundances based on observations 
#for (i in 1:length(mariaParams$species)) {
#  params@initial_n[i,] <- params@initial_n[i,] *abund_scalar[i]
#} 
 
    #and add senescence and juvenile mortality to background mortality 
    params@mu_b <- params@mu_b + sen_mort(mariaParams, params, k.sm, xsw, sen.e) + juv_mort(mariaParams, params, juv.sm, juvsw, juv.e)

  #run the model for tmax years  
  tasm_model <- run_model(params, tmax, effort) ## run the model with the updated parameters 
  
  #compare relative biomasses 
  BiomError <- biom_error(tasm_model)
  print("Scaled biomass error")
  print(BiomError)
  
  #compare relative biomasses
  RelBioError <- relat_biom_error(tasm_model)
  print("Scaled relative biomass error")
  print(RelBioError*rel.bio.mult)  
  
  #penalise parameters that lead to long term trends 
  TrendError <- trend_error(tasm_model, tmax)
  print("Scaled trend error")
  print(TrendError*trend.mult)
  
  #Sum of all errors with relevant scalars 
  error <- BiomError + TrendError*trend.mult + RelBioError*rel.bio.mult #Sum error for biomass and trends (slopes). Since slopes are small I need to find the right multiplier to give it weight
  print("Total error")
  print(error)
  #Extinction penalty 
  #if(extinction_test){
  extinct <- Extinct_test(tasm_model)
  print("extinction error")
  print(extinct)
  
  error <- error + extinct
  #}
  
  print(paste(optimizer_count, "Error:", error))
  #print(exp(startpars)) #for erepro 
  print(startpars)
  return(error)
}

## run_model function

run_model <- function(params, tmax, effort) {
  tasm_model <- project(params, t_max = tmax, effort= effort, dt = 0.2, diet_steps = 0)
 return (tasm_model)

}

## extinction test function

Extinct_test <- function(tasm_model, extinct_threshold=extinct_threshold.par){

# This test is based on an initial biomass near what the target should be. 
# If the difference is more than 2 order of magnitude than the initial biomass, then there is a strong power penalty 
# Function arguments: 
# tasm_model          - result of project() 
# extinct_threshold   - If biomass decreases below this fraction of the initial biomass, the extintion penalty is applied 

    Biomass <- getBiomass(tasm_model) #returns tmax x species matrix
    BiomassInit<- Biomass[1,]  #biomass in the first year
    BiomassEnd<-  Biomass[dim(Biomass)[1],] #biomass in the last year 
    relB<-(BiomassEnd/BiomassInit)
    extinct.temp <- sum(1/relB[relB<extinct_threshold])
    
    return(extinct.temp)
}

## trend error function 
trend_error <- function(tasm_model, tmax) {
  Biomass <- getBiomass(tasm_model) #returns tmax x species matrix
  time <- c(burnin:tmax)  # use time from year 50 to tmax, to allow for the initial burn-in
  tr_err <- 0 #initialise trend error 
  for (i in 1:length(tasm_model@params@species_params$species)) {
     bio1 <- Biomass[c(burnin:tmax),i]
     mod <- lm(bio1 ~ time)
     slope <- as.numeric(mod$coefficients[2])
     tr_err <- tr_err + slope^2
     #print(tr_err)
  }
  
  return(tr_err)
}

## biomass error function

biom_error <- function(tasm_model, meansteps = meansteps.par){
  
  #vector of mean biomasses per survey. It ranges from 400 to 3500g 
  bio_obs <-  mariaParams$meanBio/(500*10) ## biomass per m2, so divide survey area of 500m2 and depth of 10 meters, so approximately 5000m3
  
  bio_model.t <- getBiomass(tasm_model) #years x species
  bio_model <- colMeans(bio_model.t[c((tmax - meansteps):tmax),])

  #Scale model biomasses to biomasses per survey - not sure this is needed as I scale initial abundances to biomasses per survey
  #bio_model.scaled <-bio_model * some_scalar #converts biomass from g/m3 to mtons / whole Eastern Bering Sea 

  #Get observed difference between observed and predicted biomass
  bdif<- log10(bio_obs)-log10(bio_model)
        
  b_error <-sum(na.omit(bdif^2)) # I don't think I need to scale in by variance  - Julia?
  
  #/ var(na.omit(log10(Bobs[Bobs!=0]))) #Standardize the RSS by the variance of the SSB (or biomass) observations. We're doing this because combining RSS from yeild and they on a smaller scale than biomas, but we choose to weight the data types equally (e.g., see Hilborn and Walters 1992).

  return(b_error)
}

## relative biomass error function
relat_biom_error <- function(tasm_model, meansteps = meansteps.par){
  
  #vector of mean biomasses per survey. It ranges from 400 to 3500g 
  bio_obs <-  mariaParams$meanBio/(500*10) ## biomass per m2, so divide survey area of 500m2 and depth of 10 meters, so approximately 5000m3
  rel_bio_obs <- bio_obs/(max(bio_obs))
  
  bio_model.t <- getBiomass(tasm_model) #years x species
  bio_model <- colMeans(bio_model.t[c((tmax - meansteps):tmax),])
  rel_bio_mod <- bio_model/(max(bio_model))
  
  #get the difference in relative abundances
  rel_error <- sum(abs(rel_bio_obs - rel_bio_mod))

  return(rel_error)
}

```


### Calibration: parameters

```{r warning=FALSE, message=FALSE, echo=T, eval = T}

tmax = 150 #how many years to run the model for each optimiser round. I found that 50 years is not enough, because species start going extinct after 100 years
effort = 0 #fishing mortality for calibrations
extinct_threshold.par = 0.01 #relative biomass at the end for a species to be extinct in the penalty function
burnin = 80 #how many initial years to remove for the trend error function - we only penalise trends in biomasses after initial burn-in period
trend.error.mult = 1e+10 # weight of the error function pelanlising a trend in biomass through time (we want stable biomasses). This has to be adjusted 
rel.bio.mult = 1000 # weight of the error function to compare relative biomasses. The maximum error can be 17, while the overall biomass error is ca 75 to 150. The relative biomasses are more important to get right, so this error has a weight of 100
inter.mult.par <- 2 #multiplier for the interaction paramter for upper and lower bounds - how many times do we allow the interaction to be higher/lower. The lower bound is later set at a minimum of 0.1 and the upper at a maximum of 0.8
erepro.mult.par <- 3 #multiplier for the erepro parameter for upper and lower bounds. This is on logarithmic scale. The upper bound is set at a maximum of 0.5 and lower at a a minimum of 1e-7
gamma.mult.par <- 2 #same as above but for gamma parameter. No limits for the upper bound, the lower bound is set a minimum of 1
meansteps.par <- 10 # how many final years to use for the mean species biomass to be compared with the observations
numb_of_spp <- length(mariaParams$species)

```

### Calibration: optim

```{r warning=FALSE, message=FALSE, echo=T, eval = F}
optimizer_count=0 # Initialize count of function evaluations

# starting parameter values
# let's say I start with exploring erepro in my 7 species model 
#startpars <- rep(0.001,times =length(mariaParams$species)) #initial erepro values as a test

#initial values from the steady state calculations - we do optim on the log scale
#startpars <- log(stedTas@species_params$erepro)

#how many params will I calibrate? 
#erepro for all species 
#gamma  for all species 
#interaction: avail of Trachinops, Pictilabrus and urchins

startpars_erepro <- log(mariaParams$erepro)
startpars_gamma <- mariaParams$gamma
startpars_inter <- c(availTr, availPi, availUr, availUrLob)
#startpars <- c(0.35, 0.6, 0.6)

startpars <- c(startpars_erepro, startpars_gamma, startpars_inter)

# Now describe the optim function that will do the search 

optim_result <-optimParallel(startpars,                   # starting parameter values
                     fn= calibratePar_Tasm,      # function to feed parameter values to the mizer object, run the model and return the error. It calls other functions from itself
                     mariaParams = mariaParams,       # pass the original parameter file (John passed his sim object after the burn-in period but I don't run burn-in for now)
                     inter = inter,
                     meansteps = meansteps.par,      # how many last years to use for the mean biomass calculations
                     #extinction_test=TRUE,      # I use it    
  lower = lower_bounds(startpars, inter.opt = T, erepro.opt = T, gamma.opt = T, inter.mult = inter.mult.par, erepro.mult = erepro.mult.par, gamma.mult = gamma.mult.par),
  upper = upper_bounds(startpars, inter.opt = T, erepro.opt = T, gamma.opt = T, inter.mult = inter.mult.par, erepro.mult = erepro.mult.par, gamma.mult = gamma.mult.par),
                     method="L-BFGS-B",
                     control=list(maxit=10, REPORT=1, trace=6), 
                     effort = effort, 
  parallel = list(loginfo = TRUE, forward = TRUE))

library("optimParallel")
library("parallel")
cl <- makeCluster(detectCores())
cl <- makeCluster(6)

setDefaultCluster(cl = cl)
clusterEvalQ(cl, library(mizer))

#in optim function have parallel = list(loginfo = TRUE, forward = TRUE)
#optim hessian=TRUE will give you parameter correlations and 

#For erepro calibration you want to print this
#exp(optim_result$par)/exp(startpars)

#For interaction calibration print this 
optim_result$par

#save(optim_result, file = "optim_resultFeb12.RData")

```




### Pred kernel experiments 

```{r}
tmax = 100
dt = 0.2
tasm1 <- project(params, t_max = tmax, effort= 0, dt = dt, diet_steps = 0)
plot(tasm1)

params2 <- change_pred_kernel(params, params@pred_kernel)
tasm2 <- project(params2, t_max = tmax, effort= 0, dt = dt, diet_steps = 0)
plot(tasm2)

```

### Net intake through time (e)

```{r warning=FALSE, message=FALSE, echo=T}

model <- tasm1

#get the net intake values (starvation mortality is e/(0.1*w))

intake_matrix <- array(NA, dim = c(length(params@species_params$species), no_size_groups, tmax))

for (i in 1:tmax) {
  timetoget <- i

e_values <- getEReproAndGrowth(params,model@n[timetoget,,],model@n_pp[timetoget,], model@n_bb[timetoget,], model@n_aa[timetoget,], model@intTempScalar[,,(timetoget/dt)], model@metTempScalar[,,(timetoget/dt)])

intake_matrix[,,i] <- e_values

}

timetoplot = seq(from = 1, to = tmax, by = 20)

#create a colour gradient to plot
colfunc <- colorRampPalette(c("black", "lightgrey"))
mycol <- colfunc(length(timetoplot))

for (i in 1:length(mariaParams$species)) {
  
  first_slot <- max(which(params@w < mariaParams$w_min[i])) #only plot relevant size values
  last_slot <- min(which(params@w > mariaParams$w_inf[i]))
  
  plot(log(params@w[c(first_slot:last_slot)]),intake_matrix[i,c(first_slot:last_slot),1], type = 'l', xlab = 'Log size, g', ylab = 'net intake, g', main = mariaParams$species[i])
      for (j in 1:length(timetoplot)) {
        points(log(params@w[c(first_slot:last_slot)]),intake_matrix[i,c(first_slot:last_slot),j],type = 'l', col = mycol[j])
      }
  abline(v = log(mariaParams$w_mat[i]), lty = 2, col = 'orange')
  abline(h = 0, lty = 1, col = 'red')
  
}

##
starvMort <- (intake_matrix)/(0.1*params@w)
starvMort[intake_matrix>0] <- 0
starvMort <- -starvMort


```

### Background depletion 

```{r warning=FALSE, message=FALSE, echo=T}

tasm1 <- tasm1

smallcut <- max(which(tasm1@params@w_full < min_w_bb)) +1
largecut <- min(which(tasm1@params@w_full > w_bb_cutoff)) -1

largecutP <- min(which(tasm1@params@w_full > w_pp_cutoff)) -1

smallcutA <- max(which(tasm1@params@w_full < min_w_aa)) +1
largecutA <- min(which(tasm1@params@w_full > w_aa_cutoff)) -1

plot(log(tasm1@params@w_full[c(smallcut:largecut)]), (tasm1@n_bb[tmax,c(smallcut:largecut)]/tasm1@n_bb[1,c(smallcut:largecut)]),ylim = c(0,1), xlim = c( log(tasm1@params@w_full[1]), log(tasm1@params@w_full[largecut])), type = 'l', col = 'red', xlab = "log weight,g", ylab = "depletion compared to initial abundance", main = "Depletion of plankton, benthos(red) and algae(green)")
points(log(tasm1@params@w_full[c(1:largecutP)]), (tasm1@n_pp[tmax,c(1:largecutP)]/tasm1@n_pp[1,c(1:largecutP)]),type = 'l', col = 'black')
points(log(tasm1@params@w_full[c(smallcutA:largecutA)]), (tasm1@n_aa[tmax,c(smallcutA:largecutA)]/tasm1@n_aa[1,c(smallcutA:largecutA)]),type = 'l', col = 'green')


```

### Growth curves    
Plot and check individual growth curves. Note that the reference curves are based on vb_k and Linf, and vb_k is a highly uncertain parameter for many species. So black refrence curves are not good in many species. 

```{r warning=FALSE, message=FALSE, echo=T}
tasm1 <- tasm1
## growth curves
for (i in 1: length(params@species_params$species)) {
  
  plotGrowthCurves(tasm1, species = params@species_params$species[i])
  #plotGrowthCurves(tasm1temp, species = params@species_params$species[i])
  #plotFeedingLevel(tasm1, species = params@species_params$species[i])

}
```

#Plot all sizespectra at a chosen time

```{r warning=FALSE, message=FALSE, echo=T}
tplot =2

model <- tasm1

#plot size spectra
plot(log10(model@params@w_full), log10(model@n_pp[tplot,]*params@dw_full), type = 'l', ylim = c(-10,10), xlim = c(-5,5),lwd = 2, xlab = 'size, g', ylab = 'abundance')
points(log10(model@params@w_full), log10(model@n_bb[tplot,]*params@dw_full), type = 'l', lwd = 2, col = 'red')
points(log10(model@params@w_full), log10(model@n_aa[tplot,]*params@dw_full), type = 'l', lwd = 1, col = 'green')

for (i in 1: length(params@species_params$species)) {
  points(log10(model@params@w), log10(model@n[tplot,i,]*params@dw), type = 'l')
  
}

ppatmat = log(mariaParams$w_mat/mariaParams$beta) # prefered prey size at predators maturation size
abline(v = c(ppatmat), lty = 2, col = 'grey')

abline(v = log(mariaParams$w_mat), lty = 2, col = 'yellow')

```

### Run again with optim

```{r warning=FALSE, message=FALSE, echo=T}

#put optimised parameters into the param matrix
 mariaParams$erepro <- exp(optim_result$par[c(1:no_sp)])

 #Then put the gamma values
 mariaParams$gamma <- optim_result$par[c((no_sp+1):(no_sp*2))]

 ###now update the interaction matrix
 interPars <- optim_result$par[c((no_sp*2+1):(no_sp*2+4))]
  
 inter[c(which(inter[,which(mariaParams$species == "T_caudimaculatus")] >0)),which(mariaParams$species == "T_caudimaculatus")] <- interPars[1]
inter[c(which(inter[,which(mariaParams$species == "P_laticlavius")] >0)),which(mariaParams$species == "P_laticlavius")] <- interPars[2]
inter[c(which(inter[,which(mariaParams$species == "urchins")] >0)),which(mariaParams$species == "urchins")] <- interPars[3]
inter[which(mariaParams$species == "lobsters"), which(mariaParams$species == "urchins")] <- interPars[4]



### setup again with new erepro
params <- MizerParams(mariaParams, interaction = inter, no_w = no_size_groups, kappa = kappa, lambda = lambda, w_pp_cutoff = w_pp_cutoff, r_pp = r_pp, min_w_pp = min_w_pp, kappa_ben = kappa_ben, lambda_ben = lambda_ben, w_bb_cutoff = w_bb_cutoff, r_bb = r_bb, min_w_bb = min_w_bb, kappa_alg = kappa_alg, lambda_alg = lambda_alg, w_aa_cutoff = w_aa_cutoff, r_aa = r_aa, min_w_aa = min_w_aa)

## again update initial abundances based on observations 
#for (i in 1:length(mariaParams$species)) {
# get the name of the species for indexing 
#  sppname <- as.character(mariaParams$species[i])
  # multiply initial abundances by the relative abundance of this species in the observations. Also multiply it by the total_scalar applied to all species to get the abundance closer to the equilibrium conditons 
  #params@initial_n[i,] <- params@initial_n[i,] * observ_abund$scaledAbund[which(observ_abund$species == sppname)] * total_scalar
#  params@initial_n[i,] <- params@initial_n[i,] *sc[i]
#} 

## add senescence and juvenile mortality to background mortality 
params@mu_b <- params@mu_b + sen_mort(mariaParams, params, k.sm, xsw, sen.e) + juv_mort(mariaParams, params, juv.sm, juvsw, juv.e)

#params@species_params$erepro[which(params@species_params$species == "D_lewini")] <- 0.024
#params@species_params$erepro[which(params@species_params$species == "P_bachus")] <- 0.0099

## setup run time 
tmax = 250
dt = 0.2
tasm1 <- project(params, t_max = tmax, effort= 0, dt = dt, diet_steps = 1)
plot(tasm1)

```











### Mortality - FIX

The graph is not perfect, but a quick look at the assumed juvenile and senescence mortality functions and how they compare with the predation mortality. Minimum size is shown with dashed blue line, maturation size is shown with orange lines. Basically we see that juvenile mortality is very high at minimum sizes and drops quickly at around sizes 0.1g, by which time predation becomes a more important force. Note, no juvenile mortality is imposed on sharks, for which min size is assumed to be 1g. Senescence mortality kicks in just before the 90% of maximum size. The curves look steep, but the maximum values are around 1 at maximum size. 

```{r warning=FALSE, message=FALSE, echo=T}

#at which time to plot 
year = 50

m2mort <- getM2(params,model@n[year,,],model@n_pp[year,], model@n_bb[year,], model@n_aa[year,], model@intTempScalar[,,(year/dt)])

getPredMort <- function(object, n, n_pp, n_bb, n_aa, pred_rate, intakeScalar, time_range, drop = TRUE

# get predation mortality 
m2mort <- getM2(model, model@intTempScalar[,,(year/dt)], )[year,,]
#get non-predation mortality 
nonm2mort <- params@mu_b

#plot non predation mort
plot(log(params@dw), nonm2mort[1,], type = 'l', ylim = c(0, 5), col = 'grey', xlab = "log size, g", ylab = "instantaneous mortality", main = "Non-predation and predation (red) mortality")
for (i in 2:length(mariaParams$species)) {
  points(log(params@dw), nonm2mort[i,], type = 'l')  
}
#add predation 
for (i in 1:length(mariaParams$species)) {
  points(log(params@dw), m2mort[i,], type = 'l', col = 'red')  
}
#who w_min
abline(v = log(mariaParams$w_min), lty = 2, col  = 'blue')
abline(v = log(mariaParams$w_mat), lty = 2, col  = 'orange')

#totmort <- getM2(tasm1)[tmax,,] + params@mu_b

#or just explore non predation mortality
#plot(x = log(tasm1@params@w), tasm1@params@mu_b[1,], ylim = c(0,50), type = 'l')
#for (i in 1: length(params@species_params$species)) {
#  points(log(tasm1@params@w), tasm1@params@mu_b[i,], type = 'l')
#}
#abline(v = log(mariaParams$w_min), lty = 2, col  = 'red')
#abline(v = log(mariaParams$w_inf), lty = 2, col  = 'blue')


```

### Size spectra  

Initial and final abundance at size: still too many individuals at largest sizes

```{r warning=FALSE, message=FALSE, echo=T}
#plot initial and final abundance 
par(mfrow = c(1,2))

mycol = c('black','red', 'green', 'orange','blue','violet','pink','grey','yellow','black','red', 'green', 'orange','blue','violet','pink','grey','yellow')

#initial
plot(log(params@w), log(params@initial_n[1,]), type = 'l', ylim = c(-40, max(log(params@initial_n[]))), main = "initial n", xlab = "log w", ylab = "numbers")

for (i in 2:length(params@species_params$species)) {
  points(log(params@w), log(params@initial_n[i,]), type = 'l', col = mycol [i] )
}
abline (h = c(0,-10,-20,-30), lty = 2, col = 'grey')
abline (v = c(-5,0,5), lty = 2, col = 'grey')

#final
plot(log(params@w), log(tasm1ef@n[tmax,1,]), type = 'l', ylim = c(-40, max(log(params@initial_n[]))), main = "final n", xlab = "log w", ylab = "numbers")

for (i in 2:length(params@species_params$species)) {
  points(log(params@w), log(tasm1ef@n[tmax,i,]), type = 'l', col = mycol [i] )
}
abline (h = c(0,-10,-20,-30), lty = 2, col = 'grey')
abline (v = c(-5,0,5), lty = 2, col = 'grey')
par(mfrow = c(1,1))

```


### Plot initial densities 
```{r}

plot(log(params@w_full), log(params@initial_n_pp), type = 'l', xlim = c(-10,10), ylim = c(-30,30), col = 'grey')
points(log(params@w_full), log(params@initial_n_aa), type = 'l', col = 'green')
points(log(params@w_full), log(params@initial_n_bb), type = 'l', col = 'red')

for (i in 1: length(params@species_params$species)) {
  points(log(params@w), log(params@initial_n[i,]), type = 'l')
  
} 

# default initial abundances are quite far from the observed ones 

#get the minimum observed weight size group 
minobs <- min(which(params@w > 0.1))

#get mean model abundances in the observed weight groups (so they are comparable to surveys)
modelAb <- rep(NA, length(params@species_params$species))
## derive initial abundance scalar 
for (i in 1:length(mariaParams$species)) {
modelAb[i] <- (sum(params@initial_n[i,c(minobs:200)]))
}

#get a scalar between observed abundances and model abundances 
sc <- (params@species_params$meanAb/500)/modelAb

#plot them again, this time species abundances adjusted
plot(log(params@w_full), log(params@initial_n_pp), type = 'l', xlim = c(-10,10), ylim = c(-30,10), col = 'grey')
points(log(params@w_full), log(params@initial_n_aa), type = 'l', col = 'green')
points(log(params@w_full), log(params@initial_n_bb), type = 'l', col = 'red')

for (i in 1: length(params@species_params$species)) {
  points(log(params@w), log(params@initial_n[i,]*sc[i]), type = 'l')
  
} 



plot(params@w_full, params@initial_n_pp, type = 'l', xlim = c(0.1,20), ylim = c(0,200), xlab = 'weight,g', ylab = 'numbers')
points((params@w_full), (params@initial_n_aa), type = 'l', col = 'green')
points((params@w_full), (params@initial_n_bb), type = 'l', col = 'red')

biom_pp <- params@w_full * params@initial_n_pp
biom_aa <- params@w_full * params@initial_n_aa
biom_bb <- params@w_full * params@initial_n_bb

plot(params@w_full, biom_pp, type = 'l', xlim = c(0.1,20), ylim = c(0,50), xlab = 'weight,g', ylab = 'biomass,g')
points((params@w_full), biom_aa, type = 'l', col = 'green')
points((params@w_full), biom_bb, type = 'l', col = 'red')

sum(biom_bb)
sum(biom_aa)
sum(biom_pp)

```

### DATA CODES

### Code to get abundance per suvery 

Getting data for urchins, lobsters and fish. This is setup to extract either average or time series mean abundance and mean Biomass data. Currently using only data before year 2000.

```{r warning=T, message=FALSE, echo=T, eval=F}
load(file="C:/Users/astaa/Dropbox/asta/RLS/DATAsets/rlsbas.RData")
library(tidyverse)

MainLocs =c("Tinderbox Reserve", "Tinderbox External", "Maria Island Reserve", "Maria Island Vincinity", "Maria External", "Ninepin Internal","Ninepin External", "Bicheno Internal", "Bicheno External")

urchinSum <- rls_bas %>% 
  filter(Location %in% MainLocs) %>%
  filter (year < 2000) %>% 
    filter(Method==2) %>% 
      filter( GENUS == "Heliocidaris" | GENUS == "Centrostephanus" | GENUS == "Goniocidaris" ) %>%
     #   group_by (year) %>%
          summarise(meanAb = mean(Abundance)*5, sdAb = sd(Abundance)*5, maxAb = max(Abundance)*5, meanBio = meanAb * 5) #to convert to biomass, assume mean weight of ca1g (5cm across, 0.01*5^3). For biomass I assume that an average urchin weighs 5 grams

lobsterSum <- rls_bas %>% 
  filter(Location %in% MainLocs) %>%
  filter (year < 2000) %>%
    filter(Method==2) %>% 
      filter( GENUS == "Jasus") %>%
      #  group_by (year) %>%
          summarise(meanAb = mean(Abundance)*5, sdAb = sd(Abundance)*5, maxAb = max(Abundance)*5, meanBio = meanAb*300) #to convert to biomass, assume mean weight of 150g
#For method 2 remember to multiply everything by 5 to get comparable abundances. For biomass I assume that average observed lobster weighs 300g (below the maturaiton size)

## now get observed fish abundances for the required period 

ModelSpp = c("Caesioperca rasor", "Cheilodactylus spectabilis" ,"Dinolestes lewini" ,"Latridopsis forsteri" ,"Aplodactylus arctidens" ,"Notolabrus tetricus", "Notolabrus fucicola" ,"Cephaloscyllium laticeps", "Meuschenia australis" ,"Acanthaluteres vittiger" ,"Pseudophycis bachus" ,"Trachinops caudimaculatus", "Olisthops cyanomelas" ,"Pentaceropsis recurvirostris", "Meuschenia freycineti" ,"Pictilabrus laticlavius")

MatchNames <- list()
MatchNames$modelName <- c("C_rasor", "C_spectabilis" ,"D_lewini" ,"L_forsteri" ,"A_arctidens" ,"Notolabrus", "Notolabrus" ,"C_laticeps", "leatherjack" ,"leatherjack" ,"P_bachus" ,"T_caudimaculatus", "O_cyanomelas" ,"Boarfish", "M_freycineti" ,"P_laticlavius")
MatchNames$sppName <- ModelSpp

fishSumm <- rls_bas %>% 
  filter (Location %in% MainLocs) %>% 
  filter (year < 2000) %>% 
    filter (Method == 1) %>%
      filter(TAXONOMIC_NAME %in% ModelSpp) %>% 
        group_by(TAXONOMIC_NAME) %>%
          summarise(meanAb = mean(Abundance), sdAb = sd(Abundance), maxAb = max(Abundance), meanBio = mean(BioMass), sdBio = sd(BioMass), maxBio = max(BioMass))

fishSumm$modelName <- MatchNames$modelName[match(fishSumm$TAXONOMIC_NAME, MatchNames$sppName)]

fishAbund <- fishSumm %>% group_by(modelName) %>% summarise(meanAb = mean(meanAb), sdAb = mean(sdAb), maxAb = mean(maxAb), meanBio = mean(meanBio), sdBio = mean(sdBio), maxBio = mean(maxBio))

#Add urchin observations
fishAbund <- add_row(fishAbund)
fishAbund[length(fishAbund$modelName),c(2:5)] <- urchinSum
fishAbund[length(fishAbund$modelName),1] <- "urchins"

#Add lobster observations
fishAbund <- add_row(fishAbund)
fishAbund[length(fishAbund$modelName),c(2:5)] <- lobsterSum
fishAbund[length(fishAbund$modelName),1] <- "lobsters"

#for Trachinops mean abundance is 110 but it gives diffcult calibration so I set to lower number for now (note that here we only record abundance above 2cm, whereas model counts all abundance and assumes same size slopes for all species)
fishAbund$meanAb[which(fishAbund$modelName == "T_caudimaculatus")] <- 60
fishAbund$meanAb[which(fishAbund$modelName == "urchins")] <- 100 #instead of 230, which drives the system nuts I think

save(fishAbund, file = "meanAbundBefore2000.RData")

#For lobsters LW conversion is as 
#W = 0.000285*L^3.114 males
#W = 0.000271*L^3.135 females
#to convert from carapace length in mm to grams

```

### Derive benthic spectra 

In this section I use data from very small benthic invertebrates collected by Kate Fraser (PhD student), and standardised per m2. For invertebrates bigger than 2g, data from RLS is used and compiled by Freddie Heather into weight bins on log2 scale. The data is combined, binned into equal weight bins on log10 scale and then a regression slope is fitted. 
From the first regression and by plotting abundances we can see that there is a big increase in abundance at largest weight groups, which is mostly due to urchins and partly lobsters. Given that urchins and lobsters are modelled explicitly we want to exclude them from teh background spectrum. Ideally this should be done by reanalysing Freddie's data without urchins and lobsters, but for now I just test a regression with a steeper sloe. Assuming a steeper slope and actually gives a better fit to small weight groups. We have to remember that Kate's data can underestimate the abundance (critters escape), but is less likely to overestimate it


```{r, eval = F}

library(tidyverse)

Freddie <- read.csv(file = "C:/Users/astaa/Dropbox/asta/TasmanianModel/Parameters/BenthicInvertData/Freddie.csv") #data of large invertebrates

Kate <- read.csv(file = "C:/Users/astaa/Dropbox/asta/TasmanianModel/Parameters/BenthicInvertData/all_bysample_size.csv") #small invertebrates
KateCodes <- read.csv(file = "C:/Users/astaa/Dropbox/asta/TasmanianModel/Parameters/BenthicInvertData/sample_data_env_vectors.csv") #site codes for Kate's data

# coordinates for eastern Tasmanian coast
#TasmCoorLon <- c(146.5, 148.5) 
#TasmCoorLat <- c(-41.00, -43.60)
KateWgtGroups <- c(1.8E-06,4.8E-06, 1.1E-05, 2.9E-05, 7.1E-05, 1.8E-04, 4.5E-04, 1.1E-03, 2.8E-03, 6.8E-03, 1.7E-02, 4.2E-02, 1.1E-01, 2.5E-01, 6.8E-01, 1.6E+00) #weight groups used in Kate's data (size of sieves)

TasmCodes <- KateCodes %>% filter (lat > -43.60 & lat < -41.00) %>% filter (lon > 146.5 & lon < 148.5)
TasmCodes.l <- unique(TasmCodes$SampleCode)
TasmInvSmall <- Kate %>% filter (SampleCode %in% TasmCodes.l) 
TasmInvSmall <- TasmInvSmall[,-c(1:2)]

MeanAbund <- apply(TasmInvSmall, 2, mean)
MinAbund <- apply(TasmInvSmall, 2, min)
MaxAbund <- apply(TasmInvSmall, 2, max)

InvData <- list()
InvData$wgt <- KateWgtGroups
InvData$meanAb <- MeanAbund
InvData$minAb <- MinAbund
InvData$maxAb <- MaxAbund
InvData <- as.data.frame(InvData)
InvData <- InvData[-c(which(InvData$meanAb ==0)),] # remove size groups that have 0 abundance

#combine with Freddie's data about average abundance at size in Tasmania
InvData <- rbind(InvData, Freddie)
#save(InvData, file = "data/BenticInvAbund_Tasm.RData")

load(file = "data/BenticInvAbund_Tasm.RData")
## now I have to bin data into bins of equal size on log10 scale 
InvData$cat <- NA
log10(min(InvData$wgt)) # minimum size on log10 scale
log10(max(InvData$wgt))

breaks <- seq(from = -6, to = 4) #make a vector of breaks on a log10 scale

#for each empirical weight put assign the bin on the log10 scale
for (i in 1:length(breaks)) {
  temp <- which((log10(InvData$wgt) > breaks[i]) & (log10(InvData$wgt) < breaks[i+1]))
  InvData$cat[c(temp)] <- i
  
}
#get abundances for the weight groups on the equal log10 scale (normalised)
InvDataBinned <- InvData %>% group_by(cat) %>% summarise (mean_ab = sum(meanAb), min_ab = sum(minAb), max_ab = sum(maxAb), meanWgt = mean(wgt))

#add the midpoint weight for each of the log10 weight groups 
InvDataBinned$wgtBinLog10 <- seq(from = -5.5, to = 3.5, by = 1) #mean weight in the bin on log 10 scale
InvDataBinned$wgtGroup <- 10^(InvDataBinned$wgtBinLog10)
InvDataBinned$bioPerm2 <- InvDataBinned$wgtGroup * InvDataBinned$mean_ab
#save(InvDataBinned, file = "data/InvDataBinned.RData")

#Now we fit linear regression to the total abundance data 
bs <- lm(log10(InvDataBinned$mean_ab) ~ log10(InvDataBinned$wgtGroup))
summary(bs)

#Plot it 
plot(log10(InvDataBinned$wgtGroup), log10(InvDataBinned$mean_ab), type = 'l', ylim = c(-4, 6), main = "Benthos slopes in Tasmania: log10(Abund) = 1.1 - 0.62*log10(wgt_bin)", xlab = "Log10 (w,g)", ylab = "Log10 (Abundance: min, mean, max)")
points(log10(InvDataBinned$wgtGroup), log10(InvDataBinned$min_ab), type = 'l', lty = 2)
points(log10(InvDataBinned$wgtGroup), log10(InvDataBinned$max_ab), type = 'l', lty = 2)
abline(lm(log10(InvDataBinned$mean_ab) ~ log10(InvDataBinned$wgtGroup)), col = 'red')
points(log10(InvDataBinned$wgtGroup), AbNoUr, col = 'blue', type = 'l')
abline(v = log10(2), lty = 3) #this is where Freddie's data starts and Kate's data ends

## We can see that there is a big increase in abundance at largest weight groups, which is mostly due to urchins and partly lobsters. Given that urchins and lobsters are modelled explicitly we want to exclude them from teh background spectrum. Assuming a steeper slope and actually gives a better fit to small weight groups 
AbNoUr <- 0.8 - 0.85*log10(InvDataBinned$wgtGroup)
## this equation is now used a very general approximation of the benthic slope

#or if just plotting the original (nonbinned) weights 
plot(log10(InvData$wgt), log10(InvData$meanAb), type = 'l', ylim = c(-4, 6), main = "Benthos slopes in Tasmania: not normalised weight groups", xlab = "Log10 (w,g)", ylab = "Log10 (Abundance: min, mean, max)")
points(log10(InvData$wgt), log10(InvData$minAb), type = 'l', lty = 2)
points(log10(InvData$wgt), log10(InvData$maxAb), type = 'l', lty = 2)
abline(lm(log10(InvData$meanAb) ~ log10(InvData$wgt)), col = 'red')
abline(v = log10(2), lty = 3)

```


### Calibrate benthic spectrum abundance

Mizer does not have defined units and can be calibrated for m3, m2 or the total model area. I am using measures per m2, because this is how RLS data is presented. 
AbNoUr <- 0.8 - 0.85*log10(InvDataBinned$wgtGroup)

```{r, eval = F}
#inferred observed abundance (on log10 scale from Kate's and Freddie's data)
#Data using the estimated regression slope
obsab <- 0.8 - 0.85*log10(params@w_full) ## expected abundance using the equation applied to the weight groups used in mizer

#length(params@dw_full[c(min(which(params@w_full > min_w_bb)):max(which(params@w_full < w_bb_cutoff)))])
modbb <- log10(params@initial_n_bb * params@dw_full) 
modpl <- log10(params@initial_n_pp * params@dw_full)
modaa <- log10(params@initial_n_aa * params@dw_full)

plot(log10(params@w_full), obsab, type = 'l', lwd = 1.5, xlim = c(-6,2), ylim = c(-3, 5), main = "Slopes of initial benthic abundance versus observed abundance", xlab = "Log10, w, g", ylab = "Log10(Abundance per m2")
points(log10(params@w_full), modbb, type = 'l', col = 'red')
points(log10(params@w_full),modpl, type = 'l', col = 'green')
points(log10(params@w_full),modaa, type = 'l', col = 'blue')
#abline(v = log10(0.05), lty = 2)
abline(v = log10(min_w_bb), lty =2)

```

Comment on the above: there is a big dip in the observed abundance data at around the size of 1g. This could be because the two methods used are just diverging on their sampling accuracy at this size or because at this size there is indeed a strong benthos depletion and a dip in the slope. I actually do get a dip at this point. The theoretical abundance at 1g would be 10individuals per m2, while we get less than 1 in observations. 


### Get initial abundance scalar 

I want to calibrate initial model abundances with the relative abundances observed in the surveys (in the first 5 years), so that I at least start in the correct ballpark for relative abundances. For this I only use model abundances that are likely to be in the range observable by RLS divers, e.g. bigger than 4cm (the smallest size group is 2cm, but few individuals are observed in this group), or very roughly 1g

```{r}

temp <- params@initial_n * params@dw * params@w #biomass 
initBiom <- apply(temp, 1, sum) # total biomass
initBiom[1] #wrasse
initBiom[3] #trachinops

temp_ben <- params@initial_n_bb * params@dw_full * params@w_full #biomass of benthos
initBiomBen <- sum(temp_ben)
initBiomBen

#load(file = "data/InvDataBinned.RData")
obsBenBio <- sum(InvDataBinned$bioPerm2[c(1:7)])
obsBenBio <- 19

initBiomBen/obsBenBio  #initial benthos biomass is about 40 times bigger than observed

initBiom[1]/(mariaParams$meanBio[1]/500) #initial biomass of wrasse 1700 times bigger than observed
initBiom[3]/(mariaParams$meanBio[3]/500) #for trachinops it is 17K times bigger than observed

modelBiom <- mariaParams$meanBio/500

#Now this is a scalar to modify initial densities while accounting for the match of observed vs modelled initial benthos densities. I will have to divide initial densities by this scalar
scalarForInitial <- (initBiom/modelBiom)/(initBiomBen/obsBenBio)
tt <- round(as.vector(scalarForInitial),1)

scalarForInit <-c(38.5 ,  0.3, 403.4,  10.6 ,  5.0  , 3.4 ,  0.2 ,  0.01, 114.5 ,  4.9  , 1.5 , 15.9,  10.3 ,  2.5 ,  0.4  , 7.1  , 0.4)


```


### Get approximate rmax value

```{r warning=FALSE, message=FALSE, echo=T}
mariaParams$r_max <- Inf

### setup again with new erepro
params <- MizerParams(mariaParams, interaction = inter, no_w = no_size_groups, kappa = kappa, lambda = lambda, w_pp_cutoff = w_pp_cutoff, r_pp = r_pp, min_w_pp = min_w_pp, kappa_ben = kappa_ben, lambda_ben = lambda_ben, w_bb_cutoff = w_bb_cutoff, r_bb = r_bb, min_w_bb = min_w_bb, kappa_alg = kappa_alg, lambda_alg = lambda_alg, w_aa_cutoff = w_aa_cutoff, r_aa = r_aa, min_w_aa = min_w_aa)

#update initial abundances
params@initial_n <- params@initial_n/scalarForInitial

## add senescence and juvenile mortality to background mortality 
params@mu_b <- params@mu_b + sen_mort(mariaParams, params, k.sm, xsw, sen.e) + juv_mort(mariaParams, params, juv.sm, juvsw, juv.e)

#params@species_params$erepro[which(params@species_params$species == "D_lewini")] <- 0.024
#params@species_params$erepro[which(params@species_params$species == "P_bachus")] <- 0.0099

## setup run time 
tmax = 3
tasm1 <- project(params, t_max = tmax, effort= 0, dt = dt, diet_steps = 0)
plot(tasm1)

#At around the time of year 5 biomasses are still fairly stable, so I will get RDI from this time 
model <- tasm1

#this is rmax independent recruitment with the default erepro values and when rmax was set to infinity 
rdi <- getRDI(params,model@n[tmax,,],model@n_pp[tmax,], model@n_bb[tmax,], model@n_aa[tmax,], model@intTempScalar[,,(tmax/dt)], model@metTempScalar[,,(tmax/dt)])

#check the values
plot(mariaParams$w_inf, rdi)
mean(rdi)
sd(rdi)
mean(rdi)/sd(rdi)
median(rdi)

#We can see that the variation is not too huge, and at this time step the size spectra look reasonable, so there is sufficient recruitment happening. Let's set the rmax at 5 times this recruitment value, which will then set a cap on the maximum recruitment 

#now lets set rmax at 10 times this values 
mariaParams$r_max <- rdi *5

write.csv(mariaParams, file = "Tasm_parametersMar7.csv")

#Recruitment values at tmax = 3 were
#c(0.0000238 0.0012437 0.0000137 0.0001259 0.0002768 0.0001776 0.0007411 0.0003186 0.0000155 0.0001562 0.0002777 0.0000693 0.0000103 0.0000334 0.0012713 0.0000530 0.0018097)

```


### Steady state to get erepro

This is not currently done because I will start with the general equation and optimise is with optim

```{r warning=T, message=FALSE, echo=T, eval = F}

### setup  
params <- MizerParams(mariaParams, interaction = inter, no_w = no_size_groups, kappa = kappa, lambda = lambda, w_pp_cutoff = w_pp_cutoff, r_pp = r_pp, min_w_pp = min_w_pp, kappa_ben = kappa_ben, lambda_ben = lambda_ben, w_bb_cutoff = w_bb_cutoff, r_bb = r_bb, min_w_bb = min_w_bb, kappa_alg = kappa_alg, lambda_alg = lambda_alg, w_aa_cutoff = w_aa_cutoff, r_aa = r_aa, min_w_aa = min_w_aa)

#update initial abundances
params@initial_n <- params@initial_n/scalarForInitial

## add senescence and juvenile mortality to background mortality 
params@mu_b <- params@mu_b + sen_mort(mariaParams, params, k.sm, xsw, sen.e) + juv_mort(mariaParams, params, juv.sm, juvsw, juv.e)


## find a steady state by modifying erepro
stedTas <- steady(params, t_max = 300) #t_max gives the number of interations needed 

#plot erepro values and check if they are realistic 
plot(log(stedTas@species_params$w_inf), stedTas@species_params$erepro, xlab = 'log w_inf', ylab = 'erepro', pch = 19)

#just output the values
stedTas@species_params$erepro
params@species_params$species

#to do this iteratively put, the new erepro values into params and run steady again. Might need to do a few iterations
params@species_params$erepro <- stedTas@species_params$erepro

#now run steady function again
stedTas <- steady(params, t_max = 300) 

#for which species erepro is > 1
mariaParams$species[c(which(stedTas@species_params$erepro > 0.8))]

#which species have lowest and highest erepro
print("max erepro for species")
as.character(stedTas@species_params$species[which(stedTas@species_params$erepro == max(stedTas@species_params$erepro))])
max(stedTas@species_params$erepro)

print("min erepro for species")
as.character(stedTas@species_params$species[which(stedTas@species_params$erepro == min(stedTas@species_params$erepro))])
min(stedTas@species_params$erepro)

#Put this into the original parameter file
mariaParams$erepro <- stedTas@species_params$erepro

```




### old scaling attempts 

```{r, eval = F }
median(modelBiom/initBiom)


minobs <- min(which(params@w > 1)) #which size categories in the model would be observable by RLS

#get mean model abundances in the observed weight groups (so they are comparable to surveys)
modelAb <- rep(NA, length(params@species_params$species))

## derive initial abundance scalar 
for (i in 1:length(mariaParams$species)) {
  
modelAb[i] <- (sum(params@initial_n[i,c(minobs:200)]*params@dw[c(minobs:200)]))

}

#get a scalar between observed abundances and model abundances 
#Observed mean abundances are for survey area of 500m2, so they are scaled to be for m2, to correspond to model units

observ_abund <- params@species_params$meanAb/500
abund_scalar <- observ_abund/modelAb

abund_scalar

##

totalModelAb <- (params@initial_n * (abund_scalar/100))*params@dw 
observModelAb <- (params@initial_n[,c(minobs:200)]* (abund_scalar/100)) *params@dw[c(minobs:200)]

totalModelAbF <- apply(totalModelAb, 1, sum)
obserModelAbF <- apply(observModelAb, 1, sum)

totalModelAbF
obserModelAbF
modelAb*abund_scalar
observ_abund

obserModelAbF/observ_abund

##


minw <- min(which(params@w > 0.05)) #size class of minimum weight
maxw <- max(which(params@w < mariaParams$w_inf[1]))

plot(log10(params@w[c(minw:maxw)]), (params@initial_n[1,c(minw:maxw)]*params@dw[c(minw:maxw)]*abund_scalar[1]), type  = 'l', xlab = "Log10, w,g", ylab = "log10, Abundance per m2")
for (i in 2:length(mariaParams$species)) {
points(log10(params@w), (params@initial_n[i,]*params@dw*abund_scalar[i]), type = 'l')
}


```


### 
Lobsters eat small fish and urchins. They used to be the dominant predator of small fish (hypothesis)

Southern and eastern rock lobster 
Southern calamari and snapper thermal preference 

Use relative abundance and use species specific size spectra from the first 5 years 


1. Fishing mortality using data inside and outside MPA
2. Benthic resource slope and regeneration rate and how they change with temperature and what does it mean to the community 
3. How does temperature alter the community structure inside MPA using long term data

how are these systems functioning?
how does fishing and climate affect that?