TempGrowthModel_revised.Rmd

---
title: "Mechanistic temperature-size rule explanation should reconcile physiological and mortality responses to temperature "
author: "Shane A. Richards and Asta Audzijonyte"
date: "29/03/2022"
output:
  html_document:
    highlight: tango
    theme: cerulean
    toc: yes
    toc_float: no
  word_document:
    toc: yes
  pdf_document:
    toc: yes
---

# IMPORTANT!

Parameter names in this code may not fully correspond to the names used in the manuscript. 
- Length dependent mortality steepness is called zp in manuscript but it is called m1 in the code. 


```{r message=FALSE, warning=FALSE}
rm(list = ls()) # clear memory

library(tidyverse)
library(tibble)
library(ggplot2)
library(cowplot)
library(tidyr)

```

# Useful functions

```{r}
# Reserve:Structure ratio 
RSRatio <- function(age) {
  tmp <- rs1*(age-a_bar)
  tmp <- max(tmp, -20) # bound below
  tmp <- min(tmp,  20) # bound above

  return(rs_min + (rs_max - rs_min)*exp(tmp) / (1.0 + exp(tmp)))
} 

# Predator length [m]
predatorLength <- function(S) {
  return((S/l_const)^(1/3.0)) # assumes invariant growth
} 

# Mortality rate [d-1]
mortProb <- function(l) {
  m_rate <- m_min + (m_max-m_min)*exp(-m1*l) # instantaneous mortality rate
  return(1.0 - exp(-m_rate)) # probability die per day
} 

# intake rate [d d-1] 
# Egle - extended to include temperature
grossIntake <- function(S) {
  Eiscalar <- exp((-Ei/k)*(1/Temp-1/Tref))
  
  return((g0*S^g1)*Eiscalar*S^(ci*(Temp-Tref)))
}

# DEB maintenance version [g d-1 g-1]
# Egle - extended to include temperature
maintenance <- function(S, R) {
  Emscalar <- exp((-Em/k)*(1/Temp-1/Tref))
  
  # return((ms*S + mr*w*R)*Emscalar*(S+R)^(cm*(Temp-Tref)))
  return((ms*S + mr*R)*Emscalar*(S+R)^(cm*(Temp-Tref)))
}

# Energetic reproductive cost
reproCost <- function(S) {
  return(ra*S^rb)
}

# scaling factor
phi <- function(Temp, Tref, m, x) {
  if (x == "I") { # Intake
    Eiscalar <- exp((-Ei/k)*(1/Temp-1/Tref))
    Val      <- Eiscalar * m^(ci*(Temp - Tref))
  } else { # Maintenance
    Emscalar <- exp((-Em/k)*(1/Temp-1/Tref))
    Val      <- Emscalar * m^(cm*(Temp - Tref))
  }
  
  return(Val)
}
```

# Model parameters for baseline scenario

```{r message=FALSE, warning=FALSE}
# Model parameters for the baseline zero fishing scenario 
max_age_years <- 20      # years of simulation
max_age_days  <- 365*max_age_years # (d)

rs_max  <- 1.3           # maximum RS ratio
rs_min  <- 0.0           # minimum RS ratio

# length-weight conversion: 
l_const <- 1250/(0.60^3) # (g m-1) num = weight (g), denom = length (m)
# length-weight conversion uses weight of S only and assumes that 
# 1250g of S weight (ca 3000g of total) corresponds to 60cm long fish

g0      <- 0.1           # intake rate constant: intake when structural weight
                         #   is 1 g (g d-1), includes assimilation efficiency 
g1      <- 0.6667        # power to uptake rate with S weight

ms      <- 0.003         # maintenance cost of structural mass (g d-1 g-1)
mr      <- 0.0003        # maintenance cost of reversible mass (g d-1 g-1)  

s_eff   <- 0.3333        # conv. efficiency of assimilated intake to structure
r_eff   <- 0.9           # conv. efficiency of assimilated intake to reversible

# Reproductive cost function 
ra      <- 6             # reprod cost for 1 g of struct weight: (g g-1)
rb      <- 0.6           # reprod cost power w.r.t. structural weight

# Mortality parameters 
#m_min   <- 0.2/365      # background mortality rate (d-1) 

# m1                    # steepness of the length-dependent mortality rate
m_max   <- 4/365        # max length-dependent mortality rate (d-1)

s1      <- 7.0          # steepness of the condition related mortality rate
s_max   <- 4/365        # maximum condition related rate (d-1) [when R = 0]

 #Fishing mortality parameters 
# Egle - F changed to 0
Fm   <- 0               # Instantaneous fishing mortality (day-1)
Fmid <- 0.3             # length (in meters) of the 50% fishing selectivity
Fk   <- 20              # steepness of the logistic fishing function 

# Egle - temperature parameters added 

# These are important - we are looking at growth changes after 3C of warming
Temp <- 283              # ambient temperature, in Kelvin
Tref <- 280              # reference temperature, in Kelvin
k    <- 0.00008617332    # Boltzman constant 
```

# Scenarios

