diff --git a/Setup_ZooMizer_grid.R b/Setup_ZooMizer_grid.R index 5e9fd6f..4d86dea 100644 --- a/Setup_ZooMizer_grid.R +++ b/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")) diff --git a/ZooMizerChloExperiments.Rmd b/ZooMizerChloExperiments.Rmd index 4d0f873..fdb63d0 100644 --- a/ZooMizerChloExperiments.Rmd +++ b/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 + + ``` diff --git a/ZooMizerResourceFunctions.R b/ZooMizerResourceFunctions.R index 86d7855..013bb97 100644 --- a/ZooMizerResourceFunctions.R +++ b/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 diff --git a/data/TestGroups_mizer.csv b/data/TestGroups_mizer.csv index efce046..c42ebd4 100644 --- a/data/TestGroups_mizer.csv +++ b/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 diff --git a/data/enviro_grid_sst15.RDS b/data/enviro_grid_sst15.RDS new file mode 100644 index 0000000..039c0cd Binary files /dev/null and b/data/enviro_grid_sst15.RDS differ diff --git a/data/enviro_sst_15_chlox25.RDS b/data/enviro_sst_15_chlox25.RDS new file mode 100644 index 0000000..fb67ebf Binary files /dev/null and b/data/enviro_sst_15_chlox25.RDS differ diff --git a/data/hpcmeanerepros.RDS b/data/hpcmeanerepros.RDS new file mode 100644 index 0000000..5b44235 Binary files /dev/null and b/data/hpcmeanerepros.RDS differ diff --git a/data/rmaxs.RDS b/data/rmaxs.RDS new file mode 100644 index 0000000..76c3370 Binary files /dev/null and b/data/rmaxs.RDS differ diff --git a/runzoomizerparallel.R b/runzoomizerparallel.R index 5b7f566..ac983cd 100644 --- a/runzoomizerparallel.R +++ b/runzoomizerparallel.R @@ -1,122 +1,302 @@ +# Script to run ZooMizer coupled model over a range of environmental conditions +# Patrick Sykes +# 1 July 2021 +# Updated 10 March 2023 + require(mizer) require(tidyverse) require(foreach) require(assertthat) -ZooMizer_coupled <- function(ID, tmax = 25, effort = 0) { - groups <- read.csv("data/TestGroups_mizer.csv") #load in ZooMSS groups - zoo_groups <- groups[which(groups$Type == "Zooplankton"),] #just keep the zooplankton - input <- readRDS("data/enviro_test20.RDS")[ID,] - - stable_zoomizer <- readRDS("test_grid_20210317.RDS")[[ID]] - - times <- length(getTimes(stable_zoomizer)) - 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 - intercept <- mean(getCommunitySlope(stable_zoomizer, - species = which(groups$Type == "Zooplankton"), - max_w = 1e-5, - )$intercept[ceiling(times/2+1):times]) - slope <- mean(getCommunitySlope(stable_zoomizer, - species = which(groups$Type == "Zooplankton"), - max_w = 1e-5, - )$slope[ceiling(times/2+1):times]) - - - new <- stable_zoo[which(groups$Type == "Zooplankton"),] - - temp_eff <- 2.^((input$sst - 30)/10) - - fish_params <- newTraitParams(no_sp = 5, - min_w = 10^(-3), - min_w_inf = 10^3, - max_w_inf = 10^7, - no_w = (7-(-3))*10+1, # sets dx = 0.1 basically - min_w_pp = 10^(-14.5), - alpha = 0.25, # * 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 - # f0 = 0.6, - # h = 10^50, - R_factor = 1.01, #RMax = RFactor * RDI, default 4. Note RDI - 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 - ) - - # fish_params@species_params$gamma <- fish_params@species_params$gamma * temp_eff - # fish_params <- setSearchVolume(fish_params) - - zoo_params <- newZooMizerParams(groups = zoo_groups, input = input, fish_params = fish_params) - - zoo_params@initial_n <- new[which(groups$Type == "Zooplankton"),] - - fish_params <- fish_params %>% - setComponent(component = "zoo", - initial_value = zoo_params@initial_n, - dynamics_fun = "zoo_dynamics", - 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") + +library(doParallel) +cl <- makePSOCKcluster(max(1, detectCores()-1)) ## set up cores-1 machines +registerDoParallel(cl, cores = (max(1, detectCores()-1))) +clusterEvalQ(cl, lapply(c("mizer", "assertthat", "tidyverse"), require, character.only = TRUE)) %>% invisible() +clusterExport(cl, c("ZooMizer_coupled")) + +jobname <- "20230418_erepro_1_Rmax_Julia_sst15" + +sims <- foreach(ID = 1:21, + .packages = c("mizer", "assertthat", "tidyverse") + # .export = c("ZooMizer_coupled") +) %dopar% { + start <- Sys.time() + zoomssjobname <- "20230404_sst15" + ID_char <- sprintf("%04d",ID) + source("uncoupledmodel.R") + source("ZooMizerResourceFunctions.R") + + # environment(new_project_simple) <- asNamespace('mizer') + # assignInNamespace("project_simple", new_project_simple, ns = "mizer") + + # Run the ZooMSS model inside Mizer: + # phyto_fixed <- function(params, n, n_pp, n_other, rates, dt, kappa = 10^enviro$phyto_int[i], lambda= 1-enviro$phyto_slope[i], ... ) { + # npp <- kappa*params@w_full^(-lambda) #returns the fixed spectrum at every time step + # npp[params@w_full > params@resource_params$w_pp_cutoff] <- 0 + # return(npp) + # } + + # density-dependent mortality + # totalMortDD <- function(params, n, n_pp, n_other, t, f_mort, pred_mort, ...){ + # mort <- pred_mort + params@mu_b + f_mort + # # Add contributions from other components + # # for (i in seq_along(params@other_mort)) { + # # mort <- mort + + # # do.call(params@other_mort[[i]], + # # list(params = params, + # # n = n, n_pp = n_pp, n_other = n_other, t = t, + # # component = names(params@other_mort)[[i]], ...)) + # # } + # ddmort_raw <- sweep(n, "w", params@w^params@species_params$q[1] * params@dw, "*") + # ddmort <- ddmort_raw * 0 + # sp_sel <- 1:nrow(ddmort) #c(6,7,8,9) + # mu0 <- 80 # Benoit & Rochet value = 80, but no density-independent mortality in that model + # ddmort[sp_sel,] <- mu0 * ddmort_raw[sp_sel,] #* params@species_params$w_inf + # ddmort[6,] <- ddmort[6,] * 50 + # return(mort + ddmort) + # } + + # zoomss <- fZooMizer_run(groups = read.csv("data/TestGroups_mizer.csv"), input = readRDS("data/enviro_sst_15_chlox25.RDS")[ID,], no_w = (1*177+1)) # dx=0.1 + # saveRDS(zoomss, file = paste0("Output/", zoomssjobname, "_ZooMSS_", ID_char,".RDS")) + + mid <- Sys.time() - fish_params@other_params$temp_eff <- 2.^((input$sst - 30)/10) + 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_sst_15_chlox25.RDS")[ID,] + + stable_zoomizer <- readRDS(paste0("Output/", zoomssjobname, "_ZooMSS_", ID_char,".RDS")) + + + times <- length(getTimes(stable_zoomizer)) + 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 + )$intercept[ceiling(times/2+1):times]) + slope <- mean(getCommunitySlope(stable_zoomizer, + species = which(groups$Type == "Zooplankton"), + max_w = 1e-5 + )$slope[ceiling(times/2+1):times]) + new <- colMeans(stable_zoomizer@n[ceiling(times/2+1):times,,])[which(groups$Type == "Zooplankton"),] + + + ## Alternative method using existing ZooMizer run + # stable_zoomizer <- readRDS("zooplankton_sst15.RDS")[[ID]] + + # stable_zoo <- colMeans(stable_zoomizer@params@initial_n_other$zoo) #@n[ceiling(times/2+1):times,,]) + # sel <- which(stable_zoomizer@params@other_params$zoo$params@w <= 1e-5) + # fit <- lm(log(stable_zoomizer@params@other_params$zoo$params@w[sel]) ~ log(stable_zoo[sel])) + # intercept <- fit$coefficients[1] + # slope <- fit$coefficients[2] + + # new <- stable_zoomizer@params@initial_n_other$zoo[which(groups$Type == "Zooplankton"),] + + temp_eff <- 2.^((input$sst - 30)/10) + fish_params <- newTraitParams(no_sp = 13, + min_w = 10^(-3), + min_w_inf = 10^1, + max_w_inf = 10^7, + no_w = 10*(7-(-3))+1, # sets dx = 0.1 + min_w_pp = 10^(-14.5), + 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 = 640 * temp_eff, #takes care of temp effect on PredRate and Encounter + # f0 = 0.6, + h = 40, + beta = 100, + sigma = 1.3, + # R_factor = 1.01, #RMax = RFactor * RDI, default 4. Note RDI + 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 + ) + + # 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") + + + + # initialN(fish_params) <- get_initial_n(fish_params) * 10^3 # adjust initial n for fish + initialNResource(fish_params@other_params$zoo$params) <- exp(intercept)*fish_params@w_full^(slope) / fish_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) + + + # 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) #relationship from Julia's model + # fish_params@species_params$erepro <- readRDS("data/hpcmeanerepros.RDS") # mean of big HPC run - benefit of same across all runs + # fish_params@species_params$R_max <- fish_params@species_params$R_max * 10^(5+3*log10(input$chlo)) #scale R_max with chlorophyll + # fish_params <- setBevertonHolt(fish_params, erepro = 1) + # fish_params@species_params$R_max <- fish_params@species_params$R_max * 10000get + # 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)) + } - initialNResource(fish_params) <- resource_zooMizer(fish_params, fish_params@initial_n_other) #TODO: set this to be stable state of a regular ZooMizer run + environment(project_simple) <- asNamespace('mizer') + assignInNamespace("project_simple", project_simple, ns = "mizer") - plotSpectra(fish_params, wlim = c(10^-14.5, NA), total = FALSE) + sim <- ZooMizer_coupled(ID, tmax = 1000, effort = 0, zoomssjobname) + end <- Sys.time() + sim@params@other_params$time <- data.frame(zoomss = mid-start, zoomizer = end-mid) + sim +} - initialN(fish_params) <- get_initial_n(fish_params) # TODO: adjust n0_mult and a parameters - fish_params <- setParams(fish_params) - -return(project(fish_params, t_max = tmax, dt = 0.01, effort = effort)) -} - - -# PredRate_temp <- function (params, n, n_pp, n_other, t, feeding_level, ...) -# { -# no_sp <- dim(params@interaction)[1] -# no_w <- length(params@w) -# no_w_full <- length(params@w_full) -# if (length(params@ft_pred_kernel_p) == 1) { -# n_total_in_size_bins <- sweep(n, 2, params@dw, "*", -# check.margin = FALSE) -# pred_rate <- sweep(params@pred_kernel, c(1, 2), (1 - -# feeding_level) * params@search_vol * n_total_in_size_bins, -# "*", check.margin = FALSE) -# pred_rate <- colSums(aperm(pred_rate, c(2, 1, 3)), dims = 1) -# return(pred_rate * params@other_params$temp_eff) -# } -# idx_sp <- (no_w_full - no_w + 1):no_w_full -# Q <- matrix(0, nrow = no_sp, ncol = no_w_full) -# Q[, idx_sp] <- sweep((1 - feeding_level) * params@search_vol * -# n, 2, params@dw, "*") -# pred_rate <- Re(t(mvfft(t(params@ft_pred_kernel_p) * mvfft(t(Q)), -# inverse = TRUE)))/no_w_full -# pred_rate[pred_rate < 1e-18] <- 0 -# return(pred_rate * params@ft_mask * params@other_params$temp_eff) +saveRDS(sims, paste0("Output/",jobname,"_ZooMizer.RDS")) + +enviro <- readRDS("data/enviro_sst_15_chlox25.RDS") +zoo_params <- readRDS("Output/20220118_alldt_pt01_ZooMSS_0001.RDS")@params %>% removeSpecies(c("Fish_Small", "Fish_Med", "Fish_Large")) + +# sims <- foreach(ID = 1:20, +# .packages = c("mizer", "assertthat", "tidyverse") +# #.export = c("ZooMizer_coupled", "PredRate_temp") +# ) %dopar% { +# ID_char <- sprintf("%04d",ID) +# readRDS(paste0("Output/20220212_type1approx_ZooMizer_", ID_char,".RDS")) # } +# check biomass time series plots +timeseries <- foreach(ID = 1:21) %dopar% { + source("ZooMizerPlots.R") + source("ZooMizerSummaryFunctions.R") + plotBiomass_ZooMizer(sims[[ID]], zoo_params)+ + labs(title = paste0("SST = ", enviro$sst[ID], ", chlo = ", round(enviro$chlo[ID],3))) +} -library(doParallel) -cl <- makePSOCKcluster(max(1, detectCores()-1)) ## set up cores-1 machines -registerDoParallel(cl, cores = (max(1, detectCores()-1))) -clusterEvalQ(cl, lapply(c("mizer", "assertthat", "tidyverse"), require, character.only = TRUE)) %>% invisible() -clusterExport(cl, c("ZooMizer_coupled", "PredRate_temp")) +library(patchwork) +timeseriesgrid <- wrap_plots(timeseries, nrow = 7, ncol = 3) + plot_layout(guides = "collect") +timeseriesgrid +ggsave(paste0(jobname,"_timeseriesplots.png"), timeseriesgrid, width = 15, height = 28) +# biomasses <- foreach(i=1:length(sims), .combine = rbind) %dopar% sum(colMeans(tail(getBiomass(sims[[i]]),25))) +#averaged size spectra plots +spectra <- foreach(i = 1:21) %dopar% { + plotSpectra_ZooMizer(sims[[i]], zoo_params, time_range = c(501,1000), wlim = c(1e-14, NA))+ + labs(title = paste0("SST = ", enviro$sst[i], ", chlo = ", round(enviro$chlo[i],3))) +} +spectplots <- wrap_plots(spectra, nrow = 7, ncol = 3) + plot_layout(guides = "collect") +ggsave(paste0(jobname,"spectraplots.png"), spectplots, width = 15, height = 28) -sims <- foreach(ID = 1:20, - .packages = c("mizer", "assertthat", "tidyverse") - #.export = c("ZooMizer_coupled", "PredRate_temp") -) %dopar% { - source("ZooMizerResourceFunctions.R") - ZooMizer_coupled(ID, tmax = 100, effort = 0) + +#get biomass data +b <- foreach(ID=1:21, .combine = rbind) %dopar% { + bioms <- getBiomass_ZooMizer(sims[[ID]], sims[[ID]]@params@other_params$zoo$params) + rows <- ceiling(nrow(bioms)/2):nrow(bioms) + colMeans(bioms[rows,]) } +saveRDS(b, paste0(jobname,"_biomdf.RDS")) + +df <- cbind(enviro[1:21,], b) +propdf <- df +propdf[,c(3:5,11:12)] <- propdf[,c(3:5,11:12)] / rowSums(propdf[,c(3:5,11:12)]) +propdf[,13:19] <- propdf[,13:19] / rowSums(propdf[,13:19]) +propdf[,20:32] <- propdf[,20:32] / rowSums(propdf[,20:32]) + +# plot group biomass vs chlo +# Microzoo/phyto +df2 <- gather(propdf, Group, Proportion, c(pico_biom:micro_biom,Flagellates:Ciliates)) +gg1 <- ggplot(df2, aes(x=log10(chlo), y = Proportion, fill = Group)) + + geom_area() + + labs(title = "Microzoo proportion") + +#zoo +df3 <- gather(propdf, Group, Proportion, Larvaceans:Jellyfish) +gg2 <- ggplot(df3, aes(x=log10(chlo), y = Proportion, fill = Group)) + + geom_area() + + labs(title = "Zooplankton proportion") + +# fish +df4 <- gather(propdf, Group, Proportion, "1":"13") +gg3 <- ggplot(df4, aes(x=log10(chlo), y = Proportion, fill = Group)) + + geom_area() + + labs(title = "Fish proportion") -saveRDS(sims, "zoomizeroutput.RDS") \ No newline at end of file +#combine and save plots +biomchlo <- wrap_plots(gg1,gg2,gg3, nrow = 3, ncol = 1) + plot_layout(guides = "collect") +ggsave(paste0(jobname, "_biomvschlo.png"), biomchlo, width = 18, height = 15) diff --git a/tileplots.R b/tileplots.R index a8cae30..ecfec45 100644 --- a/tileplots.R +++ b/tileplots.R @@ -36,7 +36,8 @@ biom.df <- cbind(enviro, all_bioms) %>% mutate(AllFish = fish1 + fish2 + fish3 + fish4 +fish5) -tiles <- foreach(i = 9:23) %dopar% { +today <- Sys.Date() +tiles <- foreach(i = 9:ncol(biom.df)) %dopar% { column <- sym(colnames(biom.df)[i]) ggplot(biom.df, aes(x=sst, y = log10(chlo)))+ geom_tile(aes(fill=log10(!!column)))+ diff --git a/uncoupledmodel.R b/uncoupledmodel.R index da536d5..d1c41e2 100644 --- a/uncoupledmodel.R +++ b/uncoupledmodel.R @@ -332,10 +332,10 @@ fZooMizer_run <- function(groups, input, no_w = 178){ 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 ) - + mf.params@species_params$w_min <- groups$w_min #fix Mizer setting the egg weight to be one size larger for some groups. #mf.params@initial_n[] <- readRDS("data/initialn.RDS") @@ -349,7 +349,8 @@ fZooMizer_run <- function(groups, input, no_w = 178){ mf.params <- setZooMizerConstants(params = mf.params, Groups = groups, sst = input$sst) #mf.params@initial_n[] <- readRDS("data/initialn.RDS") - + mf.params@initial_n_pp <- kappa*mf.params@w_full^(1-lambda) / mf.params@dw_full #returns the fixed spectrum at every time step + mf.params@initial_n_pp[mf.params@w_full > mf.params@resource_params$w_pp_cutoff* (1 - 1e-06)] <- 0 #mf.params <- setParams(mf.params) # mf.params <- setRateFunction(mf.params, "PredRate", "new_PredRate") @@ -378,351 +379,356 @@ fZooMizer_run <- function(groups, input, no_w = 178){ ### try to fix issues with weight vectors disagreeing with ZooMSS - allow for w_full to be specified. # This involves slight editing to newMultispeciesParams and un-commenting lines in emptyParams functions -# new_newMultispeciesParams <- function( -# species_params, -# interaction = NULL, -# no_w = 100, -# #w_full = NA, #added this line -# min_w = 0.001, -# max_w = NA, -# min_w_pp = NA, -# # setPredKernel() -# pred_kernel = NULL, -# # setSearchVolume() -# search_vol = NULL, -# # setMaxIntakeRate() -# intake_max = NULL, -# # setMetabolicRate() -# metab = NULL, -# p = 0.7, -# # setExtMort -# z0 = NULL, -# z0pre = 0.6, -# z0exp = n - 1, -# # setReproduction -# maturity = NULL, -# repro_prop = NULL, -# RDD = "BevertonHoltRDD", -# # setResource -# resource_rate = NULL, -# resource_capacity = NULL, -# n = 2 / 3, -# r_pp = 10, -# kappa = 1e11, -# lambda = 2.05, -# w_pp_cutoff = 10, -# resource_dynamics = "resource_semichemostat", -# # setFishing -# gear_params = data.frame(), -# selectivity = NULL, -# catchability = NULL, -# initial_effort = NULL) { -# no_sp <- nrow(species_params) -# -# ## For backwards compatibility, allow r_max instead of R_max -# if (!("R_max" %in% names(species_params)) && -# "r_max" %in% names(species_params)) { -# names(species_params)[names(species_params) == "r_max"] <- "R_max" -# } -# -# ## Create MizerParams object ---- -# params <- new_emptyParams(species_params, -# gear_params, -# # w_full = w_full, -# no_w = no_w, -# min_w = min_w, -# max_w = max_w, -# min_w_pp = min_w_pp) -# -# ## Fill the slots ---- -# params <- params %>% -# set_species_param_default("n", n) %>% -# set_species_param_default("p", p) -# params <- set_species_param_default(params, "q", -# lambda - 2 + params@species_params$n) -# if (is.null(interaction)) { -# interaction <- matrix(1, nrow = no_sp, ncol = no_sp) -# } -# params <- -# setParams(params, -# # setInteraction -# interaction = interaction, -# # setPredKernel() -# pred_kernel = pred_kernel, -# # setSearchVolume() -# search_vol = search_vol, -# # setMaxIntakeRate() -# intake_max = intake_max, -# # setMetabolicRate() -# metab = metab, -# # setExtMort -# z0 = z0, -# z0pre = z0pre, -# z0exp = z0exp, -# # setReproduction -# maturity = maturity, -# repro_prop = repro_prop, -# RDD = RDD, -# # setResource -# resource_rate = resource_rate, -# resource_capacity = resource_capacity, -# r_pp = r_pp, -# kappa = kappa, -# lambda = lambda, -# n = n, -# w_pp_cutoff = w_pp_cutoff, -# resource_dynamics = resource_dynamics, -# # setFishing -# gear_params = gear_params, -# selectivity = selectivity, -# catchability = catchability, -# initial_effort = initial_effort) -# -# params@initial_n <- get_initial_n(params) -# params@initial_n_pp <- params@cc_pp -# params@A <- rep(1, nrow(species_params)) -# -# return(params) -# } -# -# new_emptyParams <- function(species_params, -# gear_params = data.frame(), -# no_w = 100, -# min_w = 0.001, -# w_full = NA, -# max_w = NA, -# min_w_pp = 1e-12) { -# assert_that(is.data.frame(species_params), -# is.data.frame(gear_params), -# no_w > 10) -# -# ## Set defaults ---- -# if (is.na(min_w_pp)) min_w_pp <- 1e-12 -# species_params <- set_species_param_default(species_params, "w_min", min_w) -# min_w <- min(species_params$w_min) -# -# species_params <- validSpeciesParams(species_params) -# gear_params <- validGearParams(gear_params, species_params) -# -# if (is.na(max_w)) { -# max_w <- max(species_params$w_inf) -# } else { -# if (max(species_params$w_inf) > max_w * (1 + 1e-6)) { # The fudge factor -# # is there to avoid false alerts due to rounding errors. -# too_large <- species_params$species[max_w < species_params$w_inf] -# stop("Some of your species have an maximum size larger than max_w: ", -# toString(too_large)) -# } -# } -# -# # Set up grids ---- -# if (missing(w_full)) { -# # set up logarithmic grids -# dx <- log10(max_w / min_w) / (no_w - 1) -# # Community grid -# w <- 10^(seq(from = log10(min_w), by = dx, length.out = no_w)) -# # dw[i] = w[i+1] - w[i]. Following formula works also for last entry dw[no_w] -# dw <- (10^dx - 1) * w -# # To avoid issues due to numerical imprecision -# min_w <- w[1] -# -# # For fft methods we need a constant log bin size throughout. -# # Therefore we use as many steps as are necessary so that the first size -# # class includes min_w_pp. -# x_pp <- rev(seq(from = log10(min_w), -# to = log10(min_w_pp), -# by = -dx)) - dx -# w_full <- c(10^x_pp, w) -# # If min_w_pp happened to lie exactly on a grid point, we now added -# # one grid point too much which we need to remove again -# if (w_full[2] == min_w_pp) { -# w_full <- w_full[2:length(w_full)] -# } -# no_w_full <- length(w_full) -# dw_full <- (10^dx - 1) * w_full -# } else { -# # # use supplied w_full -# no_w_full <- length(w_full) - 1 -# dw_full <- diff(w_full) -# w_full <- w_full[seq_along(dw_full)] -# # # Check that sizes are increasing -# if (any(dw_full <= 0)) { -# stop("w_full must be increasing.") -# } -# w_min_idx <- match(min_w, w_full) -# if (is.na(w_min_idx)) { -# stop("w_min must be contained in w_full.") -# } -# w <- w_full[w_min_idx:no_w_full] -# dw <- dw_full[w_min_idx:no_w_full] -# no_w <- length(w) -# min_w_pp <- w_full[1] -# } -# -# # Basic arrays for templates ---- -# no_sp <- nrow(species_params) -# species_names <- as.character(species_params$species) -# gear_names <- unique(gear_params$gear) -# mat1 <- array(0, dim = c(no_sp, no_w), -# dimnames = list(sp = species_names, w = signif(w,3))) -# ft_pred_kernel <- array(NA, dim = c(no_sp, no_w_full), -# dimnames = list(sp = species_names, k = 1:no_w_full)) -# ft_mask <- plyr::aaply(species_params$w_inf, 1, -# function(x) w_full < x, .drop = FALSE) -# -# selectivity <- array(0, dim = c(length(gear_names), no_sp, no_w), -# dimnames = list(gear = gear_names, sp = species_names, -# w = signif(w, 3))) -# catchability <- array(0, dim = c(length(gear_names), no_sp), -# dimnames = list(gear = gear_names, sp = species_names)) -# initial_effort <- rep(0, length(gear_names)) -# names(initial_effort) <- gear_names -# -# interaction <- array(1, dim = c(no_sp, no_sp), -# dimnames = list(predator = species_names, -# prey = species_names)) -# -# vec1 <- as.numeric(rep(NA, no_w_full)) -# names(vec1) <- signif(w_full, 3) -# -# # Round down w_min to lie on grid points and store the indices of these -# # grid points in w_min_idx -# w_min_idx <- as.vector(suppressWarnings( -# tapply(species_params$w_min, 1:no_sp, -# function(w_min, wx) max(which(wx <= w_min)), wx = w))) -# # function(w_min, wx) max(which(round(log10(wx), digits = 2) <= round(log10(w_min), digits = 2))), wx = w))) -# # Due to rounding errors this might happen: -# w_min_idx[w_min_idx == -Inf] <- 1 -# names(w_min_idx) <- species_names -# species_params$w_min <- w[w_min_idx] -# -# # if (any(species_params$w_min < w[w_min_idx]) || any(species_params$w_min > w[w_min_idx + 1])) { -# # msg <- "The `w_min_idx` should point to the start of the size bin containing the egg size `w_min`." -# # errors <- c(errors, msg) -# # } -# -# # Colour and linetype scales ---- -# # for use in plots -# # Colour-blind-friendly palettes -# # From http://dr-k-lo.blogspot.co.uk/2013/07/a-color-blind-friendly-palette-for-r.html -# # cbbPalette <- c("#000000", "#009E73", "#e79f00", "#9ad0f3", "#0072B2", "#D55E00", -# # "#CC79A7", "#F0E442") -# # From http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-colorblind-friendly-palette -# # cbbPalette <- c("#E69F00", "#56B4E9", "#009E73", -# # "#F0E442", "#0072B2", "#D55E00", "#CC79A7") -# # Random palette gemerated pm https://medialab.github.io/iwanthue/ -# colour_palette <- c("#815f00", -# "#6237e2", -# "#8da600", -# "#de53ff", -# "#0e4300", -# "#430079", -# "#6caa72", -# "#ee0053", -# "#007957", -# "#b42979", -# "#142300", -# "#a08dfb", -# "#644500", -# "#04004c", -# "#b79955", -# "#0060a8", -# "#dc8852", -# "#007ca9", -# "#ab003c", -# "#9796d9", -# "#472c00", -# "#b492b0", -# "#140000", -# "#dc8488", -# "#005c67", -# "#5c585a") -# # type_palette <- c("solid", "dashed", "dotdash", "longdash", -# # "twodash") -# type_palette <- c("solid") -# -# if ("linecolour" %in% names(species_params)) { -# linecolour <- species_params$linecolour -# # If any NA's first fill them with unused colours -# linecolour[is.na(linecolour)] <- -# setdiff(colour_palette, linecolour)[1:sum(is.na(linecolour))] -# # if there are still NAs, start from beginning of palette again -# linecolour[is.na(linecolour)] <- -# colour_palette[1:sum(is.na(linecolour))] -# } else { -# linecolour <- rep(colour_palette, length.out = no_sp) -# } -# names(linecolour) <- as.character(species_names) -# linecolour <- c(linecolour, "Total" = "black", "Resource" = "green", -# "Background" = "grey", "Fishing" = "red") -# -# if ("linetype" %in% names(species_params)) { -# linetype <- species_params$linetype -# linetype[is.na(linetype)] <- "solid" -# } else { -# linetype <- rep(type_palette, length.out = no_sp) -# } -# names(linetype) <- as.character(species_names) -# linetype <- c(linetype, "Total" = "solid", "Resource" = "solid", -# "Background" = "solid", "Fishing" = "solid") -# -# # Make object ---- -# # Should Z0, rrPP and ccPP have names (species names etc)? -# params <- new( -# "MizerParams", -# w = w, -# dw = dw, -# w_full = w_full, -# dw_full = dw_full, -# w_min_idx = w_min_idx, -# maturity = mat1, -# psi = mat1, -# initial_n = mat1, -# intake_max = mat1, -# search_vol = mat1, -# metab = mat1, -# mu_b = mat1, -# ft_pred_kernel_e = ft_pred_kernel, -# ft_pred_kernel_p = ft_pred_kernel, -# pred_kernel = array(), -# gear_params = gear_params, -# selectivity = selectivity, -# catchability = catchability, -# initial_effort = initial_effort, -# rr_pp = vec1, -# cc_pp = vec1, -# sc = w, -# initial_n_pp = vec1, -# species_params = species_params, -# interaction = interaction, -# other_dynamics = list(), -# other_encounter = list(), -# other_mort = list(), -# rates_funcs = list( -# Rates = "mizerRates", -# Encounter = "mizerEncounter", -# FeedingLevel = "mizerFeedingLevel", -# EReproAndGrowth = "mizerEReproAndGrowth", -# PredRate = "mizerPredRate", -# PredMort = "mizerPredMort", -# FMort = "mizerFMort", -# Mort = "mizerMort", -# ERepro = "mizerERepro", -# EGrowth = "mizerEGrowth", -# ResourceMort = "mizerResourceMort", -# RDI = "mizerRDI", -# RDD = "BevertonHoltRDD"), -# resource_dynamics = "resource_semichemostat", -# other_params = list(), -# initial_n_other = list(), -# A = as.numeric(rep(NA, no_sp)), -# linecolour = linecolour, -# linetype = linetype, -# ft_mask = ft_mask -# ) -# -# return(params) -# } +new_newMultispeciesParams <- function( + species_params, + interaction = NULL, + no_w = 100, + #w_full = NA, #added this line + min_w = 0.001, + max_w = NA, + min_w_pp = NA, + # setPredKernel() + pred_kernel = NULL, + # setSearchVolume() + search_vol = NULL, + # setMaxIntakeRate() + intake_max = NULL, + # setMetabolicRate() + metab = NULL, + p = 0.7, + # setExtMort + z0 = NULL, + z0pre = 0.6, + z0exp = n - 1, + # setReproduction + maturity = NULL, + repro_prop = NULL, + RDD = "BevertonHoltRDD", + # setResource + resource_rate = NULL, + resource_capacity = NULL, + n = 2 / 3, + r_pp = 10, + kappa = 1e11, + lambda = 2.05, + w_pp_cutoff = 10, + resource_dynamics = "resource_semichemostat", + # setFishing + gear_params = data.frame(), + selectivity = NULL, + catchability = NULL, + initial_effort = NULL) { + no_sp <- nrow(species_params) + + ## For backwards compatibility, allow r_max instead of R_max + if (!("R_max" %in% names(species_params)) && + "r_max" %in% names(species_params)) { + names(species_params)[names(species_params) == "r_max"] <- "R_max" + } + + ## Create MizerParams object ---- + params <- new_emptyParams(species_params, + gear_params, + # w_full = w_full, + no_w = no_w, + min_w = min_w, + max_w = max_w, + min_w_pp = min_w_pp) + + ## Fill the slots ---- + params <- params %>% + set_species_param_default("n", n) %>% + set_species_param_default("p", p) + params <- set_species_param_default(params, "q", + lambda - 2 + params@species_params$n) + if (is.null(interaction)) { + interaction <- matrix(1, nrow = no_sp, ncol = no_sp) + } + params <- + setParams(params, + # setInteraction + interaction = interaction, + # setPredKernel() + pred_kernel = pred_kernel, + # setSearchVolume() + search_vol = search_vol, + # setMaxIntakeRate() + intake_max = intake_max, + # setMetabolicRate() + metab = metab, + # setExtMort + z0 = z0, + z0pre = z0pre, + z0exp = z0exp, + # setReproduction + maturity = maturity, + repro_prop = repro_prop, + RDD = RDD, + # setResource + resource_rate = resource_rate, + resource_capacity = resource_capacity, + r_pp = r_pp, + kappa = kappa, + lambda = lambda, + n = n, + w_pp_cutoff = w_pp_cutoff, + resource_dynamics = resource_dynamics, + # setFishing + gear_params = gear_params, + selectivity = selectivity, + catchability = catchability, + initial_effort = initial_effort) + + params@initial_n <- get_initial_n(params) + params@initial_n_pp <- params@cc_pp + params@A <- rep(1, nrow(species_params)) + + return(params) +} + +new_emptyParams <- function(species_params, + gear_params = data.frame(), + no_w = 100, + min_w = 0.001, + w_full = NA, + max_w = NA, + min_w_pp = 1e-12) { + assert_that(is.data.frame(species_params), + is.data.frame(gear_params), + no_w > 10) + + ## Set defaults ---- + if (is.na(min_w_pp)) min_w_pp <- 1e-12 + species_params <- set_species_param_default(species_params, "w_min", min_w) + min_w <- min(species_params$w_min) + + species_params <- validSpeciesParams(species_params) + gear_params <- validGearParams(gear_params, species_params) + + if (is.na(max_w)) { + max_w <- max(species_params$w_inf) + } else { + if (max(species_params$w_inf) > max_w * (1 + 1e-6)) { # The fudge factor + # is there to avoid false alerts due to rounding errors. + too_large <- species_params$species[max_w < species_params$w_inf] + stop("Some of your species have an maximum size larger than max_w: ", + toString(too_large)) + } + } + + # Set up grids ---- + if (missing(w_full)) { + # set up logarithmic grids + dx <- log10(max_w / min_w) / (no_w - 1) + # Community grid + w <- 10^(seq(from = log10(min_w), by = dx, length.out = no_w)) + # dw[i] = w[i+1] - w[i]. Following formula works also for last entry dw[no_w] + dw <- (10^dx - 1) * w + # To avoid issues due to numerical imprecision + min_w <- w[1] + + # For fft methods we need a constant log bin size throughout. + # Therefore we use as many steps as are necessary so that the first size + # class includes min_w_pp. + x_pp <- rev(seq(from = log10(min_w), + to = log10(min_w_pp), + by = -dx)) - dx + w_full <- c(10^x_pp, w) + # If min_w_pp happened to lie exactly on a grid point, we now added + # one grid point too much which we need to remove again + if (w_full[2] == min_w_pp) { + w_full <- w_full[2:length(w_full)] + } + no_w_full <- length(w_full) + dw_full <- (10^dx - 1) * w_full + } else { + # # use supplied w_full + no_w_full <- length(w_full) - 1 + dw_full <- diff(w_full) + w_full <- w_full[seq_along(dw_full)] + # # Check that sizes are increasing + if (any(dw_full <= 0)) { + stop("w_full must be increasing.") + } + w_min_idx <- match(min_w, w_full) + if (is.na(w_min_idx)) { + stop("w_min must be contained in w_full.") + } + w <- w_full[w_min_idx:no_w_full] + dw <- dw_full[w_min_idx:no_w_full] + no_w <- length(w) + min_w_pp <- w_full[1] + } + + # Basic arrays for templates ---- + no_sp <- nrow(species_params) + species_names <- as.character(species_params$species) + gear_names <- unique(gear_params$gear) + mat1 <- array(0, dim = c(no_sp, no_w), + dimnames = list(sp = species_names, w = signif(w,3))) + ft_pred_kernel <- array(NA, dim = c(no_sp, no_w_full), + dimnames = list(sp = species_names, k = 1:no_w_full)) + ft_mask <- plyr::aaply(species_params$w_inf, 1, + function(x) w_full < x, .drop = FALSE) + + selectivity <- array(0, dim = c(length(gear_names), no_sp, no_w), + dimnames = list(gear = gear_names, sp = species_names, + w = signif(w, 3))) + catchability <- array(0, dim = c(length(gear_names), no_sp), + dimnames = list(gear = gear_names, sp = species_names)) + initial_effort <- rep(0, length(gear_names)) + names(initial_effort) <- gear_names + + interaction <- array(1, dim = c(no_sp, no_sp), + dimnames = list(predator = species_names, + prey = species_names)) + + vec1 <- as.numeric(rep(NA, no_w_full)) + names(vec1) <- signif(w_full, 3) + + # Round down w_min to lie on grid points and store the indices of these + # grid points in w_min_idx + w_min_idx <- as.vector(suppressWarnings( + tapply(species_params$w_min, 1:no_sp, + function(w_min, wx) max(which(wx <= w_min)), wx = w))) + # function(w_min, wx) max(which(round(log10(wx), digits = 2) <= round(log10(w_min), digits = 2))), wx = w))) + # Due to rounding errors this might happen: + w_min_idx[w_min_idx == -Inf] <- 1 + names(w_min_idx) <- species_names + species_params$w_min <- w[w_min_idx] + + # if (any(species_params$w_min < w[w_min_idx]) || any(species_params$w_min > w[w_min_idx + 1])) { + # msg <- "The `w_min_idx` should point to the start of the size bin containing the egg size `w_min`." + # errors <- c(errors, msg) + # } + + # Colour and linetype scales ---- + # for use in plots + # Colour-blind-friendly palettes + # From http://dr-k-lo.blogspot.co.uk/2013/07/a-color-blind-friendly-palette-for-r.html + # cbbPalette <- c("#000000", "#009E73", "#e79f00", "#9ad0f3", "#0072B2", "#D55E00", + # "#CC79A7", "#F0E442") + # From http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-colorblind-friendly-palette + # cbbPalette <- c("#E69F00", "#56B4E9", "#009E73", + # "#F0E442", "#0072B2", "#D55E00", "#CC79A7") + # Random palette gemerated pm https://medialab.github.io/iwanthue/ + colour_palette <- c("#815f00", + "#6237e2", + "#8da600", + "#de53ff", + "#0e4300", + "#430079", + "#6caa72", + "#ee0053", + "#007957", + "#b42979", + "#142300", + "#a08dfb", + "#644500", + "#04004c", + "#b79955", + "#0060a8", + "#dc8852", + "#007ca9", + "#ab003c", + "#9796d9", + "#472c00", + "#b492b0", + "#140000", + "#dc8488", + "#005c67", + "#5c585a") + # type_palette <- c("solid", "dashed", "dotdash", "longdash", + # "twodash") + type_palette <- c("solid") + + if ("linecolour" %in% names(species_params)) { + linecolour <- species_params$linecolour + # If any NA's first fill them with unused colours + linecolour[is.na(linecolour)] <- + setdiff(colour_palette, linecolour)[1:sum(is.na(linecolour))] + # if there are still NAs, start from beginning of palette again + linecolour[is.na(linecolour)] <- + colour_palette[1:sum(is.na(linecolour))] + } else { + linecolour <- rep(colour_palette, length.out = no_sp) + } + names(linecolour) <- as.character(species_names) + linecolour <- c(linecolour, "Total" = "black", "Resource" = "green", + "Background" = "grey", "Fishing" = "red") + + if ("linetype" %in% names(species_params)) { + linetype <- species_params$linetype + linetype[is.na(linetype)] <- "solid" + } else { + linetype <- rep(type_palette, length.out = no_sp) + } + names(linetype) <- as.character(species_names) + linetype <- c(linetype, "Total" = "solid", "Resource" = "solid", + "Background" = "solid", "Fishing" = "solid") + + # Make object ---- + # Should Z0, rrPP and ccPP have names (species names etc)? + params <- new( + "MizerParams", + metadata = list(), + mizer_version = packageVersion("mizer"), + extensions = vector(mode = "character"), + time_created = lubridate::now(), + time_modified = lubridate::now(), + w = w, + dw = dw, + w_full = w_full, + dw_full = dw_full, + w_min_idx = w_min_idx, + maturity = mat1, + psi = mat1, + initial_n = mat1, + intake_max = mat1, + search_vol = mat1, + metab = mat1, + mu_b = mat1, + ft_pred_kernel_e = ft_pred_kernel, + ft_pred_kernel_p = ft_pred_kernel, + pred_kernel = array(), + gear_params = gear_params, + selectivity = selectivity, + catchability = catchability, + initial_effort = initial_effort, + rr_pp = vec1, + cc_pp = vec1, + sc = w, + initial_n_pp = vec1, + species_params = species_params, + interaction = interaction, + other_dynamics = list(), + other_encounter = list(), + other_mort = list(), + rates_funcs = list( + Rates = "mizerRates", + Encounter = "mizerEncounter", + FeedingLevel = "mizerFeedingLevel", + EReproAndGrowth = "mizerEReproAndGrowth", + PredRate = "mizerPredRate", + PredMort = "mizerPredMort", + FMort = "mizerFMort", + Mort = "mizerMort", + ERepro = "mizerERepro", + EGrowth = "mizerEGrowth", + ResourceMort = "mizerResourceMort", + RDI = "mizerRDI", + RDD = "BevertonHoltRDD"), + resource_dynamics = "resource_semichemostat", + other_params = list(), + initial_n_other = list(), + A = as.numeric(rep(NA, no_sp)), + linecolour = linecolour, + linetype = linetype, + ft_mask = ft_mask + ) + + return(params) +} # setmort_test <- function(params, sst){ # #### CALCULATES CONSTANT BITS OF THE MODEL FUNCTIONS FOR EACH GROUP diff --git a/uncoupledmodel_gridrun.R b/uncoupledmodel_gridrun.R index c888008..753d44a 100644 --- a/uncoupledmodel_gridrun.R +++ b/uncoupledmodel_gridrun.R @@ -1,5 +1,5 @@ source("fZooMSS_CalculatePhytoParam.R") -enviro <- readRDS("data/enviro_grid20.RDS") +enviro <- readRDS("data/enviro_grid_sst15.RDS") enviro$dt <- 0.01 enviro$tmax <- 1000 @@ -35,4 +35,4 @@ zoomssgrid <- foreach(i=1:nrow(enviro), fZooMizer_run(groups = Groups, input = input, no_w = 177+1) } -saveRDS(zoomssgrid, file = "initial_zooplankton_20220104.RDS") \ No newline at end of file +saveRDS(zoomssgrid, file = "20230314_sst15_ZooMSS.RDS") \ No newline at end of file