From 8fac63559cef33dcf92b5e6f0969a603eddcfa42 Mon Sep 17 00:00:00 2001 From: Jim Ianelli Date: Tue, 20 Feb 2024 19:55:55 -0800 Subject: [PATCH] Step one to start organizgin working examples --- DESCRIPTION | 2 +- R/readData.R | 190 ++-- examples/atka/plot.R | 39 - examples/atka/spm.dat | 2 +- examples/atka/spm.tpl | 2509 ----------------------------------------- 5 files changed, 119 insertions(+), 2623 deletions(-) delete mode 100644 examples/atka/plot.R delete mode 100644 examples/atka/spm.tpl diff --git a/DESCRIPTION b/DESCRIPTION index ec65c1a..c91a805 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -13,5 +13,5 @@ Suggests: knitr, ggplot2, rmarkdown VignetteBuilder: knitr -RoxygenNote: 7.1.2 +RoxygenNote: 7.3.0 URL: http://afsc-assessments.github.io/spmR/ diff --git a/R/readData.R b/R/readData.R index e0d0b8d..247032c 100644 --- a/R/readData.R +++ b/R/readData.R @@ -1,116 +1,160 @@ -#' list2dat +#' Write List Object to Projection Model Data File #' -#' write list object to projection model data file +#' This function writes a list object to a file formatted for a projection model. +#' Each element of the list is written with a header. #' -#' @param D objection to be written to -#' @return written data file for spm model +#' @param D List object to be written. +#' @param fn Filename where the data will be written. +#' @param hdr Header text to be included in the file. +#' @return The function does not return a value; it writes to a file. #' @export -list2dat <- function(D,fn,hdr="a new file") { +#' @examples +#' # Example usage: +#' # list2dat(myList, "datafile.dat", "Header for new file") +list2dat <- function(D, fn, hdr="a new file") { + # Open file connection + sink(fn) + on.exit(sink()) # Ensure the connection is closed when the function exits + cat(paste0("# ", hdr, "\n")) + + for (i in seq_along(D)) { + cat(paste0("#", names(D[i]), "\n")) + write.table(D[[i]], append = TRUE, quote = FALSE, row.names = FALSE, col.names = FALSE) + } # The following writes a data file - cat(file=fn,paste0("# ",hdr,"\n")) - ol <-length(D) - for (i in 1:ol){ - cat(file=fn,paste0("#",names(D[i]),"\n"),append=TRUE) - write.table(D[[i]],file=fn,append=TRUE,quote=FALSE,row.names=FALSE,col.names=FALSE) - } + # cat(file=fn,paste0("# ",hdr,"\n")) + # ol <-length(D) + # for (i in 1:ol){ + # cat(file=fn,paste0("#",names(D[i]),"\n"),append=TRUE) + # write.table(D[[i]],file=fn,append=TRUE,quote=FALSE,row.names=FALSE,col.names=FALSE) + # } } -#' dat2list -#' Read list object to projection model data file +#' Read Data File into List Object +#' +#' Reads a data file formatted for a projection model and converts it into a list object. +#' The function handles numeric and NA values appropriately. #' -#' @return update from last year's files +#' @param fn Filename of the data file to be read. +#' @return A list object containing the data from the file. #' @export -dat2list <- function(fn) -{ - options(warn=-1) #Suppress the NA message in the coercion to double - ifile=scan(fn,what="character",flush=TRUE,blank.lines.skip=FALSE,quiet=TRUE) - #idx=sapply(as.double(ifile),is.na) - idx=substr(ifile,1,1)=="#" - vnam=ifile[idx] #list names - #vnam - nv=length(vnam) #number of objects - A=list() - ir=0 - for(i in 1:nv) - { - ir=match(vnam[i],ifile) - ##print(ir) - #print(vnam[i]) - if(i!=nv) irr=match(vnam[i+1],ifile) else irr=length(ifile)+1 #next row - dum=NA - if(irr-ir==2) dum=as.double(scan(fn,skip=ir,nlines=1,quiet=TRUE,what="")) - if(irr-ir>2) dum=as.matrix(read.table(fn,skip=ir,nrow=irr-ir-1,fill=TRUE)) - - if(is.numeric(dum))#Logical test to ensure dealing with numbers - { - A[[substr(vnam[i],2,10)]]=dum - } - if(is.na(dum))#Logical test to ensure dealing with numbers - { - A[[substr(vnam[i],2,10)]]=scan(fn,skip=ir,nlines=1,quiet=TRUE,what="") - } - } - options(warn=0) - - return(A) +#' @examples +#' # Example usage: +#' # myDataList <- dat2list("datafile.dat") +dat2list <- function(fn) { + options(warn = -1) # Suppress warnings temporarily + on.exit(options(warn = 0)) # Reset warning options on exit + + file_content = scan(fn, what = "character", flush = TRUE, blank.lines.skip = FALSE, quiet = TRUE) + header_indices = substr(file_content, 1, 1) == "#" + header_names = file_content[header_indices] + + list_data = list() + for (i in seq_along(header_names)) { + start_line = match(header_names[i], file_content) + end_line = if (i < length(header_names)) match(header_names[i + 1], file_content) - 1 else length(file_content) + + data_chunk = file_content[(start_line + 1):end_line] + if (length(data_chunk) == 1) { + list_data[[substr(header_names[i], 2)]] = as.numeric(data_chunk) + } else { + list_data[[substr(header_names[i], 2)]] = read.table(text = data_chunk, fill = TRUE) + } + } + + return(list_data) } -#' print_Tier3_tables +#' Print Tier 3 Tables +#' +#' Generates and prints HTML tables for Tier 3 projections including catch, ABC, fishing mortality, and spawning biomass for various scenarios. #' -#' @param df +#' @param df Data frame containing Tier 3 data. +#' @param modname Name of the model used, defaults to "base". +#' @param stock Name of the stock, defaults to "BSAI Atka mackerel". +#' @return HTML tables for the specified Tier 3 data. #' @export -#' -print_Tier3_tables <- function(df, modname="base",stock="BSAI Atka mackerel") { - library(xtable) +#' @importFrom xtable xtable +#' @importFrom dplyr select, group_by, summarise, spread +#' @examples +#' # Example usage: +#' # print_Tier3_tables(myDataFrame, "model1", "Some Fish Stock") +print_Tier3_tables <- function(df, modname="base", stock="BSAI Atka mackerel") { tabcap<-tablab <- c("tier3_C","tier3_ABC","tier3_F","tier3_SSB") tabcap[1]=paste0("Tier 3 projections of ",stock," catch for the 7 scenarios.") tabcap[2]=paste0("Tier 3 projections of ",stock," ABC for the 7 scenarios.") tabcap[3]=paste0("Tier 3 projections of ",stock," fishing mortality for the 7 scenarios.") tabcap[4]=paste0("Tier 3 projections of ",stock," spawning biomass for the 7 scenarios.") - + # Stock Alt Sim Yr SSB Rec Tot_biom SPR_Implied F Ntot Catch ABC OFL AvgAge AvgAgeTot SexRatio FABC FOFL bfsum <- df %>% select(Alt,Yr,SSB,F,ABC ,Catch) %>% group_by(Alt,Yr) %>% summarise(Catch=mean(Catch),SSB=mean(SSB),F=mean(F),ABC=mean(ABC)) - - tC <- bfsum %>% select(Alt,Yr,Catch) %>% spread(Alt,Catch) + + tC <- bfsum %>% select(Alt,Yr,Catch) %>% spread(Alt,Catch) names(tC) <- c("Catch","Scenario 1","Scenario 2","Scenario 3","Scenario 4","Scenario 5","Scenario 6","Scenario 7") - - tB <- bfsum %>% select(Alt,Yr,SSB) %>% spread(Alt,SSB) + + tB <- bfsum %>% select(Alt,Yr,SSB) %>% spread(Alt,SSB) names(tB) <- c("SSB","Scenario 1","Scenario 2","Scenario 3","Scenario 4","Scenario 5","Scenario 6","Scenario 7") - - tF <- bfsum %>% select(Alt,Yr,F) %>% spread(Alt,F) + + tF <- bfsum %>% select(Alt,Yr,F) %>% spread(Alt,F) names(tF) <- c("F","Scenario 1","Scenario 2","Scenario 3","Scenario 4","Scenario 5","Scenario 6","Scenario 7") - - tA <- bfsum %>% select(Alt,Yr,ABC) %>% spread(Alt,ABC) + + tA <- bfsum %>% select(Alt,Yr,ABC) %>% spread(Alt,ABC) names(tA) <- c("ABC","Scenario 1","Scenario 2","Scenario 3","Scenario 4","Scenario 5","Scenario 6","Scenario 7") - + tab <- (data.frame(tC)) rownames(tab)<-c() cap <- tabcap[1] - for (i in 2:length(tab[1,]) ) - tab[,i] <- formatC((tab[,i]), format="d", big.mark=",") + for (i in 2:length(tab[1,]) ) + tab[,i] <- formatC((tab[,i]), format="d", big.mark=",") tab <- xtable(tab, caption = cap, label=paste0("tab:",tablab[1]), digits=0, auto=TRUE, align=rep("r",(length(tab[1,])+1)) ) print(tab, "html", caption.placement = "top",include.rownames = FALSE, sanitize.text.function = function(x){x}, scalebox=.85) - + tab <- (data.frame(tB)) cap <- tabcap[2] - for (i in 2:length(tab[1,]) ) - tab[,i] <- formatC(as.numeric(tab[,i]), format="d", big.mark=",") + for (i in 2:length(tab[1,]) ) + tab[,i] <- formatC(as.numeric(tab[,i]), format="d", big.mark=",") tab <- xtable(tab, caption = cap, label=paste0("tab:",tablab[2]),digits=0, auto=TRUE, align=rep("r",(length(tab[1,])+1)) ) print(tab, "html", caption.placement = "top",include.rownames = FALSE, sanitize.text.function = function(x){x}, scalebox=.85) - + tab <- (data.frame(tF)) cap <- tabcap[3] - for (i in 2:length(tab[1,]) ) - tab[,i] <- formatC(as.numeric(tab[,i]), format="f",digits=3) + for (i in 2:length(tab[1,]) ) + tab[,i] <- formatC(as.numeric(tab[,i]), format="f",digits=3) tab <- xtable(tab, caption = cap, label=paste0("tab:",tablab[3]), digits=3, align=rep("r",(length(tab[1,])+1)) ) print(tab, "html", caption.placement = "top",include.rownames = FALSE, sanitize.text.function = function(x){x}, scalebox=.85) - + tab <- (data.frame(tA)) cap <- tabcap[4] - for (i in 2:length(tab[1,]) ) - tab[,i] <- formatC(as.numeric(tab[,i]), format="d", big.mark=",") + for (i in 2:length(tab[1,]) ) + tab[,i] <- formatC(as.numeric(tab[,i]), format="d", big.mark=",") tab <- xtable(tab, caption = cap, label=paste0("tab:",tablab[4]),digits=0, auto=TRUE, align=rep("r",(length(tab[1,])+1)) ) print(tab, "html", caption.placement = "top",include.rownames = FALSE, sanitize.text.function = function(x){x}, scalebox=.85) return(tab) +# +# +# if (!requireNamespace("dplyr", quietly = TRUE)) stop("dplyr package is required but not installed") +# if (!requireNamespace("xtable", quietly = TRUE)) stop("xtable package is required but not installed") +# +# # Compute summaries +# summary_df <- df %>% +# dplyr::select(Alt, Yr, SSB, F, ABC, Catch) %>% +# dplyr::group_by(Alt, Yr) %>% +# dplyr::summarise(Catch = mean(Catch), SSB = mean(SSB), F = mean(F), ABC = mean(ABC)) +# +# # Prepare and print tables +# table_types <- c("Catch", "SSB", "F", "ABC") +# for (type in table_types) { +# formatted_table <- create_formatted_table(summary_df, type, stock) +# print_table(formatted_table, type, stock) +# } +# } +# +# create_formatted_table <- function(df, type, stock) { +# # Table creation logic +# } +# +# print_table <- function(table, type, stock) { +# # Table printing logic +#} } diff --git a/examples/atka/plot.R b/examples/atka/plot.R deleted file mode 100644 index 8555737..0000000 --- a/examples/atka/plot.R +++ /dev/null @@ -1,39 +0,0 @@ -radian -rm(list=ls()) -ls() -#source(paste0(here::here("R"),"/prelims.R")) -source("../../R/prelims.R") -#------------------------------------------------------------------------------- -# Visual compare runs -#------------------------------------------------------------------------------- -library(ggridges) - -source("../../R/compareRuns.r") -#--Projections--------- -pdf("proj.pdf") -pfn <- "amak" -pfn <- "5yr" -pdt <- read_csv(paste0("spm_detail.csv")) -pdt$Alternative <- as.factor(pdt$Alternative) -names(pdt) - -pt <- pdt%>%filter(Yr>2021,Alternative==2) %>% group_by(Yr,Alternative) %>% summarise(Catch=mean(Catch),ABC=mean(ABC),OFL=mean(OFL),SSB=median(SSB) ,lb=quantile(SSB,.2) ,ub=quantile(SSB,.8) ) -pt -ggplot(pt,aes(x=Yr,y=SSB,fill=Alternative)) + geom_line() + mytheme + ylim(c(0,340000)) + geom_ribbon(aes(ymin=lb,ymax=ub,fill=Alternative),alpha=0.25) + labs(y="Spawning biomass (kt)",x="Year") + scale_x_continuous(breaks=seq(2022,2035,2)) -c1 <- ggplot(pt,aes(x=Yr,y=Catch,color=Alternative,size=1.)) + geom_line(size=1.5) + mytheme + labs(y="Catch (kt)",x="Year") + scale_x_continuous(breaks=seq(2015,2032,2)) -c1 <- c1 + geom_line(aes(x=Yr,y=ABC),size=1) -#c1 <- c1 + geom_line(data=pt[as.numeric(Alternative)==2,.(Yr,ABC)],aes(x=Yr,y=ABC)) -c1 -pt[as.numeric(Alternative)==2,.(Yr,ABC)] -pt <- pdt[Yr>2018,.(Catch=mean(Catch),ABC=mean(ABC),OFL=mean(OFL)),.(Yr,Alternative)] -pt -ggplot(pt,aes(x=Yr,y=OFL,color=Alternative)) + geom_line() + mytheme -pdt -pdx <-rbind(pdt) -setkey(pdx,Yr,Alternative) - -pt <- pdx[.(Yr>2018,(Alternative)==1),.(Catch=mean(Catch),ABC=mean(ABC),OFL=mean(OFL),SSB=median(SSB) ),.(Yr,config)] -pt <- pdx[Yr>2018&Alternative==1,.(Catch=mean(Catch),ABC=mean(ABC),OFL=mean(OFL),SSB=median(SSB) ,lb=quantile(Catch,.1) ,ub=quantile(Catch,.9) ),.(Yr,config)] -ggplot(pt,aes(x=Yr,y=Catch,color=config,shape=config)) + geom_line() + geom_point() + mytheme + geom_ribbon(aes(ymin=lb,ymax=ub,fill=config),alpha=0.25) + labs(x="Year")+ scale_x_continuous(breaks=seq(2015,2032,2)) + ylim(c(0,130)) -dev.off() - diff --git a/examples/atka/spm.dat b/examples/atka/spm.dat index d350203..89a5b38 100644 --- a/examples/atka/spm.dat +++ b/examples/atka/spm.dat @@ -34,7 +34,7 @@ std # Run name 5 scenarios # OY Max 2.00E+06 # data files for each species -../amak.prj +amak.prj # ABC Multipliers 1 # scalars diff --git a/examples/atka/spm.tpl b/examples/atka/spm.tpl deleted file mode 100644 index b791fc1..0000000 --- a/examples/atka/spm.tpl +++ /dev/null @@ -1,2509 +0,0 @@ - // DONT FORGET ABOU RHO!!! - //////////////////////////////////////////////////////////////////////////// - // spm.tpl - // Version 0.6 Released Oct 2022 - // - // NOTE: Draft subject to change - // Conventions: - // nspp, ispp = index for species/stock - // ngear, igear= index for gears used within age-structured species/stocks - // - // Oct 2005 add ability to compute objective function for maximizing yield while minimizing variability - // want to have the ability to compute the average first-diff squared and avg yield for each trajectory - // Alt3b_obj_fun = avg_yield - 0.5 * sqrd_1st_diff - // - //////////////////////////////////////////////////////////////////////////// -DATA_SECTION - !!CLASS ofstream means_out("means.out") - !!CLASS ofstream alts_proj("alt_proj.out") - !!CLASS ofstream percent_out("percentiles.out") - !!CLASS ofstream percent_db("percentdb.out") - // !!CLASS ofstream Alt3bstuff("alt3b.out") - !!CLASS ofstream detail_out("spm_detail.csv") - !!CLASS ofstream prof_F("F_profile.out"); - !!CLASS ofstream elasticity("elasticity.csv"); - int condition_SR - int ipro - int isim - number kink_adj - !! condition_SR = 0;// - int rnseed // Random number seed - !! rnseed = 123; - //////////////////////////////////////////////////////////////////////////// - // !! ad_comm::change_datafile_name("setup.dat"); Jim changed to spm.dat file - !! *(ad_comm::global_datafile) >> run_name; // Read in the name of this run - init_int Tier - int alt ; // Alternative specfications (relic of PSEIS) - vector rec_vector(1,300) // storage vector for historical and future recruitment... - vector wtd_rec(1,25) - vector wtd_div(1,25) // denominator of wtd recruitment - init_int nalts - init_ivector alt_list(1,nalts) - init_int TAC_ABC // Flag to set TAC equal to ABC (1 means true, otherwise false) - init_int SrType // Type of recruitment curve (1=Ricker, 2 Bholt) - // Specify form of recruitment generator (1 = use observed mean and std - // 2 = use estimated SRR and estimated sigma R - // 3 = use input parameters from srecpar.dat : - // 4 = use input parameters from srecpar.dat : - init_int Rec_Gen - init_int Fmsy_F35 // Specify if conditioned so that Fmsy = F35 (affects SRR fitting) may need to be species specific... - init_number Rec_Cond // Specify prior condition that recruitment at half and double average SSB is similar to average historical Rec - init_int Write_Big // Flag to write big file (of all simulations rather than a summary, 0 means don't do it, otherwise do it) - init_int npro // Number of projection years - init_int nsims // Number of simulaions - init_int styr // First year of projection - !! cout<< "First year:\t"<> spp_file_name(i); - - write_log(nyrs_catch_in); - write_log(nspp); - write_log(OY_min); - write_log(OY_max); - write_log(spp_file_name); - - END_CALCS - !! cout <<"OYMax "<> bzero_in; - *(ad_comm::global_datafile) >> phizero_in; - *(ad_comm::global_datafile) >> alpha_in; - *(ad_comm::global_datafile) >> sigmar_in; - *(ad_comm::global_datafile) >> rho_in; - // rho_in=0.82; - } - // Open up tac-model parameters - ad_comm::change_datafile_name("tacpar.dat"); - END_CALCS - init_int nntmp - init_int nnodes - init_vector maxabc(1,ntacspp) - init_matrix theta(0,nnodes,1,ntacspp) - !! cout<<"read tacpar"<> spname(i); // 1 - cout<> SSL_spp(i); // 2 - *(ad_comm::global_datafile) >> Const_Buffer(i); // 3 - *(ad_comm::global_datafile) >> ngear(i); // 4 - *(ad_comm::global_datafile) >> nsexes(i); // 6 - *(ad_comm::global_datafile) >> avg_5yrF(i); // 7 - *(ad_comm::global_datafile) >> FABC_Adj(i); // 8 - *(ad_comm::global_datafile) >> SPR_abc(i); // 9 - *(ad_comm::global_datafile) >> SPR_ofl(i); // 10 - *(ad_comm::global_datafile) >> spawnmo(i); // 11 - *(ad_comm::global_datafile) >> nages(i); // 12 - cout<<"nages: "<> Fratiotmp(i,j); // 13 - write_log( Fratiotmp(i)); - - for (int k=1;k<=nages(i);k++) - *(ad_comm::global_datafile) >> M_Ftmp(i,k); // 14 - if (nsexes(i)==2) - for (int k=1;k<=nages(i);k++) - *(ad_comm::global_datafile) >> M_Mtmp(i,k); // 15 - write_log( M_Ftmp(i)); - - for (int k=1;k<=nages(i);k++) - *(ad_comm::global_datafile) >> pmaturetmp_F(i,k); // 16 - write_log( pmaturetmp_F(i)); - - if (nsexes(i)==2) - for (int k=1;k<=nages(i);k++) - *(ad_comm::global_datafile) >> pmaturetmp_M(i,k); // 15 - write_log( pmaturetmp_M(i)); - cout << "Mature: "<< pmaturetmp_F(i)(1,nages(i)) <> wt_Ftmp(i,k); // 17 - write_log( wt_Ftmp(i)); - for (int j=1;j<=ngear(i);j++) - for (int k=1;k<=nages(i);k++) - *(ad_comm::global_datafile) >> wt_gear_Ftmp(i,j,k); // 18 - write_log( wt_gear_Ftmp(i)); - - if (nsexes(i)==2) - for (int j=1;j<=ngear(i);j++) - for (int k=1;k<=nages(i);k++) - *(ad_comm::global_datafile) >> wt_gear_Mtmp(i,j,k);// 19 - write_log( wt_gear_Mtmp(i)); - - for (int j=1;j<=ngear(i);j++) - for (int k=1;k<=nages(i);k++) - *(ad_comm::global_datafile) >> sel_Ftmp(i,j,k); // 20 - write_log( sel_Ftmp(i)); - - if (nsexes(i)==2) - for (int j=1;j<=ngear(i);j++) - for (int k=1;k<=nages(i);k++) - *(ad_comm::global_datafile) >> sel_Mtmp(i,j,k); // 21 - write_log( sel_Mtmp(i)); - for (int k=1;k<=nages(i);k++) - *(ad_comm::global_datafile) >> n0_Ftmp(i,k); // 22 - write_log( n0_Ftmp(i)); - - if (nsexes(i)==2) - for (int k=1;k<=nages(i);k++) - *(ad_comm::global_datafile) >> n0_Mtmp(i,k); // 23 - write_log( n0_Mtmp(i)); - - cout<<"N: "<> nrec(i); // 24 - cout<<"nrec: "<> Rtmp(i,j); // 25 - cout<<"Rec: "<> SSBtmp(i,j); // 26 - cout<<"SSB: "<> Rsim(ispp,i,j); - Rsim(ispp,i,j) *= exp(rnorms(i,j)*.375); // about 15% CV to get historical mean - } - if (nsims<=5) Rsim(ispp,i,j) = AMeanRec(ispp); // XXX constant recruitment - } - } - envin.close(); - AMeanRec(ispp) *= .5; // Arithmetic mean - cout <<"recruits"< 0 ) - F_begin_yr(k,ispp) = SolveF2(n0_F(ispp), n0_M(ispp), Obs_Catch(k,ispp),ispp); // F_yr_one(ispp) = SolveF2(n0_F(ispp), n0_M(ispp), yr_one_catch(ispp),ispp); - else - F_begin_yr(ispp) = 0; // F_yr_one(ispp) = 0; // cout< 1) // Debugging if LP to be done or not - { - int on=0; - if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-nolp"))>-1) - dolp=0; - else - dolp=1; - } - END_CALCS - - -PRELIMINARY_CALCS_SECTION - double tmp1; - write_alts_hdr(); - get_SB100(); - cout<<"SSB, Biomass at unfished: "<1) - nyrs_catch = nyrs_catch_in; - else - nyrs_catch = 1;// NOTA BUENO: this is changed under the new EIS Alternatives (May 06) - Do_Sims(); - // if (alt==2) write_alts(); - write_alts(); - } - -FUNCTION compute_obj_fun - dvariable tmp1; - dvariable tmp2; - tmp1.initialize(); - for (int ispp=1;ispp<=nspp;ispp++) - { - Get_Bzero(ispp); - // write_srec();exit(1); - if(Fmsy_F35>0) - { - get_msy(ispp); - switch (current_phase()) - { - case 1 : tmp1 = 1.e1*square(log(Fmsy(ispp))-log(F35(ispp))); break; - case 2 : tmp1 = 1.e2*square(log(Fmsy(ispp))-log(F35(ispp))); break; - case 3 : tmp1 = 1.e2*square(log(Fmsy(ispp))-log(F35(ispp))); break; - case 4 : tmp1 = 1.e3*square(log(Fmsy(ispp))-log(F35(ispp))); break; - default: tmp1 = 1.e3*square(log(Fmsy(ispp))-log(F35(ispp))); break; - } - if(Fmsy_F35==2) - switch (current_phase()) - { - case 1 : tmp1 += 1.e1*square(log(Bmsy(ispp))-log(0.35 * SB100(ispp))); break; - case 2 : tmp1 += 1.e2*square(log(Bmsy(ispp))-log(0.35 * SB100(ispp))); break; - case 3 : tmp1 += 1.e2*square(log(Bmsy(ispp))-log(0.35 * SB100(ispp))); break; - case 4 : tmp1 += 1.e3*square(log(Bmsy(ispp))-log(0.35 * SB100(ispp))); break; - default : tmp1+= 1.e3*square(log(Bmsy(ispp))-log(0.35 * SB100(ispp))); break; - } - obj_fun += tmp1; - } - if (Rec_Cond>0.) - { - // More conditioning here--make mean recruitment consistent with half and double - // recruitment at avg spawning biomass... - double vartmp = 2.*Rec_Cond*Rec_Cond; - /* - double ssb1 = AMeanSSB(ispp) * 0.5; - double ssb2 = AMeanSSB(ispp) ; - double ssb3 = AMeanSSB(ispp) * 2.0; - // double ssb4 = value(Bzero(ispp)) ; - tmp2.initialize(); - tmp2 += square(log(AMeanRec(ispp)) - log(SRecruit( ssb1, ispp)))/ vartmp; - tmp2 += square(log(AMeanRec(ispp)) - log(SRecruit( ssb2, ispp)))/ vartmp; - tmp2 += square(log(AMeanRec(ispp)) - log(SRecruit( ssb3, ispp)))/ vartmp; - // tmp2 += 24.* square(log(AMeanRec(ispp)) - log(SRecruit( ssb4, ispp))); - // tmp2 += 24.* square(log(AMeanSSB(ispp)) - log( ssb4 )); - */ - tmp2 = square(log(Bzero(ispp)) - log(SB100(ispp)))/ vartmp; - obj_fun += tmp2; - // cout << "SSB "<OY_max) - Actual_Catch = TAC/sum(TAC)*OY_max; - else - Actual_Catch = TAC; - } - /* ////////////////////////////////////////////////////////////////////////////// //int use_max = 1; // Need to move this into the setup file... double max_catch=1500.;// EBS POllock special case (CHANGE THIS) //if (use_max==1) for (int ispp=1;ispp<=nspp;ispp++) Actual_Catch(ispp) = min(TAC(ispp),max_catch); // cout<<"AC: "<2) - write_sim("Alternative ",ispp); - } - write_by_time(); - - -FUNCTION Alt4_TAC - // Return vector of Alt4 TACs (given Alt4_Fabc, and the SSL prey condition to be at least up to SSL MaxPerm) - // Get all ABCs under Alt4_Fabc - // Sum and test if below OY_min - double sumabc; - sumabc = sum(ABC); - if (sumabc3) - { - Alt4_Fcalc = 1; - double diff; - double ssl_sum=0.; - diff = OY_min-sumabc; - // if so, sum up SSL ABCs - ssl_sum = SSL_spp*ABC; - for (int ispp=1;ispp<=nspp;ispp++) - { - // then apportion them to difference between current ABC and target (OY_min) difference - if (SSL_spp(ispp)) - { - TAC(ispp) = ABC(ispp) + diff * ABC(ispp)/ssl_sum; - // Special for BSAI Pcod to subtract off 3%...and get the totals to match up - // if (ispp==2) TAC(ispp) /= 0.97; For EIS work - } - else - { - TAC(ispp) = ABC(ispp); - } - } - // cout << ipro<<" "<0.01) ftmp = SolveF2(N_F(ispp),N_M(ispp),TAC(ispp),ispp); - // cout <= 0.2 *SBzero(ispp) & SBtmp < SBKink(ispp) ) // NOTE Same as Am 56 here (until 20% bzero reached) - Ftmp = (Fabc(ispp)*(1/(1-alpha ))*(SBtmp/SBKink(ispp) - alpha )); - if (SBtmp > SBKink(ispp) ) - Ftmp = (Fabc(ispp)); - SBtmp = N_females * elem_prod( wt_mature_F(ispp), mfexp( -yrfrac(ispp)*(M_F(ispp) + Ftmp*F_age ))); - } - return(Ftmp); - -FUNCTION double Get_F_Am56(const dvector& F_age, const dvector& N_females, const int ispp ) - double Ftmp; // cout<< spname(ispp)<< endl; - { - for (ii=1;ii<=3;ii++) // Iterate to get month of spawning correct - { - if (SBtmp < alpha*SBKink(ispp)) - Ftmp = 0.; - if (SBtmp >= alpha*SBKink(ispp) & SBtmp < SBKink(ispp) ) - Ftmp = (Fabc(ispp)*(1/(1-alpha))*(SBtmp/SBKink(ispp) - alpha)); - if (SBtmp > SBKink(ispp) ) - Ftmp = (Fabc(ispp)); - SBtmp = N_females * elem_prod( wt_mature_F(ispp), mfexp( -yrfrac(ispp)*(M_F(ispp) + Ftmp*F_age ))); - } - } - return(Ftmp); - -FUNCTION double Get_Fofl_t2(const dvector& F_age, const dvector& N_females, const int ispp ) - double Ftmp; // cout<< spname(ispp)<< endl; - { - for (ii=1;ii<=3;ii++) // Iterate to get month of spawning correct - { - if (SBtmp < alpha*SBKink(ispp)) - Ftmp = 0.; - if (SBtmp >= alpha*SBKink(ispp) & SBtmp < SBKink(ispp) ) - Ftmp = (Fofl(ispp)*(1/(1-alpha))*(SBtmp/SBKink(ispp) - alpha)); - if (SBtmp > SBKink(ispp) ) - Ftmp = (Fofl(ispp)); - SBtmp = N_females * elem_prod( wt_mature_F(ispp), mfexp( -yrfrac(ispp)*(M_F(ispp) + Ftmp*F_age ))); - } - } - return(Ftmp); - -FUNCTION double Get_Fofl_t(const dvector& F_age, const dvector& N_females, const int ispp ) - double Ftmp; // cout<< spname(ispp)<< endl; - { - for (ii=1;ii<=3;ii++) // Iterate to get month of spawning correct - { - if (SBtmp < alpha*SBKink(ispp)) - Ftmp = 0.; - if (SBtmp >= alpha*SBKink(ispp) & SBtmp < SBKink(ispp) ) - Ftmp = (Fofl(ispp)*(1/(1-alpha))*(SBtmp/SBKink(ispp) - alpha)); - if (SBtmp > SBKink(ispp) ) - Ftmp = (Fofl(ispp)); - SBtmp = N_females * elem_prod( wt_mature_F(ispp), mfexp( -yrfrac(ispp)*(M_F(ispp) + Ftmp*F_age ))); - } - } - return(Ftmp); - - -FUNCTION double Get_F_t(const dvector& F_age, const dvector& N_females, const int ispp ) - double Ftmp; // cout<< spname(ispp)<< endl; - if (SSL_spp(ispp)) - { - Ftmp = Get_F_SSL_prey(F_age, N_females, ispp); - } - else - { - Ftmp = Get_F_Am56(F_age, N_females, ispp); - } - return(Ftmp); - -FUNCTION void Project_Pops(const int& isim, const int& i) - double ctmp; - if ( i == npro && isim%int(nsims/4)==0) cout << "Year "< 0) - { - if (i <= nyrs_catch || TAC_ABC==0) // Use TAC setting algorithm for alt 2 only, for all others, set TAC==ABC - { - Ftmp = SolveF2(N_F(ispp),N_M(ispp),Actual_Catch(ispp),ispp); - } - else - { - if (alt==7 && i<=3) - Ftmp = Get_F(1,ispp); // Set to the F rather than solving every time... - else - if (alt==77 && i<=2) - Ftmp = Get_F(1,ispp); // Set to the F rather than solving every time... - else - { - // if (alt==2 && sum(TAC)>OY_max) - if ((alt==2 ||alt==98||alt==97) && sum(TAC)>OY_max) - { - Ftmp = SolveF2(N_F(ispp),N_M(ispp),OY_max*TAC(ispp)/sum(TAC),ispp); - // cout<<"HERE "<< Ftmp<<" "< 0.) - { - for (m=1;m<=ngear(ispp);m++) - { - Ftottmp_F = Ftmp*Frat(ispp,m)*sel_F(ispp,m); - Ftottmp_spr += Ftottmp_F ; - Ftottmp_M = Ftmp*Frat(ispp,m)*sel_M(ispp,m); - ctmp += (wt_gear_F(ispp,m) * elem_prod(elem_div(Ftottmp_F, Z_F(ispp)),elem_prod(1.-S_F(ispp),N_F(ispp)))); // Catch equation (vectors) - ctmp += (wt_gear_M(ispp,m) * elem_prod(elem_div(Ftottmp_M, Z_M(ispp)),elem_prod(1.-S_M(ispp),N_M(ispp)))); // Catch equation (vectors) - } - SPRsim(ispp,isim,i) = Implied_SPR(Ftottmp_spr,ispp); - // cout <NUL "); - ifstream tac("tac.dat"); - tac>>TAC; - tac.close(); /* */ - for (int ispp=1;ispp<=nspp;ispp++) - Actual_Catch(ispp) = min(TAC(ispp),ABC(ispp)); - } - - -FUNCTION Avg_Age - Avg_Age_End.initialize(); - Avg_Age_Mat.initialize(); - for (int ispp=1;ispp<=nspp;ispp++) - { - // if(!isit_const(ispp) ) - { - dvector age_seq(1,nages(ispp)); - age_seq.initialize(); - age_seq.fill_seqadd(1,1); - dvector ntmp(1,nages(ispp)); - ntmp = N_F(ispp); - Avg_Age_End(ispp) += age_seq * ntmp /sum(ntmp); - Avg_Age_sum(ispp) += Avg_Age_End(ispp) ; - ntmp = elem_prod(ntmp,pmature_F(ispp)) ; - Avg_Age_Mat(ispp) += elem_prod(ntmp,pmature_F(ispp)) * age_seq/sum(ntmp); - // cout< 1e-6) - { - iter++; - ftmp += (TACin-cc) / btmp; - Ftottmp_F.initialize(); - Ftottmp_M.initialize(); - for (m=1;m<=ngear(ispp);m++) - { - Ftottmp_F += ftmp*Fratsel_F(m); - Ftottmp_M += ftmp*Fratsel_M(m); - } - Z_F = Ftottmp_F + M_F(ispp); - Z_M = Ftottmp_M + M_M(ispp); - S_F = mfexp( -Z_F ); - S_M = mfexp( -Z_M ); - cc = 0.0; - for (m=1;m<=ngear(ispp);m++) - { - cc += (wt_gear_F(ispp,m) * elem_prod(elem_div(ftmp*Fratsel_F(m), Z_F),elem_prod(1.-S_F,N_F))); // Catch equation (vectors) - cc += (wt_gear_M(ispp,m) * elem_prod(elem_div(ftmp*Fratsel_M(m), Z_M),elem_prod(1.-S_M,N_M))); // Catch equation (vectors) - } - dd = cc / TACin - 1.; - //cout << ispp<<" "<< ftmp << " "<< cc << " "<100) {cerr<<"Bombed on catch solver--check scales for "<nages(ispp)) - { - Tg += double(ii) * wt_mature_F(ispp,nages(ispp)) * ntmp; - tmp += wt_mature_F(ispp,nages(ispp)) * ntmp; - ntmp *= exp(-M_F(ispp,nages(ispp))); - } - else - { - Tg += double(ii) * wt_mature_F(ispp,ii) * ntmp; - tmp += wt_mature_F(ispp,ii) * ntmp; - ntmp *= exp(-M_F(ispp,ii)); - } - } - // Tg /= value(phizero(ispp)); - Tg /= tmp; - report << Tg<< " "; - } - report << endl; - report << "Options SR_Type Project_Recr_Assmption SR_Condition"< 1e-6) - means_out << mean_value <<" "; - else - means_out << " NA "; - } - means_out << endl; - } - means_out << endl; - -FUNCTION void write_TACs(const adstring& Title) - // This one prints out species over time (means only), but without headings - means_out <<"Alternative "< 1e-6) - means_out << mean_value <<" "; - else - means_out << " NA "; - } - means_out << endl; - } - means_out << endl; - - - -FUNCTION void write_sim(const adstring& Title, d3_array& Outtmp) - // This one prints out species over time (means only), but without headings - d3_array mtmp(1,nspp,1,npro,1,nsims); - for (int ispp=1;ispp<=nspp;ispp++) for (int i=1;i<=npro;i++) for (int k=1;k<=nsims;k++) - mtmp(ispp,i,k)=Outtmp(ispp,k,i); - means_out <<"Alternative "< 1e-6) - means_out << mean_value <<" "; - else - means_out << " NA "; - } - means_out << endl; - } - means_out << endl; - -FUNCTION void write_sim(const adstring& Title,const d3_array& Outtmp,const dvar_vector& bb) - // This one prints out species over time (means only) - d3_array mtmp(1,nspp,1,npro,1,nsims); - for (int ispp=1;ispp<=nspp;ispp++) for (int i=1;i<=npro;i++) for (int k=1;k<=nsims;k++) - mtmp(ispp,i,k)=Outtmp(ispp,k,i); - means_out <<"Alternative "<1e-6) - means_out << mean_value <<" "; - else - means_out << " NA "; - } - means_out << endl; - } - means_out << endl; - -FUNCTION void write_sim(const adstring& Title,const int& ispp) - // For each species separately - percent_out <<"Alternative "< - adstring xspname; - adstring_array targsppname(1,20); - adstring_array spp_file_name(1,20); - adstring_array gearname(1,8); - adstring_array spname(1,90); - adstring_array areaname(1,8); - adstring run_name; - - ofstream write_log("Input_Log.rep"); - #undef write_log - #define write_log(object) write_log << #object "\n" << object << endl; - - // A routine to get transpose, sort and return a matrix ---- - dmatrix TranSort (const dmatrix m1) - { - RETURN_ARRAYS_INCREMENT(); - dmatrix vtmp=trans(m1); - int npro=m1.colmax(); - int nsim=m1.rowmax(); - for (int i=1;i<=npro;i++) - vtmp(i) = sort(vtmp(i),nsim); - - RETURN_ARRAYS_DECREMENT(); - return(vtmp); - } - // #include - -FUNCTION dvariable SRecruit(const dvariable& Stmp,const int& ispp) - RETURN_ARRAYS_INCREMENT(); - dvariable RecTmp; - switch (SrType) - { - case 1: - RecTmp = (Stmp / phizero(ispp)) * mfexp( sr_alpha(ispp) * ( 1. - Stmp / Bzero(ispp) )) ; //Ricker form from Dorn - break; - case 2: - RecTmp = Stmp / ( sr_alpha(ispp) + beta(ispp) * Stmp); //Beverton-Holt form - break; - case 3: - // RecTmp = mfexp(log_avgrec); //Avg recruitment - break; - case 4: - RecTmp = Stmp * mfexp( sr_alpha(ispp) - Stmp * beta(ispp)) ; //old Ricker form - break; - } - - RETURN_ARRAYS_DECREMENT(); - return RecTmp; - - -FUNCTION dvar_vector SRecruit(const dvar_vector& Stmp,const int& ispp) - RETURN_ARRAYS_INCREMENT(); - dvar_vector RecTmp(Stmp.indexmin(),Stmp.indexmax()); - switch (SrType) - { - case 1: - RecTmp = elem_prod((Stmp / phizero(ispp)) , mfexp( sr_alpha(ispp) * ( 1. - Stmp / Bzero(ispp) ))) ; //Ricker form from Dorn - break; - case 2: - RecTmp = elem_prod(Stmp , 1. / ( sr_alpha(ispp) + beta(ispp) * Stmp)); //Beverton-Holt form - break; - case 3: - // RecTmp = mfexp(log_avgrec); //Avg recruitment - break; - case 4: - RecTmp = elem_prod(Stmp ,mfexp( sr_alpha(ispp) - Stmp * beta(ispp))); //old Ricker form - break; - } - RETURN_ARRAYS_DECREMENT(); - //cout <5||F1<0.01)) - // { - // ii=5; - // F1=M(ispp,nages(ispp)); // When things bomb (F <0 or F really big then just set it to F35...) - // } - // else - { - F2 = F1 + df*.5; - F3 = F2 - df; //F1 = double(ii)/100; - yld1 = yield(F1,Stmp,Rtmp,ispp); //cout <