Skip to content

Commit

Permalink
Merge pull request #244 from NEFSC/dev_branch
Browse files Browse the repository at this point in the history
v2.2.0 Age Structure and v6681 Code Update
  • Loading branch information
andybeet authored Mar 18, 2024
2 parents ebc792b + e139728 commit 1c94a94
Show file tree
Hide file tree
Showing 89 changed files with 24,581 additions and 21,918 deletions.
7 changes: 6 additions & 1 deletion .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,6 @@
*.Rproj
/Atlantis_Runs/
/Run_Code/
/Setup_Files/
/Figures/
/currentVersion/archive/

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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*

Original file line number Diff line number Diff line change
@@ -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')

Expand All @@ -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,
Expand Down
44 changes: 44 additions & 0 deletions R/Calibration_Tools/compare_init_numbers.R
Original file line number Diff line number Diff line change
@@ -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))

}
112 changes: 112 additions & 0 deletions R/Calibration_Tools/compare_weight_age.R
Original file line number Diff line number Diff line change
@@ -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)
12 changes: 8 additions & 4 deletions R/Calibration_Tools/edit_param_BH.R
Original file line number Diff line number Diff line change
Expand Up @@ -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){
Expand Down
44 changes: 35 additions & 9 deletions R/Calibration_Tools/edit_param_C_age.R
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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

}
}

Expand Down
47 changes: 47 additions & 0 deletions R/Calibration_Tools/edit_param_FSP.R
Original file line number Diff line number Diff line change
@@ -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 )
}
}
4 changes: 2 additions & 2 deletions R/Calibration_Tools/edit_param_FSPB.R
Original file line number Diff line number Diff line change
Expand Up @@ -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)){
Expand Down
Loading

0 comments on commit 1c94a94

Please sign in to comment.