Skip to content

Commit

Permalink
copy and source proj functions locally
Browse files Browse the repository at this point in the history
  • Loading branch information
@mkapur-noaa
mkapur-noaa committed Apr 10, 2024
1 parent 9ef3e55 commit 5d18cad
Show file tree
Hide file tree
Showing 5 changed files with 399 additions and 8 deletions.
49 changes: 41 additions & 8 deletions 2024/R/2024_analysis.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,11 +45,7 @@ prev_mdl_fldr = "18.2c_2020"

afscdata::bsai_fhs(year)
source(here::here(year,'r','bsai_fhs_wrangle_data.R'))

# run base model ----
## you will need to manually add the data to the DAT
## file in curr_mdl_fldr


# run retrospective ----
## NOT DOING THIS FOR AN UPDATE

Expand All @@ -59,8 +55,45 @@ source(here::here(year,'r','bsai_fhs_wrangle_data.R'))
# run projections ----
## takes less than one minute; only run this if model and/or projected catches change

# setwd(here::here(year,'model_runs','03b_projection'))
# shell('spm')
lapply(list.files(here::here(year,'r',"proj_functions/"), full.names = T, pattern = ".r$"), source)

# Write proj files ----

mod <- SS_output(here::here(year,'mgmt',curr_mdl_fldr), verbose = F)
## passed to write_proj function
NSEX=2 # number of sexes used in assessment model
Nfishery=1 # number of fisheries(fleets) #This was set equal to 2
fleets=1 # fleet index number (associated with commercial fishery)
rec_age=3 # assumed age at recruitment
max_age=21 # maximum age in model
NAGE=length(rec_age:max_age) # number of ages
FY=1964 # first year used to subset SSB, per memo this is always 1977, 1964 for consistency
rec_FY=1964 # first year used to subset recruitment
rec_LY_decrement=0 # value subtracted from assessment final year to subset recruitment vector
spawn_month=1 # spawning month
Fratios=1 # Proportion F per fishery
#passed to write_proj_spcat
ct_yrs=3 #Number of future catch years given to projection model
## passed to setup function
nsims=1000 # number of projection model simulations
nproj=14 # number of projection years ALSO USED BY get_proj_res
## passed to get_proj_res
spp="BSAI_flathead"

## this will give you the general values, but mightn't be perfectly formatted
## simply paste the right values into projection_data.dat
## and ensure the catches are specified in spm.dat
write_proj(dir=here::here(year,'model_runs','03b_projection'),
# sdir =x,
data_file="projection_data_dump.dat",
data= mod,
NSEX=NSEX, NAGE=NAGE, Nfishery=Nfishery,
fleets=fleets, rec_age=rec_age, max_age=max_age, FY=FY,
rec_FY=rec_FY, rec_LY_decrement=rec_LY_decrement,
spawn_month=spawn_month, Fratios=Fratios)

setwd(here::here(year,'model_runs','03b_projection'))
shell('spm')