```{r}
scen_par <- read.csv(file="scenarioParameters.csv", header=TRUE, sep=",",  dec=".")
nscen <- length(scen_par$Sc)
nscen_used <- nscen 

# create global scenario variables extracted from scenario 
Em <- 0.0
Ei <- 0.0
cm <- 0.0
ci <- 0.0
m1 <- 0.0
m_min <- 0.0
# create global optimal LH-variables extracted from scenario 
rs1   <- 0.0
a_bar <- 0.0
w     <- 0.0

GetScenario <- function(Sc) {
  rw <- which(scen_par$Sc == Sc)

  Em <<- scen_par$Em[rw[1]]
  Ei <<- scen_par$Ei[rw[1]]
  cm <<- scen_par$cm[rw[1]]
  ci <<- scen_par$ci[rw[1]]
  m1 <<- scen_par$m1[rw[1]]
  m_min <<- scen_par$m_min[rw[1]]

  rs1   <<- scen_par$rs1[rw[1]]
  a_bar <<- scen_par$a_bar[rw[1]]
  w     <<- scen_par$w[rw[1]]

  return(length(rw) == 0) # found a valid scenario
}

#this function will need to be modified to update for the fact that scenarios are just listed in order 

# GetScenario_new <- function(Sc_new) {
#   rw <- which(scen_par$Sc == scen_used$Sc[Sc_new])
#   
#   Em <<- scen_par$Em[rw[1]]
#   Ei <<- scen_par$Ei[rw[1]]
#   cm <<- scen_par$cm[rw[1]]
#   ci <<- scen_par$ci[rw[1]]
#   m1 <<- scen_par$m1[rw[1]]
#   m_min <<- scen_par$m_min[rw[1]]
#   
#   rs1   <<- scen_par$rs1[rw[1]]
#   a_bar <<- scen_par$a_bar[rw[1]]
#   w     <<- scen_par$w[rw[1]]
# 
#   return(length(rw) == 0) # found a valid scenario
# }
```

## turn off life-history optimisation

```{r eval = TRUE}
# This code turns off selection and uses the base-line allocation strategy 
#   for all scenarios

# do_not_optimise_new <- c(4,15,17,28)
# do_not_optimise_old <- scen_used$Sc_int[do_not_optimise_new]

## if you want to run scnenarios without life-history optimisation, uncomment lines below
# do_not_optimise_old <- 2:nrow(scen_par)
# scen_par$rs1[do_not_optimise_old] <- scen_par$rs1[1]
# scen_par$a_bar[do_not_optimise_old] <- scen_par$a_bar[1]
# scen_par$w[do_not_optimise_old] <- scen_par$w[1]
```

# Fig. 3: scaling of rates

```{r}
GetScenario(17)

TC     <- 7:10
Temp_v <- TC + 273

# Egle - extended to include temperature
Eiscalar_v <- exp((-Ei/k)*(1/Temp_v-1/Tref))
Emscalar_v <- exp((-Em/k)*(1/Temp_v-1/Tref))

df_temperature <- tibble(
  Temperature = TC,
	Intake      = Eiscalar_v,
	Metabolism  = Emscalar_v
) %>%
gather(key = "Process", value = "Scalar", 2:3)

p1 <- ggplot(df_temperature) +
	geom_line(aes(x = Temperature, y = Scalar, color = Process)) +
 # scale_color_manual(values = c("#ef8a62","#67a9cf")) +
  scale_color_manual(values = c("black","grey")) +
  ylim(1,1.4) +
  labs(
    x = expression("Temperature ("*~degree*C*")"),
    y = expression("Scaling value, "*phi)) +
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    legend.position = "none")
p1
```

```{r}
mass_min <- 1    # minimum mass to plot (g)
mass_max <- 1000 # maximum mass to plot (g)
mass_n   <- 20   # masses to plot

m <- mass_min*((mass_max/mass_min)^((1:mass_n - 1)/(mass_n-1)))

#scenarios <- c(17, 19, 28) # old scenario numbers
scenarios <- c(17, 19, 41) # final scenario sequence

v_scen <- NULL
v_type <- NULL
v_mass <- NULL
v_phi  <- NULL
# v_tmet <- NULL
# tot_int <- NULL

for (i in scenarios) {
  GetScenario(i)
  
  v_scen <- c(v_scen, rep(i, mass_n))
  v_type <- c(v_type, rep("Intake", mass_n))
  v_mass <- c(v_mass, m)
  v_phi  <- c(v_phi, phi(Temp, Tref, m, "I"))
  
  v_scen <- c(v_scen, rep(i, mass_n))
  v_type <- c(v_type, rep("Metabolism", mass_n))
  v_mass <- c(v_mass, m)
  v_phi  <- c(v_phi, phi(Temp, Tref, m, "M"))
#  v_tmet <- c(v_tmet, maintenance())
}

df_phi <- tibble(
    Sc = v_scen, type = v_type, mass = v_mass, phi = v_phi
  ) %>% mutate(
    grp = paste(Sc, type, sep = "_")
  )

df_phi$scen_name <- NA
df_phi$scen_name[which(df_phi$Sc == 17)] <- "III-1"
df_phi$scen_name[which(df_phi$Sc == 19)] <- "III-3"
df_phi$scen_name[which(df_phi$Sc == 41)] <- "V-6"

```

