v2.2.0 Age Structure and v6681 Code Update

.gitignore

@@ -5,7 +5,6 @@
 *.Rproj
 /Atlantis_Runs/
 /Run_Code/
-/Setup_Files/
 /Figures/
 /currentVersion/archive/
 
@@ -39,6 +38,7 @@ outputFolder
 # except csv
 !/currentVersion/neus_fisheries.csv
 !/currentVersion/neus_groups.csv
+!/currentVersion/neus_migrations.csv
 # except bgm
 !/currentVersion/neus_tmerc_RM2.bgm
 # except rmd
@@ -70,3 +70,8 @@ diagnostics/Obs_Hindcast_Group_Progress.xlsx
 #data files
 data/fishing_sensitivity_extended_constant_2/*
 
+
+#Setup Files for Cloud
+/Setup_Files/*
+!/Setup_Files/cloud*
+

R/Calibration_Tools/scale_mumC_age.R → R/Calibration_Tools/call_mumC_age.R

@@ -1,12 +1,12 @@
 #Script to scale Mum_C for age groups
 library(dplyr)
-source(here::here('R','edit_param_mum_age.R'))
-source(here::here('R','edit_param_C_age.R'))
+source(here::here('R','Calibration_Tools','edit_param_mum_age.R'))
+source(here::here('R','Calibration_Tools','edit_param_C_age.R'))
 
-group.names = 'RED'
+group.names = 'HER'
 
-mum.scale = 2
-C.scale = 1.5
+mum.scale = 1
+C.scale = 0.8
 
 bio.prm = here::here('currentVersion','at_biology.prm')
 
@@ -19,8 +19,10 @@ c.base.age = get_param_C_age(bio.prm)%>%
 
 group.id = which(group.names == mum.base.age[,1])
 
-new.mum = mum.base.age[group.id,2:11]*mum.scale
-new.C = c.base.age[group.id,2:11]*C.scale
+new.mum = as.numeric(mum.base.age[group.id,2:ncol(mum.base.age)])*mum.scale
+new.mum = new.mum[!is.na(new.mum)]
+new.C = as.numeric(c.base.age[group.id,2:ncol(mum.base.age)])*C.scale
+new.C = new.C[!is.na(new.C)]
 
 edit_param_mum_age(bio.prm ,
                    overwrite = T,

R/Calibration_Tools/compare_init_numbers.R

@@ -0,0 +1,44 @@
+#Script to compare initial numbers weight between reference and run
+library(ncdf4)
+dev.init = nc_open('/net/work3/EDAB/atlantis/dev_branch/currentVersion/neus_init.nc')
+master.init = nc_open('/net/work3/EDAB/atlantis/master_branch/currentVersion/neus_init.nc')
+new.6536.init = nc_open(here::here('currentVersion','neus_init_rescale_age.nc'))
+nre.6665.init = nc_open('/net/work3/EDAB/atlantis/Rob_proj/currentVersion/neus_init.nc')
+
+
+fgs.dev = read.csv('/net/work3/EDAB/atlantis/dev_branch/currentVersion/neus_groups.csv',as.is = T)
+fgs.master = read.csv('/net/work3/EDAB/atlantis/master_branch/currentVersion/neus_groups.csv',as.is = T)
+fgs.6536 = read.csv(here::here('currentVersion','neus_groups.csv'),as.is = T)
+fgs.6665 = read.csv('/net/work3/EDAB/atlantis/Rob_proj/currentVersion/neus_groups.csv',as.is = T)
+
+vert.names = fgs.dev$Name[which(fgs.dev$NumCohorts >2)]
+
+out.df = data.frame(Name = vert.names, dev = NA,master = NA, new.6536 = NA, new.6665 = NA)
+
+i=1
+for(i in 1:length(vert.names)){
+
+  ncohort.dev = fgs.dev$NumCohorts[which(fgs.dev$Name == vert.names[i])]
+  ncohort.master = fgs.master$NumCohorts[which(fgs.master$Name == vert.names[i])]
+  ncohort.6536 = fgs.6536$NumCohorts[which(fgs.6536$Name == vert.names[i])]
+  ncohort.6665 = fgs.6665$NumCohorts[which(fgs.6665$Name == vert.names[i])]
+
+  out.dev = lapply(1:ncohort.dev,function(x){
+    return(sum(ncvar_get(dev.init,paste0(vert.names[i],x,'_Nums')),na.rm=T))
+  })
+  out.master = lapply(1:ncohort.master,function(x){
+    return(sum(ncvar_get(master.init,paste0(vert.names[i],x,'_Nums')),na.rm=T))
+  })
+  out.6536 = lapply(1:ncohort.6536,function(x){
+    return(sum(ncvar_get(new.6536.init,paste0(vert.names[i],x,'_Nums')),na.rm=T))
+  })
+  out.6665 = lapply(1:ncohort.6665,function(x){
+    return(sum(ncvar_get(nre.6665.init,paste0(vert.names[i],x,'_Nums')),na.rm=T))
+  })
+
+  out.df$dev[i] = sum(unlist(out.dev))
+  out.df$master[i] = sum(unlist(out.master))
+  out.df$new.6536[i] = sum(unlist(out.6536))
+  out.df$new.6665[i] = sum(unlist(out.6665))
+
+}

R/Calibration_Tools/compare_weight_age.R

@@ -0,0 +1,112 @@
+#Script to compare weight at age for species between reference data and run
+#Determine if mum needs to be rescaled
+library(ggplot2)
+library(dplyr)
+
+#get run info
+run.name = '6536_new_age_param_init_rescale_2'
+run.dir = here::here('Atlantis_Runs',run.name,'')
+
+fgs = read.csv(here::here('currentVersion','neus_groups.csv'),as.is = T)%>%
+  select(Code,LongName,NumCohorts,NumAgeClassSize)
+fgs.old = read.csv(here::here('diagnostics','neus_groups_v2_0_1.csv'),as.is = T)%>%
+  select(Code,LongName,NumCohorts,NumAgeClassSize)
+fgs$spp.change= fgs.old$NumCohorts != fgs$NumCohorts
+
+#Read in reference survey length weigth age: indwt = kg, len = cm
+ref.stats.survdat = readRDS(here::here('data','survey_lenagewgt.rds'))%>%
+  group_by(Code,AGE)%>%
+  summarise(LENGTH = mean(LENGTH,na.rm=T),
+            WEIGHT = mean(INDWT,na.rm=T))%>%
+  rename(ref.len = 'LENGTH',
+         ref.wgt = 'WEIGHT')
+
+ref.stats.fishbase = read.csv(here::here('diagnostics','fishbase_len_wgt_ref.csv'),as.is = T)%>%
+  mutate(ref.wgt = Weight_g/1000)%>%
+  rename(ref.len = 'Length',
+         AGE = 'Age')%>%
+  select(Code,AGE,ref.len,ref.wgt)
+
+ref.stats.all = ref.stats.survdat %>%
+  bind_rows(ref.stats.fishbase)
+
+ref.spp = unique(ref.stats.all$Code)
+i=1
+ref.stats.spp.ls = list()
+for(i in 1:length(ref.spp)){
+
+  fgs.match = which(fgs$Code == ref.spp[i])  
+
+  spp.agecl = 1:fgs$NumCohorts[fgs.match]
+  spp.age = c(0,(1:fgs$NumCohorts[fgs.match])*fgs$NumAgeClassSize[fgs.match],Inf)
+
+  ref.stats.spp = ref.stats.all %>%
+    filter(Code == ref.spp[i] & !is.na(AGE))
+
+  spp.age.match = cut(ref.stats.spp$AGE+1,spp.age,labels = F,right =F)
+  spp.age.match[spp.age.match>fgs$NumCohorts[fgs.match]] = fgs$NumCohorts[fgs.match]
+
+  ref.stats.spp$agecl = spp.age.match
+  ref.stats.spp.ls[[i]] = ref.stats.spp
+}
+ref.stats.all2 = bind_rows(ref.stats.spp.ls)
+
+
+#Get run length weight stats
+run.wgt = readRDS(paste0(run.dir,'Post_Processed/Data/max_weight.rds'))%>%
+  group_by(species,agecl)%>%
+  summarise(maxMeanWeight = max(maxMeanWeight,na.rm=T)/1000)
+run.len = readRDS(paste0(run.dir,'Post_Processed/Data/length_age.rds'))%>%
+  group_by(species,agecl)%>%
+  summarise(length = mean(atoutput,na.rm=T))
+
+run.stats = run.wgt %>% 
+  left_join(run.len)%>%
+  left_join(fgs,by = c(species = 'LongName'))%>%
+  # mutate(AGE = agecl * NumAgeClassSize)%>%
+  rename(run.wgt = 'maxMeanWeight',
+         run.len = 'length')
+
+#Combine and create scalar based on length and weight
+combined.stats = run.stats%>%
+  left_join(ref.stats.all2)%>%
+  select(Code,species,spp.change,agecl,AGE,run.len,ref.len,run.wgt,ref.wgt)%>%
+  mutate(len.scale = ref.len / run.len,
+         wgt.scale = ref.wgt / run.wgt)
+
+mean.scale.wgt = combined.stats %>%
+  group_by(Code,spp.change,agecl)%>%
+  summarise(mum.scale = mean(wgt.scale,na.rm=T))%>%
+  # group_by(Code,spp.change)%>%
+  # summarise(mum.scale = mean(mum.scale,na.rm=T))%>%
+  filter(spp.change == T)
+  # filter(spp.change == T & is.na(mum.scale))
+
+source(here::here('R','Calibration_Tools','edit_param_mum_age.R'))
+mum.orig = get_param_mum_age(bio.prm = here::here('currentVersion','at_biology_rescale_mumC.prm'))
+new.mum = mum.orig
+
+spp.names = sort(unique(mean.scale.wgt$Code))
+
+j = 1
+for(j in 1:length(spp.names)){
+
+  which.mum = which(mum.orig$group == spp.names[j] )
+
+  mum.age.scale = filter(mean.scale.wgt,Code == spp.names[j]) %>% arrange(agecl)
+
+  mean.scale = mean(mum.age.scale$mum.scale,na.rm=T)
+
+  mum.scale.spp = mum.age.scale$mum.scale
+  mum.scale.spp[is.na(mum.scale.spp)] = mean.scale
+
+  mum.orig.spp = as.numeric(mum.orig[which.mum,2:ncol(mum.orig)])
+
+  new.mum[which.mum,2:ncol(new.mum)] = signif(mum.orig.spp * mum.scale.spp,2)
+
+}
+
+edit_param_mum_age(bio.prm = here::here('currentVersion','at_biology_rescale_mumC.prm'),
+                   single.group = F,
+                   new.mum = new.mum,
+                   overwrite = T)

R/Calibration_Tools/edit_param_BH.R

@@ -9,11 +9,15 @@ edit_param_BH = function(bio.prm,group.name,alpha,beta,overwrite = F, new.file.n
     alpha.line = grep(paste0('BHalpha_',group.name[i]),bio.lines)
     beta.line = grep(paste0('BHbeta_',group.name[i]),bio.lines)
 
-    new.alpha = paste0('BHalpha_',group.name[i],' ',alpha[i])
-    new.beta = paste0('BHbeta_',group.name[i],' ',beta[i])
+    if(!is.na(alpha)){
+      new.alpha = paste0('BHalpha_',group.name[i],' ',alpha[i])
+      bio.lines[alpha.line] = new.alpha
+    }
 
-    bio.lines[alpha.line] = new.alpha
-    bio.lines[beta.line] = new.beta
+    if(!is.na(beta)){
+      new.beta = paste0('BHbeta_',group.name[i],' ',beta[i])  
+      bio.lines[beta.line] = new.beta
+    }
   }
 
   if(overwrite){

R/Calibration_Tools/edit_param_C_age.R

@@ -9,15 +9,29 @@
 get_param_C_age = function(bio.prm, write.output = F, output.dir, out.name ){
 
   bio.lines = readLines(bio.prm)
-  bio.lines.id = grep('^C_.*10',bio.lines)
+  bio.lines.id = grep('^C_.*',bio.lines)
   bio.lines.vals1 = bio.lines[bio.lines.id]
+  which.invert = grepl('_T15',bio.lines.vals1)
+  bio.lines.vals1 = bio.lines.vals1[!which.invert]
+  bio.lines.id = bio.lines.id[!which.invert]
+
+  group.names =unname(sapply(bio.lines.vals1,function(x) strsplit(x,'C_|\t10.00|\t| ')[[1]][2]))
+  max.age = max(as.numeric(sapply(bio.lines.vals1,function(x) strsplit(x,'C_|\t| ')[[1]][3])))
+
+  C.mat = matrix(NA, nrow = length(group.names),ncol = max.age)
+  colnames(C.mat) = paste0('C',1:max.age)
+  out.df = data.frame(C.mat)
+  out.df = cbind(data.frame(group = group.names),out.df)
 
-  group.names =unname(sapply(bio.lines.vals1,function(x) strsplit(x,'C_|\t|10.00')[[1]][2]))
-  out.df = data.frame(group = group.names,
-                      C1 = NA, C2 = NA, C3 = NA, C4 = NA, C5 = NA, C6 = NA, C7 = NA, C8 = NA, C9 = NA, C10 = NA)
   for(i in 1:length(bio.lines.id)){
     C.group = bio.lines[bio.lines.id[i] + 1 ]
-    out.df[i,2:11] = strsplit(C.group,split = "\t| ")[[1]]
+    C.split = strsplit(C.group,split = "\t| |  ")[[1]]
+
+    if(length(C.split)>max.age){print(paste0(group.names[i],' has ',length(C.split)-10,' trailing tabs ',i))}
+
+    C.out = rep(NA,max.age)
+    C.out[1:length(C.split)] = C.split
+    out.df[i,2:ncol(out.df)] = C.out
   }
   if(write.output){
     write.csv(out.df, file = paste0(output.dir,out.name,'.csv'),row.names = F)
@@ -32,21 +46,33 @@ edit_param_C_age = function(bio.prm, new.C, overwrite = F,new.file.name,single.g
 
   #Get C_XXX bio.prm lines
   bio.lines = readLines(bio.prm)
-  bio.lines.id = grep('^C.*10',bio.lines)
+  bio.lines.id = grep('^C_',bio.lines)
   bio.lines.vals = bio.lines[bio.lines.id]
-  group.names =unname(sapply(bio.lines.vals,function(x) strsplit(x,'C_|\t10| 10.00')[[1]][2]))
+  which.invert = grepl('_T15',bio.lines.vals)
+  bio.lines.vals = bio.lines.vals[!which.invert]
+  bio.lines.id = bio.lines.id[!which.invert]
+
+  group.names =unname(sapply(bio.lines.vals,function(x) strsplit(x,'C_|\t10.00|\t| ')[[1]][2]))
+  max.age = max(as.numeric(sapply(bio.lines.vals,function(x) strsplit(x,'C_|\t| ')[[1]][3])))
 
   if(single.group){
 
     ind = which(group.name == group.names)
+    new.C = new.C[!is.na(new.C)]
     C.string = paste(new.C,collapse = '\t')
     bio.lines[bio.lines.id[ind]+1] = C.string
   }else{
     for(i in 1:nrow(new.C)){
 
-      ind = which(new.C$group[i] == group.names)
-      C.string = paste(new.C[i,2:11],collapse='\t')
+      ind = which(group.names == new.C$group[i])
+
+      C.string = new.C[i,2:ncol(new.C)]
+      which.na = which(is.na(C.string))
+      if(length(which.na > 0)){C.string = C.string[-which.na]}
+      C.string = paste(C.string,collapse='\t')
+
       bio.lines[bio.lines.id[ind]+1] = C.string
+
     }
   }
 

R/Calibration_Tools/edit_param_FSP.R

@@ -0,0 +1,47 @@
+#Changes FSP parameter in at_biology.prm
+
+
+get_param_FSP = function(bio.prm, group.name = NA){
+
+  bio.lines = readLines(bio.prm)
+
+  FSP.line = grep(paste0('^FSP_'),bio.lines)
+
+  groups = sapply(bio.lines[FSP.line],function(x) return(strsplit(x,'\t| |_')[[1]][2]),USE.NAMES = F)
+
+  group.vals = sapply(bio.lines[FSP.line],function(x) strsplit(x,' |\t')[[1]][2], USE.NAMES = F)
+
+  if(is.na(group.name)){
+
+    return(data.frame(Code = groups, FSP = group.vals))
+
+  }else{
+
+    return(group.vals[which(groups == group.name)])
+  }
+
+  return(out.df)
+}
+
+edit_param_FSP = function(bio.prm,group.name,value,unit = 'value',overwrite = F, new.file.name){
+
+  bio.lines = readLines(bio.prm)
+
+  FSP.line = grep(paste0('FSP_',group.name),bio.lines)
+
+  if(unit == 'value'){
+    new.val = value
+  }else{
+    old.val = as.numeric(strsplit(bio.lines[FSP.line],' |\t')[[1]][2])
+    new.val = old.val * value
+  }
+
+  bio.lines[FSP.line] = paste0('FSP_',group.name,' ',new.val)
+
+  if(overwrite){
+    writeLines(bio.lines, con = bio.prm)
+  }else{
+    file.copy(bio.prm, new.file.name, overwrite = T)
+    writeLines(bio.lines, con = new.file.name )
+  }
+}

R/Calibration_Tools/edit_param_FSPB.R

@@ -8,8 +8,8 @@ get_param_FSPB = function(bio.prm){
 
   groups = sapply(bio.lines[FSPB.line],function(x) return(strsplit(x,'\t| |_')[[1]][2]),USE.NAMES = F)
 
-  out.df = as.data.frame(matrix(NA,nrow = length(groups), ncol = 11))
-  colnames(out.df) = c('group',paste0('age.',1:10))
+  out.df = as.data.frame(matrix(NA,nrow = length(groups), ncol = 13))
+  colnames(out.df) = c('group',paste0('age.',1:12))
   out.df$group = groups
 
   for(i in 1:length(groups)){