rec_table1 <-
read.table(here::here(year,'model_runs','03b_projection','percentdb.out')) %>%
Expand Down Expand Up @@ -96,7 +129,7 @@ write.csv(rec_table, file = here::here(year,'model_runs','03b_projection',paste0

# process results ----

# create figures ----
# misc figures ----


## catch with inset
Expand Down
52 changes: 52 additions & 0 deletions 2024/R/proj_functions/get_proj_res.r
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,52 @@
#get projection results

get_proj_res<-function(dir1="foo",dir2="foo",proj_out="Model_Proj_out/",data_file="proj_res_summ.out",FY=2018, k=1,nproj=2,species="foo")
{

#Proj output files
alt2=read.table(file=file.path(dir1,proj_out,"alt2_proj.out"),header=TRUE)
print("Read alt2_proj.out")

#Check skip lines equation. Not sure if percentiles prints the same for everyone...it should though
ssb_ref=as.data.frame(matrix(scan(file=file.path(dir1,proj_out,"percentiles.out"),skip=(23+(nproj*5)),nlines=nproj),nrow=nproj,
ncol=9,byrow=TRUE))
print("Read ssb from percentiles .out")

names(ssb_ref)=c("Year", "SSB100", "SSBabc", "SSBofl", "LowCI_SSB", "Median_SSB", "Mean_SSB", "UpperCI_SSB", "Stdev_SSB")
print("attaching names to ssb_ref")
print("first row of ssb_ref")
print(ssb_ref[1,])

#See above comment
f_ref=as.data.frame(matrix(scan(file=file.path(dir1,proj_out,"percentiles.out"),skip=(26+(nproj*6)),nlines=nproj),nrow=nproj,ncol=9,byrow=TRUE))
print("Read f from percentiles .out")

names(f_ref)=c("Year", "F0", "Fabc", "Fofl", "LowCI_F", "Median_F", "Mean_F", "UpperCI_F", "Stdev_F")
print("attaching names to f_ref")
print("first row of f_ref")
print(f_ref[1,])


res=data.frame(Spp=species,
Year=round(alt2$Year[alt2$Year>=FY & alt2$Year<=(FY+15)],0),
TotBiom=round(alt2$TotBiom[alt2$Year>=FY & alt2$Year<=(FY+15)],0),
SSB=round(alt2$SSB[alt2$Year>=FY & alt2$Year<=(FY+15)],0),
SSB100=round(ssb_ref[ssb_ref$Year>=FY & ssb_ref$Year<=(FY+15),"SSB100"],0),
SSBabc=round(ssb_ref[ssb_ref$Year>=FY & ssb_ref$Year<=(FY+15),"SSBabc"],0),
SSBofl=round(ssb_ref[ssb_ref$Year>=FY & ssb_ref$Year<=(FY+15),"SSBofl"],0),
Fofl=round(f_ref$Fofl[f_ref$Year>=FY & f_ref$Year<=(FY+15)],2),
Fabc=round(f_ref$Fabc[f_ref$Year>=FY & f_ref$Year<=(FY+15)],2),
OFL=round(alt2$OFL[alt2$Year>=FY & alt2$Year<=(FY+15)],0),
ABC=round(alt2$ABC[alt2$Year>=FY & alt2$Year<=(FY+15)],0))
names(res)=c("Spp","Year","TotBiom","SSB","SSB100","SSBabc","SSBofl","Fofl","Fabc","OFL","ABC")
res$SSB_SSBofl=round(alt2$SSB[alt2$Year>=FY & alt2$Year<=(FY+15)],0)/round(ssb_ref[ssb_ref$Year>=FY & ssb_ref$Year<=(FY+15),"SSBofl"],0)

print("Summary of results")
print(res)

print("writing summary projection results")
write.table(res,file.path(dir1,data_file), row.names=FALSE, col.names=FALSE)
print("finished writing summary projection results")

return(res)
}
216 changes: 216 additions & 0 deletions 2024/R/proj_functions/write_proj.r
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,216 @@
write_proj<-function(dir="foo",
sdir="data/",
data_file="projection_data.dat",
data=Models[[1]],NSEX=2,NAGE=30,Nfishery=2,fleets=1:2
,rec_age=3, max_age=36,FY=1977,rec_FY=1977,rec_LY_decrement=2,
spawn_month=1, Fratios=0.5)
{
#Writing data file required by proj

LY=data$Retro_year # last year in the report file is the last year of full model;
# need this calc to get end year for the retros
Y5<-LY-5 ## start year to calculate 5-year mean apical F

rec_LY=LY-rec_LY_decrement ## last year of recruitment used to get mean rec

write_file=paste0(dir,"/",data_file)

T1<-noquote(paste(data_file))
write(T1,paste(write_file),ncolumns = 1 )

T1<-noquote(" 0 # SSL Species???")
write(T1,paste(write_file),ncolumns = 1,append=T)

T1<-noquote(" 0 # Constant Buffer Dorn")
write(T1,paste(write_file),append = T)

T1<-noquote(paste(Nfishery, "# Number of fisheries", sep=" "))
write(T1,paste(write_file),append = T)

T1<-noquote(paste(NSEX, "# Number of Sexes",sep=" "))
write(T1,paste(write_file),append = T)

T1<-noquote(paste(mean(data$sprseries$sum_Apical_F[data$sprseries$Yr>Y5 & data$sprseries$Yr<=LY]),"# 5 year average F ",sep=""))
write(T1,paste(write_file),append = T)

T1<-noquote("1 # Author f")
write(T1,paste(write_file),append = T)

T1<-noquote("0.4 # SPR ABC")
write(T1,paste(write_file),append = T)

T1<-noquote("0.35 # SPR MSY")
write(T1,paste(write_file),append = T)

T1<-noquote(paste(spawn_month, "# Spawning month",sep=" "))
write(T1,paste(write_file),append = T)

T1<-noquote(paste(NAGE,"# number of ages",sep=" "))
write(T1,paste(write_file),append = T)

T1<-noquote(paste(Fratios,"# Fratio",sep=" "))
write(T1,paste(write_file),append = T)

if(NSEX<2)
{
print("In NSEX <2")
T1<-noquote("# natural mortality at age females, males")
write(T1,paste(write_file),append = T)
print("Writing natural mortality")
if(rec_age==0)
{
print("In rec_age==0")
write(as.numeric(c(subset(data$M_at_age,data$M_at_age$Yr==(LY-1)& data$M_at_age$Sex==NSEX)[,as.character(seq(rec_age,(max_age-1),1))],
subset(data$M_at_age,data$M_at_age$Yr==(LY-1)&data$M_at_age$Sex==NSEX)[,as.character((max_age-1))])),
paste(write_file),append = T,ncolumns = 45)
}
if(rec_age>0){
print("In rec_age>0")
write(as.numeric(subset(data$M_at_age,data$M_at_age$Yr==(LY-1)&data$M_at_age$Sex==NSEX)[,as.character(seq((rec_age-1),(max_age-1),1))]),
paste(write_file),append = T,ncolumns = 45)
}

T1<-noquote("# Maturity females")
write(T1,paste(write_file),append = T)
print("Writing maturity")
if(spp=="GT")write(rep(1,NAGE),paste(write_file),append = T,ncolumns = 45)
if(spp!="GT") write(round(as.numeric(subset(data$endgrowth,data$endgrowth$Sex==NSEX & data$endgrowth$int_Age %in% rec_age:max_age)[,"Age_Mat"]),4),paste(write_file),append = T,ncolumns = 45)

T1<-noquote("# wt spawn females")
write(T1,paste(write_file),append = T)
print("Writing wt spawn females")
if(spp=="GT")write(round(as.numeric(subset(data$endgrowth,data$endgrowth$Sex==NSEX & data$endgrowth$int_Age %in% rec_age:max_age)[,"Mat*Fecund"]),4),paste(write_file),append = T,ncolumns = 45)
if(spp!="GT")write(round(as.numeric(subset(data$endgrowth,data$endgrowth$Sex==NSEX & data$endgrowth$int_Age %in% rec_age:max_age)[,"Wt_Beg"]),4),paste(write_file),append = T,ncolumns = 45)

T1<-noquote("# WtAge females by fishery")
write(T1,paste(write_file),append = T)
print("Writing WtAge by sex and fishery")
for (i in 1:length(fleets)){
write(round(as.numeric(subset(data$ageselex,data$ageselex$Fleet==fleets[i] & data$ageselex$Yr==LY &
data$ageselex$Factor=="bodywt")[NSEX,as.character(seq(rec_age,max_age,1))]),4),
paste(write_file),append = T,ncolumns = 45)

}

T1<-noquote("# Selectivity by fishery")
write(T1,paste(write_file),append = T)
print("Writing Selectivity by sex and fishery")
for (i in 1:length(fleets))
{
write(round(as.numeric(subset(data$ageselex,data$ageselex$Fleet==fleets[i] & data$ageselex$Yr==LY & data$ageselex$Factor=="Asel2")[NSEX,as.character(seq(rec_age,max_age,1))]),4),paste(write_file),append = T,ncolumns = 45)
}

T1<-noquote(paste("# Numbers at age",rec_age,"-",max_age,"females males",LY))
write(T1,paste(write_file),append = T)
print("Writing Numbers at age by sex")
write(as.numeric(subset(data$natage,data$natage[,"Beg/Mid"]=="B"&data$natage$Yr==LY)[NSEX,as.character(seq(rec_age,max_age,1))]),paste(write_file),append = T,ncolumns = 45)

Rec_1<-as.numeric(data$natage[,as.character(rec_age)][data$natage$Yr<=rec_LY & data$natage$Yr>=rec_FY & data$natage$Sex==(NSEX) & data$natage[,"Beg/Mid"]=="B"])
N_rec<-length(Rec_1)
T1<-noquote("# No Recruitments")
write(T1,paste(write_file),append = T)
print("Writing number of recruit observations")
write(N_rec,paste(write_file),append = T,ncolumns = 45)

T1<-noquote("# Recruitment")
write(T1,paste(write_file),append = T)
print("Writing recruits")
write(round(Rec_1,1),paste(write_file),append = T,ncolumns = 45)

}

if(NSEX>1)
{
T1<-noquote("# natural mortality")
write(T1,paste(write_file),append = T)
print("Writing natural mortality")
for(i in 1:NSEX)
{
if(rec_age==0)
{
write(as.numeric(c(subset(data$M_at_age,data$M_at_age$Yr==(LY-1)&data$M_at_age$Sex==NSEX)[,as.character(seq(rec_age,(max_age-1),1))],
subset(data$M_at_age,data$M_at_age$Yr==(LY-1)&data$M_at_age$Sex==NSEX)[,as.character((max_age-1))])),
paste(write_file),append = T,ncolumns = 45)
}
if(rec_age>0){
write(as.numeric(subset(data$M_at_age,data$M_at_age$Yr==(LY-1)&data$M_at_age$Sex==NSEX)[,as.character(seq((rec_age-1),(max_age-1),1))]),
paste(write_file),append = T,ncolumns = 45)
}
}

T1<-noquote("# Maturity females")
write(T1,paste(write_file),append = T)
print("Writing female maturity")
if(spp=="GT")write(rep(1,NAGE),paste(write_file),append = T,ncolumns = 45)
if(spp!="GT") write(round(as.numeric(subset(data$endgrowth,data$endgrowth$Sex==1 & data$endgrowth$int_Age %in% rec_age:max_age)[,"Age_Mat"]),4),paste(write_file),append = T,ncolumns = 45)

T1<-noquote("# Maturity males")
write(T1,paste(write_file),append = T)
print("Writing male maturity")
if(spp=="GT")write(rep(1,NAGE),paste(write_file),append = T,ncolumns = 45)
if(spp!="GT") write(round(as.numeric(subset(data$endgrowth,data$endgrowth$Sex==2 & data$endgrowth$int_Age %in% rec_age:max_age)[,"Age_Mat"]),4),paste(write_file),append = T,ncolumns = 45)

T1<-noquote("# wt spawn females")
write(T1,paste(write_file),append = T)
print("Writing wt spawn females")
if(spp=="GT")write(round(as.numeric(subset(data$endgrowth,data$endgrowth$Sex==1 & data$endgrowth$int_Age %in% rec_age:max_age)[,"Mat*Fecund"]),4),paste(write_file),append = T,ncolumns = 45)
if(spp!="GT")write(round(as.numeric(subset(data$endgrowth,data$endgrowth$Sex==1 & data$endgrowth$int_Age %in% rec_age:max_age)[,"Wt_Beg"]),4),paste(write_file),append = T,ncolumns = 45)

T1<-noquote("# WtAge by sex and fishery")
write(T1,paste(write_file),append = T)
print("Writing wtAge by sex and fishery")
for(j in 1:NSEX)
{
for (i in 1:length(fleets))
{
write(round(as.numeric(subset(data$ageselex,data$ageselex$Fleet==fleets[i] & data$ageselex$Yr==LY & data$ageselex$Factor=="bodywt")[j,as.character(seq(rec_age,max_age,1))]),4),paste(write_file),append = T,ncolumns = 45)
}
}

T1<-noquote("# Selectivity by sex and fishery")
write(T1,paste(write_file),append = T)
print("Writing selectivity by sex and fishery")
for(j in 1:NSEX)
{
for (i in 1:length(fleets))
{
write(round(as.numeric(subset(data$ageselex,data$ageselex$Fleet==fleets[i] & data$ageselex$Yr==LY & data$ageselex$Factor=="Asel2")[j,as.character(seq(rec_age,max_age,1))]),4),paste(write_file),append = T,ncolumns = 45)
}
}

T1<-noquote(paste("# Numbers at age",rec_age,"-",max_age,"females males",LY))
write(T1,paste(write_file),append = T)
print("Writing Numbers at age females males")
for(j in 1:NSEX)
{
write(as.numeric(subset(data$natage,data$natage[,"Beg/Mid"]=="B"&data$natage$Yr==LY)[j,as.character(seq(rec_age,max_age,1))]),paste(write_file),append = T,ncolumns = 45)
}

print("Getting number of recruits")
Rec_1<-as.numeric(data$natage[,as.character(rec_age)][data$natage$Yr<=rec_LY & data$natage$Yr>=rec_FY & data$natage$Sex==(NSEX-1) & data$natage[,"Beg/Mid"]=="B"]
+data$natage[,as.character(rec_age)][data$natage$Yr<=rec_LY & data$natage$Yr>=rec_FY & data$natage$Sex==NSEX & data$natage[,"Beg/Mid"]=="B"])
print("Getting number of recruit obs")
N_rec<-length(Rec_1)
T1<-noquote("# No Recruitments")
write(T1,paste(write_file),append = T)
print("Writing number of recruit obs")
write(N_rec,paste(write_file),append = T,ncolumns = 45)

T1<-noquote(paste("# Recruitment age ",rec_age,"+ ", rec_FY,"-",rec_LY,sep=""))
write(T1,paste(write_file),append = T)
print("Writing number of recruits")
write(round(Rec_1,1),paste(write_file),append = T,ncolumns = 45)


}

T1<-noquote(paste("# SSB ", FY,"-",LY,sep=""))
write(T1,paste(write_file),append = T)
print("Writing SSB")
if(spp!="EBS_Pcod") write(as.numeric(data$sprseries$SSB[data$sprseries$Yr<=LY&data$sprseries$Yr>=FY]),paste(write_file),append = T,ncolumns = 45)
if(spp=="EBS_Pcod") write(as.numeric(data$sprseries$SSB[data$sprseries$Yr<=LY&data$sprseries$Yr>=FY])/2,paste(write_file),append = T,ncolumns = 45)
# file.copy(write_file, file.path(dir, sdir),overwrite=TRUE)

}


29 changes: 29 additions & 0 deletions 2024/R/proj_functions/write_setup.r
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,29 @@
write_setup <-function(dir="foo",sdir="data/",data_file="setup.dat",data=Models[[1]],nsims=100,nproj=20)
{
LY=data$Retro_year

write_file=paste0(dir,"/",data_file)
#write set-up.dat required to run projection
write(paste(noquote("std"), noquote("#No idea what this is"),sep=" "),file=write_file)
write(paste(7, noquote("#Number of alternatives???"),sep=" "),file=write_file,append=TRUE)
write(paste(1, noquote("#Vector of alt numbers"),sep=" "),file=write_file,append=TRUE)
write(2,file=write_file,append=TRUE)
write(3,file=write_file,append=TRUE)
write(4,file=write_file,append=TRUE)
write(5,file=write_file,append=TRUE)
write(6,file=write_file,append=TRUE)
write(7,file=write_file,append=TRUE)
write(paste(1,noquote("#TAC_ABC"),sep=" "),file=write_file,append=TRUE)
write(paste(1,noquote("#SRType"),sep=" "),file=write_file,append=TRUE)
write(paste(1,noquote("#Rec_Gen"),sep=" "),file=write_file,append=TRUE)
write(paste(0,noquote("#Fmsy_F35"),sep=" "),file=write_file,append=TRUE)
write(paste(0,noquote("#Rec_Cond"),sep=" "),file=write_file,append=TRUE)
write(paste(1,noquote("#Write_Big"),sep=" "),file=write_file,append=TRUE)
write(paste(nproj,noquote("#NProj"),sep=" "),file=write_file,append=TRUE)
write(paste(nsims,noquote("#Nsims"),sep=" "),file=write_file,append=TRUE)
write(paste(LY,noquote("#Styr"),sep=" "),file=write_file,append=TRUE)

#create copy in data folder
file.copy(write_file, file.path(dir, sdir),overwrite=TRUE)
}

Loading

0 comments on commit 5d18cad

Please sign in to comment.