```{r}
p2 <- ggplot(df_phi) +
 # scale_color_manual(values = c("#ef8a62","#67a9cf")) +
  scale_color_manual(values = c("black","grey")) +
  geom_line(aes(x = mass, y = phi, color = type, 
    linetype = factor(scen_name), group = grp)) +
  labs(
    x = "Mass (g)", y = expression("Scaling value, "*phi), 
    color = "Process", linetype = "Scenario") +
  ylim(0.9,1.8) +
  theme_bw() +
  theme(panel.grid = element_blank())
p2
```

```{r}
plot_grid(p1, p2, labels = c("A", "B"), rel_widths = c(1,1.5))
```

# Main loop: scenario predictions

```{r}
# useful age-dependent values
Res               <- array(data=0,c(max_age_years,365, nscen_used))
Str               <- array(data=0,c(max_age_years,365, nscen_used))
dayIntake         <- array(data=0,c(max_age_years,365, nscen_used))
dayMaintenance    <- array(data=0,c(max_age_years,365, nscen_used))
dayNetIntake      <- array(data=0,c(max_age_years,365, nscen_used))
dayLambda         <- array(data=0,c(max_age_years,365, nscen_used))
dayPredatorLength <- array(data=0,c(max_age_years,365, nscen_used))
daySurvival       <- array(data=0,c(max_age_years,365, nscen_used))
yearSpawn         <- array(data=0,c(max_age_years, nscen_used))
yearRepCost       <- array(data=0,c(max_age_years, nscen_used))
yearFitness       <- array(data=0,c(max_age_years, nscen_used))
natmortality      <- array(data=0,c(max_age_years,365, nscen_used))
fishmortality     <- array(data=0,c(max_age_years,365, nscen_used))
dayGrowth         <- array(data=0,c(max_age_years,365, nscen_used))
relGrowth         <- array(data=0,c(max_age_years,365, nscen_used))
strGrowth         <- array(data=0,c(max_age_years,365, nscen_used))

for (aa in 1: nscen_used) { # aa is the scenario number
  
  #a1 <- which(scen_par$Sc == scen_used$Sc[aa])[1] #old code when scenario numbers were changing
  # extract scenario patrameters and optimal LH
 a1 <- aa
  Em <- scen_par$Em[a1]
  Ei <- scen_par$Ei[a1] 
  cm <- scen_par$cm[a1]
  ci <- scen_par$ci[a1] 
  m1 <- scen_par$m1[a1]
  m_min <- scen_par$m_min[a1]
  
  Eiscalar <- exp((-Ei/k)*(1/Temp-1/Tref))
  Emscalar <- exp((-Em/k)*(1/Temp-1/Tref))
  
  rs1   <- scen_par$rs1[a1]  
  a_bar <- scen_par$a_bar[a1]
  w     <- scen_par$w[a1]
  
  # set initial weight
  S0 <- 1             # initial weight 
  R0 <- S0*RSRatio(0) # enforce correct initial reserve-structural ratio
  S0 <- 1/(1+RSRatio(0))
  R0 <- RSRatio(0) *S0

	# set day 1 year 1 mass
  Res[1,1,aa] <- R0
  Str[1,1,aa] <- S0

  daySurvival[1,1,aa] <- 1.0 # initially all individuals are alive  

  # perform the simulation
  for (yr in 1:max_age_years) { # for each year
     for (day in 1:364) { # for each day
      Rstart <- Res[yr,day,aa] # starting mass (reserve)
      Sstart <- Str[yr,day,aa] # starting mass (structure)
    
      # calc length using str
      dayPredatorLength[yr,day,aa] <- (Sstart/l_const)^(1/3.0) 
      fish_mort <- Fm / (1+exp(-Fk*(dayPredatorLength[yr,day,aa]-Fmid)))
      mort_rate <- m_min + 
      	(m_max-m_min)*exp(-m1*dayPredatorLength[yr,day,aa]) + # predation
        s_max*exp(-s1*Rstart/Sstart)   +                      # starve
        fish_mort                                             # fishing
      
      fishmortality[yr,day+1,aa] <- fish_mort   
      mort_prob <- 1.0 - exp(-mort_rate) # probability die this day
      # prob alive
      daySurvival[yr,day+1,aa] <- (1-mort_prob)*daySurvival[yr,day,aa] 
      natmortality[yr,day+1,aa] <- m_min + 
      	(m_max-m_min)*exp(-m1*dayPredatorLength[yr,day,aa]) + # predation
        s_max*exp(-s1*Rstart/Sstart)                          # starve

      # temperature-dependent intake
      intake <- g0*Sstart^g1 * 
      	Eiscalar * Sstart^(ci*(Temp-Tref))
    
      respiration <- (ms*Sstart + mr*Rstart) * 
      	Emscalar * (Sstart + Rstart)^(cm*(Temp-Tref)) # (g d-1)
      net_intake  <- intake - respiration             # (g d-1)

      dayIntake[yr, day,aa]      <- intake
      dayMaintenance[yr, day,aa] <- respiration
      dayNetIntake[yr, day,aa]   <- net_intake
      age                        <- 365*(yr-1) + day # age of animal (days)
      
      # add bounds to prevent numerical issues when calculating lambda
      tmp <- rs1*(age-a_bar)
      tmp <- max(tmp, -20) # bound below
      tmp <- min(tmp,  20) # bound above
      dayLambda[yr, day,aa] <- rs_min + (rs_max - rs_min)*exp(tmp) /
        (1.0 + exp(tmp)) # RS ratio = strategy
    
      if (net_intake >= 0) {
        dR <- r_eff*net_intake # maximum R allocation
        dS <- s_eff*net_intake # maximum S allocation
        # use Lambdamax instead of dayLambda and set w=1  
        if (dayLambda[yr, day,aa]*Sstart > Rstart) { # need to bump up reserves
          r_take <- min(dayLambda[yr, day,aa]*Sstart - Rstart, dR) 
          Rstart <- Rstart + r_take
          net_intake <- net_intake - r_take/r_eff
        } else { # need to bump up structure
          s_take <- min(Rstart/dayLambda[yr, day,aa] - Sstart, dS) 
          Sstart <- Sstart + s_take
          net_intake <- net_intake - s_take/s_eff
        }
        
        # partition remaining mass to keep desired ratio
        Res[yr,day+1,aa] <- Rstart + dayLambda[yr, day,aa] * r_eff * s_eff *
        	net_intake / (r_eff + dayLambda[yr, day,aa]*s_eff) 
        Str[yr,day+1,aa] <- Sstart + r_eff*s_eff*net_intake / 
        	(r_eff +     dayLambda[yr, day,aa]*s_eff)       
      } else {
      	# what should be taken from R given the conversion inefficiencies 
        dR <- (1/r_eff)*net_intake 
        newR <- (Rstart + dR)
        if (newR < 0) {
          newR <- 0
        }
     
        Res[yr,day+1,aa] <- newR
        Str[yr,day+1,aa] <- Sstart
      }  

      # percentage change in weight over the day
      Growth <- (Res[yr,day+1,aa]+Str[yr,day+1,aa]) -
      	(Res[yr,day,aa]+Str[yr,day,aa]) 
      percGrowth           <- (Growth/(Res[yr,day,aa]+Str[yr,day,aa])) * 100
      sGrowth              <- Growth/Str[yr,day,aa] * 100
      dayGrowth[yr,day,aa] <-  Growth
      relGrowth[yr,day,aa] <- percGrowth
      strGrowth[yr,day,aa] <- sGrowth
    }

    dayPredatorLength[yr,365,aa] <- dayPredatorLength[yr,364,aa]
    dayGrowth[yr,365,aa]         <- dayGrowth[yr,364,aa]
    relGrowth[yr,365,aa]         <- relGrowth[yr,364,aa]
    strGrowth[yr,365,aa]         <- strGrowth[yr,364,aa]

    # perform spawning
    repro_cost <- ra*Str[yr,365,aa]^rb # fixed cost of reproduction
    # spawning mass after cost
    spawn_mass <- max(0, w*Res[yr,365,aa] - repro_cost) 
  
    if (spawn_mass > 0) { # enough to spawn?
      Res[yr,365,aa]     <- Res[yr,365,aa] - spawn_mass - repro_cost
      yearSpawn[yr,aa]   <- spawn_mass
      yearRepCost[yr,aa] <- repro_cost
      yearFitness[yr,aa] <- spawn_mass*daySurvival[yr,365,aa]
    } 
  
    if (yr < max_age_years) {
      Str[yr+1,1,aa]         <- Str[yr,365,aa]
      Res[yr+1,1,aa]         <- Res[yr,365,aa]
      daySurvival[yr+1,1,aa] <- daySurvival[yr,365,aa]
      natmortality[yr+1,1,aa]  <- natmortality[yr,365,aa]
      fishmortality[yr+1,1,aa]  <- fishmortality[yr,365,aa]
    }
  }
}
```

