Updates to latest model parameters.

Setup_ZooMizer_grid.R

@@ -40,18 +40,27 @@ saveRDS(zoomss, file = paste0("Output/", jobname, "_ZooMSS_", ID_char,".RDS"))
 environment(project_simple) <- asNamespace('mizer')
 assignInNamespace("project_simple", project_simple, ns = "mizer")
 
-ZooMizer_coupled <- function(ID, tmax = 1000, effort = 0) {
+ZooMizer_coupled <- function(ID, tmax = 1000, effort = 0, zoomssjobname) {
   groups <- read.csv("data/TestGroups_mizer.csv") #load in ZooMSS groups
+  ID_char <- sprintf("%04d",ID)
+  #lower "starting proportion" of salps and chaetognaths
+  # groups$Prop[c(7,8)] <- 0.1*groups$Prop[c(7,8)]
+  # #lower "starting proportion" of euphausiids
+  # groups$Prop[6] <- 0.1*groups$Prop[c(6)]
+
+  # change start of senescence for salps and euphs
+  # groups$w_mat[c(6,8)] <- groups$w_mat[c(6,8)] - 2  #note these are in log10(w)
+
   zoo_groups <- groups[which(groups$Type == "Zooplankton"),] #just keep the zooplankton
 
-  input <- readRDS("data/enviro_grid20210705.RDS")[ID,]
+  input <- readRDS("data/enviro_grid20.RDS")[ID,]
 
-  stable_zoomizer <- readRDS(paste0("Output/", jobname, "_ZooMSS_", ID_char,".RDS"))
+  stable_zoomizer <- readRDS(paste0("Output/", zoomssjobname, "_ZooMSS_", ID_char,".RDS"))
 
   times <- length(getTimes(stable_zoomizer))
-  stable_zoo <- colMeans(stable_zoomizer@n[ceiling(times/2+1):times,,])
+  stable_zoo <- colMeans(stable_zoomizer@n[ceiling(times/2+1):times,,]) 
 
-  stable_zoomizer@n <- sweep(stable_zoomizer@n, 2, 2.5 * stable_zoomizer@params@species_params$Carbon, "*") #convert to carbon
+  # stable_zoomizer@n <- sweep(stable_zoomizer@n, 2, 2.5 * stable_zoomizer@params@species_params$Carbon, "*") #convert to carbon
   intercept <- mean(getCommunitySlope(stable_zoomizer,
                                       species = which(groups$Type == "Zooplankton"),
                                       max_w = 1e-5,
@@ -65,21 +74,23 @@ ZooMizer_coupled <- function(ID, tmax = 1000, effort = 0) {
   new <- stable_zoo[which(groups$Type == "Zooplankton"),]
 
   temp_eff <- 2.^((input$sst - 30)/10)
-  R_factor = 1.01
-  fish_params <- newTraitParams(no_sp = 5,
+  fish_params <- newTraitParams(no_sp = 13,
                                 min_w = 10^(-3),
-                                min_w_inf = 10^3,
+                                min_w_inf = 10^1,
                                 max_w_inf = 10^7,
-                                no_w = (7-(-3))*10+1, # sets dx = 0.1
+                                no_w = 10*(7-(-3))+1, # sets dx = 0.1
                                 min_w_pp = 10^(-14.5),
-                                alpha = 0.25, # * temp_eff, #takes care of assimilation efficiency & temp effect on Encounter
+                                alpha = 0.6, # * temp_eff, #takes care of assimilation efficiency & temp effect on Encounter
                                 kappa = exp(intercept), #10^(input$phyto_int) * 10^-1,
                                 lambda = 1 - slope, # 1 - input$phyto_slope,
-                                gamma = 1280 * temp_eff, #takes care of temp effect on PredRate and Encounter
+                                gamma = 640 * temp_eff, #takes care of temp effect on PredRate and Encounter
                                 # f0 = 0.6,
-                                # h = 10^50,
+                                h = 40,
+                                beta = 100,
+                                sigma = 1.3,
                                 # R_factor = 1.01, #RMax = RFactor * RDI, default 4. Note RDI 
-                                reproduction_level = 1/R_factor,
+                                reproduction_level = 0,
+                                fc = 0,
                                 ks = 0, #set metabolism to zero
                                 ext_mort_prop = 0, #currently zeroed since this is fitted. No need for temp effect here, calculates from PredMort
                                 knife_edge_size = 10  # min size for fishing
@@ -88,32 +99,91 @@ ZooMizer_coupled <- function(ID, tmax = 1000, effort = 0) {
   # fish_params@species_params$gamma <- fish_params@species_params$gamma * temp_eff
   # fish_params <- setSearchVolume(fish_params)
 
+  fish_params@other_params$temp_eff <- 2.^((input$sst - 30)/10)
+
   zoo_params <- newZooMizerParams(groups = zoo_groups, input = input, fish_params = fish_params)
 
   zoo_params@initial_n <- new[which(groups$Type == "Zooplankton"),]
 
+  zoo_params@species_params$interaction_resource[9] <- 1 # make jellies eat resource (i.e. fish)
+
+  ## Zooplankton Mortality
+
+  #U-shaped mortality
+  M <- getExtMort(zoo_params) # pulls out senescence
+  z0pre <- rep(0, nrow(M))
+  z0pre[c(6,8)] <- 1 * (zoo_params@species_params$w_inf[c(6,8)])^(1/4) * fish_params@other_params$temp_eff #choose groups to apply to
+  allo_mort <- outer(z0pre, w(zoo_params)^(-1/4)) # mortality of z0pre * w^(-1/4)
+  zoo_params <- setExtMort(zoo_params, ext_mort = M + allo_mort)
+
+  # add fixed level of background mortality (so in total have mizer fixed background mortality + senescence)
+  # M <- getExtMort(zoo_params) # pulls out senescence
+  # z0pre <- rep(0, nrow(M))
+  # z0pre[c(6,8)] <- 0.6 #choose groups to apply to; default Mizer value 0.6
+  # mu_b <- z0pre * zoo_params@species_params$w_inf^(-1/4) * fish_params@other_params$temp_eff
+  # zoo_params <- setExtMort(zoo_params, ext_mort = M + mu_b) # adds 0.6 * w_inf^(-1/4)
+
+  #Both
+  # M <- getExtMort(zoo_params) # pulls out senescence
+  # z0pre <- rep(0, nrow(M))
+  # z0pre[c(6,8)] <- 1 * (zoo_params@species_params$w_inf[c(6,8)])^(1/4) * fish_params@other_params$temp_eff #choose groups to apply to
+  # larv_mort <- outer(z0pre, w(zoo_params)^(-1/4)) # mortality of z0pre * w^(-1/4)
+  # 
+  # z0pre1 <- rep(0, nrow(M))
+  # z0pre1[c(6,8)] <- 6 * fish_params@other_params$temp_eff #choose groups to apply to; default Mizer value 0.6
+  # back_mort <- z0pre * zoo_params@species_params$w_inf^(-1/4)
+  # zoo_params <- setExtMort(zoo_params, ext_mort = M + larv_mort + back_mort) # adds 0.6 * w_inf^(-1/4)!
+
+  # #type II feeding for zoo
+  # zoo_params <- setRateFunction(zoo_params, "PredRate", "mizerPredRate")
+  # zoo_params <- setRateFunction(zoo_params, "FeedingLevel", "mizerFeedingLevel")
+  # zoo_params@species_params$h <- 40
+  # zoo_params <- setMaxIntakeRate(zoo_params)
+
+  #Density-dependent mortality
+  zoo_params@species_params$mu0DD <- rep(80, 9) * zoo_params@species_params$w_inf
+  zoo_params@species_params$mu0DD[6] <- zoo_params@species_params$mu0DD[6] * 50 #turn up density-dependent mortality on krill
+  zoo_params <- setRateFunction(zoo_params, "Mort", "totalMortDD")
+
+
   fish_params <- fish_params %>%
     setComponent(component = "zoo",
                  initial_value = zoo_params@initial_n,
                  dynamics_fun = "zoo_dynamics",
+                 mort_fun = "zoo_predation",
                  component_params = list(params = zoo_params)
     ) %>%
     setResource(resource_dynamics = "resource_zooMizer")
 
+
+
   # this line not needed if gamma is set to be temperature-dependent
   # fish_params <- setRateFunction(fish_params, "PredRate", "PredRate_temp")
 
-  fish_params@other_params$temp_eff <- 2.^((input$sst - 30)/10)
 
-  initialNResource(fish_params) <- resource_zooMizer(fish_params, fish_params@initial_n_other) #TODO: set this to be stable state of a regular ZooMizer run
 
-  plotSpectra(fish_params, wlim = c(10^-14.5, NA), total = FALSE)
+  # initialN(fish_params) <- get_initial_n(fish_params) * 10^3 # adjust initial n for fish
+  initialNResource(fish_params@other_params$zoo$params) <- kappa*params@w_full^(1-lambda) / params@dw_full #returns the fixed spectrum at every time step
+  initialNResource(fish_params@other_params$zoo$params)[fish_params@w_full > fish_params@other_params$zoo$params@resource_params$w_pp_cutoff* (1 - 1e-06)] <- 0
+  fish_idx <- (length(fish_params@w_full)-length(fish_params@w) + 1):length(fish_params@w_full)
+  initialNResource(fish_params@other_params$zoo$params)[fish_idx] <- initialNResource(fish_params@other_params$zoo$params)[fish_idx] + colSums(initialN(fish_params))
+  # fish_params <- setParams(fish_params)
+
+  initialNResource(fish_params) <- resource_zooMizer(params = fish_params, n_other = fish_params@initial_n_other) #TODO: set this to be stable state of a regular ZooMizer run
+  # plotSpectra(fish_params, wlim = c(10^-14.5, NA), total = FALSE)
+
 
-  initialN(fish_params) <- get_initial_n(fish_params)
-  fish_params <- setParams(fish_params)
+  # fish_params@species_params$R_max <- readRDS("data/rmaxs.RDS")/10
+  fish_params@species_params$erepro <- 1
+  # fish_params@species_params$R_max <- fish_params@resource_params$kappa * fish_params@species_params$w_inf ^(-1.5)
+  # fish_params@species_params$R_max <- fish_params@species_params$R_max * 10^(5+3*log10(input$chlo))
+  # fish_params <- setBevertonHolt(fish_params, erepro = 1)
+  # fish_params@species_params$R_max <- fish_params@species_params$R_max * 10000
+  # fish_params <- setRateFunction(fish_params, "FeedingLevel", "FeedingLevel_type3")
+  # fish_params <- setExtMort(fish_params, ext_mort = array(data=0.1, dim=dim(fish_params@mu_b)))
 
   return(project(fish_params, t_max = tmax, dt = 0.01, effort = effort))
-}  
+} 
 
 out <- ZooMizer_coupled(ID, tmax = 1000, effort = 0)
 saveRDS(out, file = paste0("Output/", jobname, "_ZooMizer_", ID_char,".RDS"))

ZooMizerChloExperiments.Rmd

@@ -343,9 +343,9 @@ gg3 <- ggplot(df4, aes(x=log10(chlo), y = Proportion, fill = Group)) +
   labs(title = "Fish proportion")
 
 library(patchwork)
-gg1
-gg2
-gg3
+gg1 / gg2 / gg3
+
+
 
 ```
 

ZooMizerResourceFunctions.R

@@ -98,6 +98,12 @@ setZooMizerConstants <- function(params, Groups, sst){
   return(params)
 }
 
+phyto_fixed <- function(params, n, n_pp, n_other, rates, dt, kappa=params@resource_params$kappa, lambda=params@resource_params$lambda, ... ) {
+  npp <- kappa*params@w_full^(1-lambda) / params@dw_full #returns the fixed spectrum at every time step
+  npp[params@w_full > params@resource_params$w_pp_cutoff* (1 - 1e-06)] <- 0
+  return(npp)
+}
+
 #' Set up MizerParams object for zooplankton model
 #' 
 #' @param groups TODO Patrick: document this
@@ -108,12 +114,13 @@ newZooMizerParams <- function(groups, input, fish_params) {
   # This is essentially the existing function fZooMizer_run() but without
   # actually running the model
 
-  kappa = 10^(input$phyto_int)
-  lambda = 1-input$phyto_slope
-  chlo = input$chlo
-  sst = input$sst
-  dt = input$dt
-  tmax = input$tmax
+  kappa <- 10^(input$phyto_int)
+  lambda <- 1-input$phyto_slope
+  w_pp_cutoff <- 10^(input$phyto_max)* (1 + 1e-06)
+  chlo <- input$chlo
+  sst <- input$sst
+  dt <- input$dt
+  tmax <- input$tmax
 
   #data
   groups$w_min <- 10^groups$w_min #convert from log10 values
@@ -128,7 +135,13 @@ newZooMizerParams <- function(groups, input, fish_params) {
 
     #todo - ramp up constant repro for coexistence
 
-  params <- new_newMultispeciesParams(species_params=groups,
+  phyto_fixed <- function(params, n, n_pp, n_other, rates, dt, kappa=resource_params(params)[["kappa"]], lambda=resource_params(params)[["lambda"]], ... ) {
+    npp <- kappa*params@w_full^(1-lambda) / params@dw_full #returns the fixed spectrum at every time step
+    npp[params@w_full > params@resource_params$w_pp_cutoff* (1 - 1e-06)] <- 0
+    return(npp)
+  }
+
+  params <- new_ZooplanktonParams(species_params=groups,
                                          interaction=NULL, #NULL sets all to 1, no strict herbivores
                                          #min_w = 10^(-10.7),
                                          #max_w = 10^7* (1 + 1e-06),
@@ -142,13 +155,22 @@ newZooMizerParams <- function(groups, input, fish_params) {
                                          resource_dynamics = "phyto_fixed",
                                          kappa = kappa,
                                          lambda = lambda,
-                                         RDD = constantRDD(species_params = groups), #first go at this
+                                         RDD = constantRDD(species_params = groups) #first go at this
                                          #pred_kernel = ... #probably easiest to just import this/pre-calculate it, once dimensions are worked out
   )
 
   params@other_params$sst <- sst
   params@other_params$chlo <- chlo
 
+
+
+
+  params@resource_dynamics<-"phyto_fixed"
+  params <- setResource(params, resource_dynamics = "phyto_fixed", kappa = kappa, lambda = lambda, w_pp_cutoff = w_pp_cutoff, n=0.7)
+
+  initialNResource(params)<-kappa*params@w_full^(1-lambda) / params@dw_full #returns the fixed spectrum at every time step
+  initialNResource(params)[params@w_full > params@resource_params$w_pp_cutoff* (1 - 1e-06)] <- 0
+
   params@species_params$w_min <- groups$w_min  #fix Mizer setting the egg weight to be one size larger for some groups.
   #params@initial_n[] <- readRDS("data/initialn.RDS")
 
@@ -289,12 +311,6 @@ resource_zooMizer <- function(params, n_other, ...) {
   return(total_n_eff)
 }
 
-phyto_fixed <- function(params, n, n_pp, n_other, rates, dt, kappa=params@resource_params$kappa, lambda=params@resource_params$lambda, ... ) {
-  npp <- kappa*params@w_full^(1-lambda) / params@dw_full #returns the fixed spectrum at every time step
-  npp[params@w_full > params@resource_params$w_pp_cutoff* (1 - 1e-06)] <- 0
-  return(npp)
-}
-
 setassim_eff <- function(groups){
   assim_eff = matrix(groups$GrossGEscale * groups$Carbon, nrow = nrow(groups), ncol = nrow(groups))
   #get_filterfeeders <- which(groups$FeedType == "FilterFeeder")
@@ -499,7 +515,7 @@ FeedingLevel_type3 <- function (params, n, n_pp, n_other, t, encounter, ...)
   return(encounter^2 / (encounter^2 + params@intake_max^2)) #sigmoidal type 3 feeding
 }
 
-new_newMultispeciesParams <- function(
+new_ZooplanktonParams <- function(
   species_params,
   interaction = NULL,
   no_w = 100,
@@ -532,7 +548,7 @@ new_newMultispeciesParams <- function(
   kappa = 1e11,
   lambda = 2.05,
   w_pp_cutoff = 10,
-  resource_dynamics = "resource_semichemostat",
+  resource_dynamics = "phyto_fixed",
   # setFishing
   gear_params = data.frame(),
   selectivity = NULL,
@@ -598,6 +614,14 @@ new_newMultispeciesParams <- function(
               selectivity = selectivity,
               catchability = catchability,
               initial_effort = initial_effort)
+  setResource(params, resource_rate = resource_rate,
+              resource_capacity = resource_capacity,
+              n = n,
+              r_pp = r_pp,
+              kappa = kappa,
+              lambda = lambda,
+              w_pp_cutoff = w_pp_cutoff,
+              resource_dynamics = resource_dynamics)
 
   params@initial_n <- get_initial_n(params)
   params@initial_n_pp <- params@cc_pp

data/TestGroups_mizer.csv

@@ -1,13 +1,13 @@
-species,Type,FeedType,Prop,w_min,w_inf,w_mat,gamma,q,PPMRscale,PPMR,FeedWidth,GrossGEscale,Carbon,Repro,PlotColour,interaction_resource
-Flagellates,Zooplankton,Heterotroph,0.1,-10.7,-6.8,-8.8,640,0.8,1.5,NA,0.36,2.5,0.15,0,seagreen3,1
-Ciliates,Zooplankton,Heterotroph,0.1,-9.3,-6.3,-8.3,640,0.8,0.04,NA,0.47,2.5,0.15,0,seagreen4,1
-Larvaceans,Zooplankton,FilterFeeder,0.1,-6.4,-3.2,-5.2,640,0.8,-3,NA,0.7,2.5,0.02,0,lightgoldenrod4,1
-OmniCopepods,Zooplankton,Omnivore,0.04,-7.5,-3.5,-5.5,640,0.8,-0.5,NA,0.57,2.5,0.12,0,cornflowerblue,1
-CarnCopepods,Zooplankton,Carnivore,0.06,-7.5,-2.5,-4.5,640,0.8,1.5,NA,0.36,2.5,0.12,0,darkred,0
-Euphausiids,Zooplankton,Omnivore,0.1,-4.19,0.2,-1.8,640,0.8,-2,NA,0.7,2.5,0.12,0,darkblue,1
-Chaetognaths,Zooplankton,Carnivore,0.1,-5.9,-0.9,-2.9,640,0.8,1,NA,0.46,2.5,0.04,0,firebrick1,0
-Salps,Zooplankton,FilterFeeder,0.01,-4.69,1.4,-0.6,640,0.8,-2.68,NA,0.7,2.5,0.02,0,lightgoldenrod3,1
-Jellyfish,Zooplankton,Carnivore,0.01,-3,2,0,640,0.8,0.73,NA,0.52,2.5,0.005,0,orange,0
-Fish_Small,Fish,Carnivore,NA,-3,2,0,640,0.8,NA,100,1.3,2.5,0.1,0,lightcoral,0
-Fish_Med,Fish,Carnivore,NA,-3,4,2,640,0.8,NA,100,1.3,2.5,0.1,0,coral3,0
-Fish_Large,Fish,Carnivore,NA,-3,7,4,640,0.8,NA,100,1.3,2.5,0.1,0,coral4,0
+"","species","Type","FeedType","Prop","w_min","w_inf","w_mat","gamma","q","PPMRscale","PPMR","FeedWidth","GrossGEscale","Carbon","Repro","PlotColour","interaction_resource"
+"1","Flagellates","Zooplankton","Heterotroph",0.1,-10.7,-6.8,-8.8,640,0.8,1.5,NA,0.36,2.5,0.15,0,"seagreen3",1
+"2","Ciliates","Zooplankton","Heterotroph",0.1,-9.3,-6.3,-8.3,640,0.8,0.04,NA,0.47,2.5,0.15,0,"seagreen4",1
+"3","Larvaceans","Zooplankton","FilterFeeder",0.1,-6.4,-3.2,-5.2,640,0.8,-3,NA,0.7,2.5,0.02,0,"lightgoldenrod4",1
+"4","OmniCopepods","Zooplankton","Omnivore",0.04,-7.5,-3.5,-5.5,640,0.8,-0.5,NA,0.57,2.5,0.12,0,"cornflowerblue",1
+"5","CarnCopepods","Zooplankton","Carnivore",0.06,-7.5,-2.5,-4.5,640,0.8,1.5,NA,0.36,2.5,0.12,0,"darkred",0
+"6","Euphausiids","Zooplankton","Omnivore",0.1,-4.19,0.2,-1.8,640,0.8,-2,NA,0.7,2.5,0.12,0,"darkblue",1
+"7","Chaetognaths","Zooplankton","Carnivore",0.1,-5.9,-0.9,-2.9,640,0.8,1,NA,0.46,2.5,0.04,0,"firebrick1",0
+"8","Salps","Zooplankton","FilterFeeder",0.01,-4.69,1.4,-0.6,640,0.8,-2.68,NA,0.7,2.5,0.02,0,"lightgoldenrod3",1
+"9","Jellyfish","Zooplankton","Carnivore",0.01,-3,2,0,640,0.8,0.73,NA,0.52,2.5,0.005,0,"orange",0
+"10","Fish_Small","Fish","Carnivore",NA,-3,2,0,640,0.8,NA,100,1.3,2.5,0.1,0,"lightcoral",0
+"11","Fish_Med","Fish","Carnivore",NA,-3,4,2,640,0.8,NA,100,1.3,2.5,0.1,0,"coral3",0
+"12","Fish_Large","Fish","Carnivore",NA,-3,7,4,640,0.8,NA,100,1.3,2.5,0.1,0,"coral4",0