# Fitness

```{r}
## Fitnesses of all scenarios
v_LifetimeFitness <- apply(yearFitness, 2, sum) # expected life-time fitness
df_LifetimeFitness <- tibble(
	Sc = scen_par$Sc,
	LifetimeFitness = v_LifetimeFitness
) %>%
left_join(scen_par, by = "Sc")

ggplot(df_LifetimeFitness) +
	geom_point(aes(x = sequence, y = LifetimeFitness)) +
	theme_bw() +
	labs(y = "Expected lifetime spawning biomass", x = "Scenario sequence (Table S1)") +
	theme()
```

# Fig. 4 - absolute rates

```{r}
#TSR consistent scenarios (sequence)
#sc= c(17, 19, 41)

days = c(2:364)
sc = 1

str.v <- NULL
mass.v <- NULL
maint.v <- NULL
int.v <- NULL

mass.t <- Res[,,sc] + Str[,,sc]

#go through 20 years and turn arrays into vectors for plotting
for (i in 1:20) {
str.v <- c(str.v, Str[i,days,sc])
len.v <- predatorLength(str.v)
mass.v <- c(mass.v, mass.t[i,days])
maint.v <- c(maint.v, dayMaintenance[i,days,sc])
int.v <- c(int.v,dayIntake[i,days,sc])
}

#plot(mass.v, maint.v, type = 'l', ylim = c(0, 30), xlim = c(0, 6000))
# plot(str.v, maint.v, type = 'l', ylim = c(0, 30), xlim = c(0, 4500))
# points(str.v, int.v, type = 'l', lty = 2)

plot(len.v, maint.v, type = 'l', ylim = c(0, 25), xlab = "Fish length (m)", ylab = "Daily intake and maintenance (g)")
points(len.v, int.v, type = 'l', lty = 2)

sc = 17

str.v <- NULL
mass.v <- NULL
maint.v <- NULL
int.v <- NULL

mass.t <- Res[,,sc] + Str[,,sc]

for (i in 1:20) {
str.v <- c(str.v, Str[i,days,sc])
len.v <- predatorLength(str.v)
mass.v <- c(mass.v, mass.t[i,days])
maint.v <- c(maint.v, dayMaintenance[i,days,sc])
int.v <- c(int.v,dayIntake[i,days,sc])
}

# points(str.v, maint.v, type = 'l', col = 'lightgray', lwd = 2)
# points(str.v, int.v, type = 'l', lty = 2, col = 'lightgray', lwd = 2)
points(len.v, maint.v, type = 'l', col = 'lightgray', lwd = 2)
points(len.v, int.v, type = 'l', lty = 2, col = 'lightgray', lwd = 2)


sc = 19

str.v <- NULL
mass.v <- NULL
maint.v <- NULL
int.v <- NULL

mass.t <- Res[,,sc] + Str[,,sc]

for (i in 1:20) {
str.v <- c(str.v, Str[i,days,sc])
mass.v <- c(mass.v, mass.t[i,days])
maint.v <- c(maint.v, dayMaintenance[i,days,sc])
int.v <- c(int.v,dayIntake[i,days,sc])
}

# points(str.v, maint.v, type = 'l', col = 'darkgray', lwd = 2)
# points(str.v, int.v, type = 'l', lty = 2, col = 'darkgray', lwd = 2)
points(len.v, maint.v, type = 'l', col = 'darkgray', lwd = 2)
points(len.v, int.v, type = 'l', lty = 2, col = 'darkgray', lwd = 2)

legend("topleft", 
  legend = c("base maint", "base intake", "III-1 maint", "III-1 intake", "III-3 maint", "III-3 intake"), 
  #col = c(rgb(0.2,0.4,0.1,0.7), rgb(0.8,0.4,0.1,0.7)), 
 # pch = c(17,19), 
 lty = c(1, 2, 1, 2, 1, 2), 
 lwd = c(1, 1, 2, 2, 2, 2), 
 col = c("black","black","lightgray","lightgray","darkgray","darkgray"),
  bty = "n", 
  pt.cex = 1, 
  cex = 1, 
  text.col = "black", 
  horiz = F , 
  inset = c(0.05, 0.05))


```

# make scenario dataframe

```{r fig.width=9, fig.height=7}
param_value <- NULL
param_type  <- NULL
age         <- NULL
scenario    <- NULL
baseline    <- NULL

#if plotting for first 8 years uncomment this
#max_age_years <- 8

# calculate year of maturation
maturation_age <- rep(1, nscen_used)
for (i in 1:nscen_used) {
  for (j in 1:max_age_years) {
   	if (yearSpawn[j,i] == 0) {
   	  maturation_age[i] <- maturation_age[i] + 1
   	}
  }
}

# create some useful summary statistics
scen_par <- scen_par %>%
	mutate(
		Em_gt_Ei = Em > Ei,
		cm_gt_ci = cm > ci
	)

df_scenario <- tibble(
  scenario = scen_par$sequence,
	maturation_age = maturation_age
)
df_scenario$Em_gt_Ei <- scen_par$Em_gt_Ei#[scen_used$Sc_int]
df_scenario$cm_gt_ci <- scen_par$cm_gt_ci#[scen_used$Sc_int]
df_scenario$m1 <- scen_par$m1#[scen_used$Sc_int]
df_scenario$m1 <- factor(df_scenario$m1)

# create a large age-length relation for each scenario
for (aa in 1: nscen_used) { # for each scenario
  for (yr in 1:max_age_years) { # for each year
    param_value <- c(param_value, dayPredatorLength[yr, ,aa]) # save reserves
    age         <- c(age, 1:365 + 365*(yr-1))            # add next year
    param_type  <- c(param_type, rep("Reversible", 365)) # add parameter
    #scenario    <- c(scenario, rep(scen_used$Sc_int[aa],365)) 
    scenario    <- c(scenario, rep(scen_par$sequence[aa],365)) # add scenario
  }
}

# place values into a data frame
df_length <- tibble(
	age = age, 
	type = param_type, 
	value = param_value, 
	scenario = scenario
) %>%
left_join(df_scenario, by = "scenario")

df_length$maturation_age <- factor(df_length$maturation_age)

# create a reference data set and scenario data set
df_baseline <- filter(df_length, scenario == "1") %>%
	dplyr::select(-scenario)
df_other    <- filter(df_length, scenario != "1")
names(df_other)[4] <- "sequence"
df_other <- left_join(df_other, scen_par, by = "sequence")

rm(df_length) # don't need this now
```

### ### 
# Plots

# Figure A1: allocation strategy

code and plot suggested by the reviewer Jan Kozlowski (thank you)

```{r}
#Eq. 1
#parameters
lambda.min<-0
lambda.max<-1.3

r<-c(0.001,0.002)
a.dash<-c(200,500)

lambda<-function(a,r,a.dash)
{temp<-lambda.max-lambda.min
temp1<-exp(r*(a-a.dash))
temp1<-temp1/(1+temp1)
lambda<-lambda.min+temp*temp1}

a<-seq(0,3650,1)
lambda.a<-lambda(a,r[1],a.dash[1])
a.year<-a/365
plot(a.year,lambda.a, type="l",xlab="Age (years)",ylab="Target reversible:structural ratio")
lambda.a<-lambda(a,r[2],a.dash[1])
lines(a.year,lambda.a, type="l",lty='dashed')
lambda.a<-lambda(a,r[1],a.dash[2])
lines(a.year,lambda.a, type="l",col="grey", lwd = 2)
lambda.a<-lambda(a,r[2],a.dash[2])
lines(a.year,lambda.a, type="l",col="grey",lty='dashed', lwd = 2)
legend(6.5,0.8,legend=c("r=0.001, a.dash=200", "r=0.002, a.dash=200","r=0.001, a.dash=500", "r=0.002, a.dash=500"),
       col=c("black", "black","grey","grey"), title=paste("Lambda.max=1.3; Lambda.min=0"),lty=c(1,2,1,2), cex=0.8)

```

### old version, not used now

```{r fig.height = 3, fig.width=5}
# bounds for plotting
l_min   <- 0.0           # minimum length for plotting (m)
l_max   <- 0.45          # maximum length for plotting (m)
w_min   <- 1             # minimum mortality rate for plotting (d-1)
w_max   <- 2000          # maximum mortality rate for plotting (d-1)
  
# prepare data frame to display RS ratio
res <- GetScenario(1)

vec_age <- seq(from = 0, to = max_age_days, by = 1)
df_RS   <- tibble(
	Age = vec_age, 
	RSratio = sapply(vec_age, FUN = RSRatio)
)

p1 <- ggplot(df_RS, aes(x = Age, y = RSratio)) + 
	geom_line() + 
  ylim(0,rs_max) + 
	labs(
	  x = "Age (days)", 
	  y = "Ratio (Reserve:Structure)", 
    subtitle = "Strategy: desired R:S ratio"
	) +
  theme_bw()
p1

```

```{r}
# prepare data frame to display reproductive costs
vec_weight <- seq(from = w_min, to = w_max, length.out = 101)
df_weight <- tibble(
	Weight = vec_weight, 
	Cost = reproCost(vec_weight)
)

p2 <- ggplot(df_weight, aes(x = Weight, y = Cost)) + 
	geom_line() + 
  geom_abline(intercept = 0, slope = 1, color = "grey") +
	labs(
	  x = "Structural weight (g)", 
	  y = "Reproductive cost (g)", 
    subtitle = "Reproduction cost (non-spawn mass)"
	) +
    xlim (0, 500) +
	theme_bw()
p2
```

# Fig. A2: Mortality 

code and plot suggested by the reviewer Jan Kozlowski (thank you)

```{r}

#Eq. 3a
zp<-c(4,6,8,10)
#m.min recalculated for a day basis; originally 0.2 and 0.4
m.min<-c(0.0005479,0.00109589)
m.max<-0.0109589
#length, in meters, from 0 to 1, by 1 cm
L<-seq(0,1, by = 0.01)

m.l<-function(zp,m.min,m.max,L)
{m.l<-m.min+(m.max-m.min)*exp(-zp*L)}

m.l1<-m.l(zp[1],m.min[1],m.max,L)
m.l2<-m.l(zp[1],m.min[2],m.max,L)
m.l3<-m.l(zp[2],m.min[1],m.max,L)
m.l4<-m.l(zp[2],m.min[2],m.max,L)
m.l5<-m.l(zp[3],m.min[1],m.max,L)
m.l6<-m.l(zp[3],m.min[2],m.max,L)
m.l7<-m.l(zp[4],m.min[1],m.max,L)
m.l8<-m.l(zp[4],m.min[2],m.max,L)

#c1<-expression("M['P,min']")

df.l<-data.frame(L,m.l1,m.l2,m.l3,m.l4)
fig2<-ggplot()+
  theme_classic()+
  theme(legend.position = c(0.7,0.8))+
  xlab("Fish length (m)")+
  ylab("Instantaneous mortality rate (per day)")+
  geom_blank(aes(show.legend=TRUE))+
  
  geom_line(data=df.l, aes(x=L,y=m.l1,color="zp = 4"))+
  geom_line(data=df.l, aes(x=L,y=m.l2,color="zp = 4"),linetype="dashed")+
  
  geom_line(data=df.l, aes(x=L,y=m.l3,color="zp = 6"))+
  geom_line(data=df.l, aes(x=L,y=m.l4,color="zp = 6"),linetype="dashed")+

  geom_line(data=df.l, aes(x=L,y=m.l5,color="zp = 8"))+
  geom_line(data=df.l, aes(x=L,y=m.l6,color="zp = 8"),linetype="dashed")+
  
  geom_line(data=df.l, aes(x=L,y=m.l7,color="zp = 10"))+
  geom_line(data=df.l, aes(x=L,y=m.l8,color="zp = 10"),linetype="dashed")+
  
  scale_color_manual(name=bquote("solid "~M[P_min]~"= 0.2; dashed"~M[P_min]~"= 0.4"),values=c("zp = 4"="black",
                          "zp = 6"="blue","zp = 8"="orange","zp = 10"="red"))
  
fig2


```

### old version, not used now

```{r}
# Predator length [m]
predatorLength <- function(S) {
  return((S/l_const)^(1/3.0)) # assumes invariant growth
} 

# Mortality rate [d-1]
mortProb <- function(l) {
  m_rate <- m_min + (m_max-m_min)*exp(-m1*l) # instantaneous mortality rate
  return(1.0 - exp(-m_rate)) # probability die per day
} 

m_max   <- 4/365 
# m_min <- 0.2/365 # 0.4
# m1 <- 8 # 6, 4, 10

S <- (1:1000)
l_const = 1250/(0.60^3)
length <- predatorLength(S)

m_min <- 0.2/365 # 0.4
m1 <- 8 # 6, 4, 10
baseMort <- mortProb(length)

#only m_min increases - scen 15
m_min <- 0.4/365 # 0.4
m1 <- 8 # 6, 4, 10
scen16 <- mortProb(length)

#only m1 (or Zp in the manusript) decreases
m_min <- 0.2/365 # 0.4
m1 <- 4 # 6, 4, 10
scen14 <- mortProb(length)

#m1 (or Zp in the manusript) decreases, but not quite as much
m_min <- 0.2/365 # 0.4
m1 <- 6 # 6, 4, 10
scen25 <- mortProb(length)

#m_min increases and m1 decreases 
m_min = 0.4/365
m1 = 4
scen15 <- mortProb(length)

#m_min increases, but m1 also increases, so that mortality declines faster 
m_min = 0.4/365
m1 = 10
scen33 <- mortProb(length)

plot(length, baseMort, type ="l", lwd = 2, xlab = "Length, m", ylim = c(0, 0.01), ylab = "Instantaneous mortality rate (per day)")
points(length, scen16, type = 'l', lty = 2, lwd = 2)
points(length, scen14, type = 'l', lty = 3, lwd = 2)
points(length, scen25, type = 'l', lwd = 2, lty = 5, col = 'gray')
points(length, scen15, type = 'l', lty = 6, lwd = 2)
points(length, scen33, type = 'l', col  = 'gray', lwd =2)

legend("topright", 
  legend = c("m_min = 0.2, zp = 8", "m_min = 0.4, zp = 8", "m_min = 0.2, zp = 4", "m_min = 0.2, zp = 6", "m_min = 0.4, zp = 4", "m_min = 0.4, zp = 10"), 
  #col = c(rgb(0.2,0.4,0.1,0.7), rgb(0.8,0.4,0.1,0.7)), 
 # pch = c(17,19), 
 lty = c(1, 2, 3, 5, 6, 1), 
 lwd = c(2, 2, 2, 2, 2, 2), 
 col = c("black","black","black","gray","black","gray"),
  bty = "n", 
  pt.cex = 2, 
  cex = 1.2, 
  text.col = "black", 
  horiz = F , 
  inset = c(0.1, 0.1))


```

### Growth (Length)

# make growth plot dataframe

```{r}
param_value <- NULL
param_type  <- NULL
age         <- NULL
scenario    <- NULL
baseline    <- NULL

# create a large age-length relation for each scenario
for (aa in 1: nscen_used) { # for each scenario
  for (yr in 1:max_age_years) { # for each year
    param_value <- c(param_value, Str[yr, ,aa] + Res[yr, , aa]) # save reserves
    age         <- c(age, 1:365 + 365*(yr-1))            # add next year
    param_type  <- c(param_type, rep("Mass", 365)) # add parameter
    scenario    <- c(scenario, rep(scen_par$sequence[aa],365))              # add scenario
  }
}

# place values into a data frame
df_mass <- tibble(
	age = age, 
	type = param_type, 
	value = param_value, 
	scenario = scenario
) %>%
left_join(df_scenario, by = "scenario")

df_mass$maturation_age <- factor(df_mass$maturation_age)

# create a reference data set and scenario data set
df_baseline_mass <- filter(df_mass, scenario == "1") %>%
	dplyr::select(-scenario)
df_other_mass    <- filter(df_mass, scenario != "1")
names(df_other_mass)[4] <- "sequence"
df_other_mass <- left_join(df_other_mass, scen_par, by = "sequence")
rm(df_mass) # don't need this now
```

# growth in all scenarios

## by length

```{r fig.width=9, fig.height=6}
# colour by age at maturation
ggplot() + 
#  scale_color_manual(values = c("#ef8a62","#67a9cf")) +
  scale_color_manual(values = c("gray","black")) +
  geom_line(data = df_other, 
  	aes(x = age, y = value, color = maturation_age), size = 1) +
  geom_line(data = df_baseline, 
  	aes(x = age, y = value), color = "black", size = 0.5, linetype = 2) + 
  labs(
  	x ="Age (d)",
	  y = "Length (m)",
  	color = "Age at\nmaturation (y)",
  	subtitle = "Dashed line is baseline, solid line is scenario"
  ) + 
	# scale_y_log10() +
  facet_wrap( ~ sequence, ncol = 5) +   
  theme_bw() + 
  theme(
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(),
    axis.text=element_text(size=8)
  )
```

## by weight

```{r fig.width=9, fig.height=6}
# colour by age at maturation
ggplot() + 
#  scale_color_manual(values = c("#ef8a62","#67a9cf")) +
  scale_color_manual(values = c("gray","black")) +
  geom_line(data = df_other_mass, 
  	aes(x = age, y = value, color = maturation_age), size = 1) +
  geom_line(data = df_baseline_mass, 
  	aes(x = age, y = value), color = "black", size = 0.5, linetype = 2) + 
  labs(
  	x ="Age (d)",
	  y = "Weight (g)",
  	color = "Age at\nmaturation (y)",
  	subtitle = "Black is baseline, coloured is scenario"
  ) + 
	scale_y_log10() +
  facet_wrap( ~ sequence, ncol = 6) +   
  theme_bw() + 
  theme(
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank()
  )
```

#Figure 2: four scenarios weight 

```{r}
# colour by age at maturation
p1 <- ggplot() + 
#  scale_color_manual(values = c("#ef8a62","#67a9cf")) +
  scale_color_manual(values = c("gray","black")) +
  geom_line(data = filter(df_other_mass, sequence %in% c(12,17,19,41)),
  	aes(x = age, y = value, color = maturation_age), size = 1) +
  geom_line(data = df_baseline_mass, 
  	aes(x = age, y = value), color = "black", size = 0.5, linetype = 2) + 
  labs(
  	x ="Age (d)",
	  y = "Weight (g)",
  	color = "Age at\nmaturation (y)"
  ) + 
	scale_y_log10() +
    xlim(0,3000) +
  facet_wrap( ~ sequence) +   
  theme_bw() + 
  theme(
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(),
    legend.position = "none"
  )
p1
# colour by age at maturation
p2 <- ggplot() + 
#  scale_color_manual(values = c("#ef8a62","#67a9cf")) +
  scale_color_manual(values = c("gray","black")) +
  geom_line(data = filter(df_other_mass, sequence %in% c(12,17,19,41), age < 3*365),
  	aes(x = age, y = value, color = maturation_age), size = 1) +
  geom_line(data = filter(df_baseline_mass, age < 3*365),
  	aes(x = age, y = value), color = "black", size = 0.5) + 
  labs(
  	x ="Age (d)",
	  y = "Weight (g)",
  	color = "Age at\nmaturation (y)"
  ) + 
  facet_wrap( ~ sequence) +   
  theme_bw() + 
  theme(
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(),
    legend.position = "none"
  )

plot_grid(p1, p2, labels = c("A", "B"), ncol = 1)
```

### ### 

## old stuff - do not run

## old code to sort out scenario numbers

```{r}
#different scenarios with already optimised parameters (in excel)  
scen_par <- read.csv(file="Scenarios3.csv", header=TRUE, sep=",",  dec=".")

nscen <- length(scen_par$Sc)

scen_used <- read.csv(file="FinalScenarioNos.csv", header=FALSE, fileEncoding="UTF-8-BOM", sep=",",  dec=".")
names(scen_used) <- "Sc_int"
# number of scenarios
scen_used <- scen_used %>%
  mutate(Sc = paste("Sc", as.character(Sc_int), sep = ""))
nscen_used <- length(scen_used$Sc)
scen_used$Sc_new <- 1:nrow(scen_used)

# ## save scenario parameter values 
scen_par$scen_new <- scen_used$Sc_new[match(scen_par$Sc, scen_used$Sc)]
scenarioParams <- scen_par[which(is.na(scen_par$scen_new) == F),-c(1,11)]
scenarioParams$scen_numb <- scenarioParams$scen_new
scenarioParams <- scenarioParams %>%
  mutate(scen_new = paste("Sc", as.character(scen_new), sep = ""))
scenarioParams <- scenarioParams %>% arrange(scen_numb)
 write.csv(scenarioParams, file = "scenarioParameters.csv")

scen_used1 <- read.csv(file="RevisedScenOrder.csv", header=TRUE, sep=",",  dec=".")
scenarioParams$sequence <- scen_used1$sequence[match(scenarioParams$scen_numb, scen_used1$scen_numb)]
scenarioParams$finalScName <- scen_used1$ï..name[match(scenarioParams$scen_numb, scen_used1$scen_numb)]
scenarioParams <- scenarioParams %>% arrange(sequence)
```