diff --git a/2023_runs/.README.md.swo b/2023_runs/.README.md.swo deleted file mode 100644 index d89600d..0000000 Binary files a/2023_runs/.README.md.swo and /dev/null differ diff --git a/2023_runs/addatsage/pm.dat b/2023_runs/addatsage/pm.dat index e457392..4363a4a 100644 --- a/2023_runs/addatsage/pm.dat +++ b/2023_runs/addatsage/pm.dat @@ -1,5 +1,5 @@ Updated_2022_ats_age -../dat/pm_23.2.dat +../dat/pm_23.3.dat ../dat/selvar23.dat control.dat ../dat/pm_fmsy_alt.dat diff --git a/2023_runs/addavo23/pm b/2023_runs/addavo23/pm deleted file mode 100755 index 94221cd..0000000 Binary files a/2023_runs/addavo23/pm and /dev/null differ diff --git a/2023_runs/addavo23/pm b/2023_runs/addavo23/pm new file mode 120000 index 0000000..5134594 --- /dev/null +++ b/2023_runs/addavo23/pm @@ -0,0 +1 @@ +../../src/pm \ No newline at end of file diff --git a/2023_runs/addavo23/pm.dat b/2023_runs/addavo23/pm.dat index a810fdb..eab7b07 100644 --- a/2023_runs/addavo23/pm.dat +++ b/2023_runs/addavo23/pm.dat @@ -1,5 +1,5 @@ Updated_2023_avo_value -../dat/pm_23.1.dat +../dat/pm_23.2.dat ../dat/selvar23.dat control.dat ../dat/pm_fmsy_alt.dat diff --git a/2023_runs/addavo23/pm.tpl b/2023_runs/addavo23/pm.tpl deleted file mode 100644 index f6a77bf..0000000 --- a/2023_runs/addavo23/pm.tpl +++ /dev/null @@ -1,6648 +0,0 @@ -///////////////////////////////////////////////////////////// -// Conventions: // -// fsh - refers to fshery // -// bts - refers to trawl survey // -// ats - refers to hydroacoustic survey // -// eac - refers to expected age composition // -// oac - refers to observed age composition // -// et - refers to expected total numbers // -// ot - refers to observed total numbers // -/////////////////////////////////////////////////// -// ngears - number of gear types (including srv \\ -// n_ - refers to the number of observations \\ -// styr - first calendar year of model \\ -// endyr - last calendar year of model \\ -// styr_ - first year of data (suffix) \\ -// endyr_r - last year for retrospective run \\ -//\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ -// Aug 25 2009 add projection to next year to get survey biomass estimates -// -DATA_SECTION - !! *(ad_comm::global_datafile) >> model_name; - !! *(ad_comm::global_datafile) >> datafile_name; - !! *(ad_comm::global_datafile) >> selchng_filename; - !! *(ad_comm::global_datafile) >> control_filename; - !! *(ad_comm::global_datafile) >> Alt_MSY_File; - !! *(ad_comm::global_datafile) >> Cov_Filename; - !! *(ad_comm::global_datafile) >> Wtage_file; - !! *(ad_comm::global_datafile) >> RawSurveyCPUE_file; - // !! *(ad_comm::global_datafile) >> control_temp_pred_filename; - !! *(ad_comm::global_datafile) >> Temp_Cons_Dist_file; - !! *(ad_comm::global_datafile) >> GenGamm_Filename; - // !! *(ad_comm::global_datafile) >> endyrn_file; - !! write_log(model_name); - !! write_log(datafile_name); - !! write_log(selchng_filename); - !! write_log(control_filename); - !! write_log(Alt_MSY_File); - !! write_log(Cov_Filename); - !! write_log(Wtage_file); - !! write_log(RawSurveyCPUE_file); - // !! write_log(control_temp_pred_filename); - !! write_log(Temp_Cons_Dist_file); - !! write_log(GenGamm_Filename); - // !! write_log(endyrn_file); - int count_Ffail; - int count_mcmc; - int count_mcsave; - !! count_mcmc=0; - !! count_mcsave=0; - int pflag; - int do_EIT1; - int q_amin - int q_amax - !! q_amin = 3; q_amax= 15; // age range overwhich q applies (for prior specifications) - vector selages(1,15) - !! selages=1.0;selages(1)=0;selages(2)=0; - vector sel_avo_in(1,15) - !! sel_avo_in(1)=0.0; sel_avo_in(2)=1; sel_avo_in(3)=1; sel_avo_in(4)=0.85; sel_avo_in(5)=0.7; sel_avo_in(6)=0.55; sel_avo_in(7)=0.3; sel_avo_in(8)=0.15; sel_avo_in(9)=0.05; sel_avo_in(10)=0.01; sel_avo_in(11)=0.01; sel_avo_in(12)=0.01; sel_avo_in(13)=0.01; sel_avo_in(14)=0.01; sel_avo_in(15)=0.01; - vector Cat_Fut(1,10); - !! do_EIT1=1; // flag to carry EIT out in future year (for simulations only) - !! pflag=0; - !!CLASS ofstream srecout("srec_Ass_out.rep") - !!CLASS ofstream projout("pm.prj") - !!CLASS ofstream nofish("nofish.rep") - !!CLASS ofstream projout2("pmsr.prj") - // writes when -mceval invoked (which uses pm.psv file from -mcmc 10000000 -mcsave 200) - !!CLASS ofstream eval("pm_eval.rep") - // Control file read from here-------------------------------------- - !! ad_comm::change_datafile_name(control_filename); cout<<"Opening "<0) // Use logistic selectivity - { - phase_selcoffs_fsh = -2; - if(phase_seldevs_fsh>0) - phase_selcoffs_fsh_dev = -phase_seldevs_fsh; - else - phase_selcoffs_fsh_dev = phase_seldevs_fsh; - - phase_logist_fsh_dev = phase_seldevs_fsh; - } - else // Use coefficient selectivities... - { - if(phase_seldevs_fsh>0) - phase_logist_fsh_dev = -phase_seldevs_fsh; - - phase_selcoffs_fsh_dev = phase_seldevs_fsh; - } - // Trawl Survey selectivity................ - if (phase_logist_bts>0) // Use logistic selectivites... - { - phase_selcoffs_bts = -2; - if(phase_seldevs_bts>0) - phase_selcoffs_bts_dev = -phase_seldevs_bts; - else - phase_selcoffs_bts_dev = phase_seldevs_bts; - - phase_logist_bts_dev = phase_seldevs_bts; - } - else // Use coefficient selectivites... - { - if(phase_seldevs_bts>0) - phase_logist_bts_dev = -phase_seldevs_bts; - - phase_selcoffs_bts_dev = phase_seldevs_bts; - } - - END_CALCS - init_int phase_selcoffs_ats - init_int phase_selcoffs_ats_dev - !! cout <<"Phase fsh coef: "<1978) - styr_est = styr; - else - styr_est = 1978; - } - else - { - styr_est = styr+1; - } - - // Set bounds for SRR params (depending on estimation and type) - if (phase_sr>0) - { - phase_nosr=-1; - phase_Rzero=3; - if (SrType==3) - phase_Rzero=-2; - } - else - { - phase_nosr=1; - phase_Rzero=-1; - } - - if (SrType==4) - Steepness_UB=4.; - else - Steepness_UB=.99; - - n_fsh_r=0; - n_bts_r=0; - n_ats_r=0; - n_avo_r=0; - endyr_r = endyr - int(ctrl_flag(28)); // where retrospective peels are set - endyr_est = endyr_r - int(ctrl_flag(29)); // lop off last couple of years - cout <<"Last yr of estimation..."<0) {ii++;} nch_ats=ii; - ii = 0; for(i=styr;i<=endyr_r;i++) if(fsh_ch_in(i)>0) {ii++;} nch_fsh=ii; - if (ctrl_flag(28)==0.) - { - n_fsh_r=n_fsh; - n_bts_r=n_bts; - n_ats_r=n_ats; - n_avo_r=n_avo; - } - else - { - int ii=0; - for (i=styr;i0) - phase_nosr = -1; - else - phase_nosr = 1; - END_CALCS - !! ad_comm::change_datafile_name(Cov_Filename);cout<<"Opening "<0) {ii++;sel_ch_sig_ats(ii)=ats_ch_in(i);yrs_ch_ats(ii)=i;} - ivector yrs_ch_fsh(1,nch_fsh); - vector sel_ch_sig_fsh(1,nch_fsh); - !! ii=0;for (i=styr;i<=endyr_r;i++) if(fsh_ch_in(i)>0) {ii++;sel_ch_sig_fsh(ii)=fsh_ch_in(i);yrs_ch_fsh(ii)=i;} - // Critical variables: yrs_ch_ats, sel_ch_sig_ats, and nch_ats - !! write_log(nch_fsh); write_log(yrs_ch_fsh); write_log(sel_ch_sig_fsh); - !! write_log(nch_ats); write_log(yrs_ch_ats); write_log(sel_ch_sig_ats); - matrix oac_fsh(1,n_fsh_r,1,nbins)//--Observed proportion at age in Fishery - matrix oac_bts(1,n_bts_r,1,nbins)//--Observed proportion at age in Survey - - // Set up whether last survey value is used or not (if ALK is from BTS instead of EIT) - int n_ats_ac_r; - !! if (use_last_ats_ac>0) n_ats_ac_r = n_ats_r; else n_ats_ac_r = n_ats_r-1; - !! if (use_last_ats_ac>0) nagecomp(3) = n_ats_r; else nagecomp(3) = n_ats_r-1; - - matrix oac_ats(1,n_ats_ac_r,1,nbins)//--Observed proportion at age in Survey - vector oa1_ats(1,n_ats_ac_r) //--Observed age 1 index from hydro survey - vector ot_fsh(1,n_fsh_r) //--Observed total numbers in Fishery - vector ot_bts(1,n_bts_r) //--Observed total numbers in BTS - vector std_ob_bts(1,n_bts_r) //--Observed total biomass in BTS - vector ob_bts(1,n_bts_r) //--Observed total biomass in BTS - - vector ot_ats(1,n_ats_r) //--Observed total numbers in HydroSurvey - vector ob_ats(1,n_ats_r) //--Observed total biomass in Survey - vector std_ob_ats(1,n_ats_r) //--Observed total biomass in Survey - int ignore_last_ats_age1; - vector lseb_ats(1,n_ats_r); - vector lvarb_ats(1,n_ats_r); - - // Stuff for passing seed via command line and for initializing new sequence - int iseed; - !! iseed=0; - int do_fmort; - !! do_fmort=0; - vector adj_1(1,10) - vector adj_2(1,10) - vector SSB_1(1,10) - vector SSB_2(1,10) - int do_check // Placeholder to have debug checker flag... - LOCAL_CALCS - do_check=0; - if (ad_comm::argc > 1) - { - int on=0; - if ( (on=option_match(argc,argv,"-io"))>-1) - do_check = 1; - if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-iseed"))>-1) - { - if (on>ad_comm::argc-2 | ad_comm::argv[on+1][0] == '-') - { - cerr << "Invalid number of iseed arguements, command line option -iseed ignored" << endl; - } - else - { - iseed = atoi(ad_comm::argv[on+1]); - cout<< "Currently using "<-1) - do_fmort=1; - } - // Some defaults------------ - adj_1=1.0; - adj_2=1.0; - SSB_1=1.0; - SSB_2=1.0; - - END_CALCS - int phase_cpue_q - int phase_avo_q - !!long int lseed=iseed; - !!CLASS random_number_generator rng(iseed); - - // ------------------------------------------------------ - // Conditionally read in generalized Gamma Q values for BTS - matrix GenGamData(1,n_bts,1,2); - vector q_GenGam(1,n_bts); - vector sd_GenGam(1,n_bts); - - LOCAL_CALCS - // Flag to use covariance for bottom trawl survey, 0=use vector, 1=use cov matrix, 2=use GenGamma - if (DoCovBTS==2) - { - ad_comm::change_datafile_name(GenGamm_Filename);cout<<"Opening "<> GenGamData ; - } - q_GenGam = extract_column(GenGamData,1); - sd_GenGam = extract_column(GenGamData,2); - write_log(q_GenGam); - write_log(sd_GenGam); - write_log(GenGamData); //exit(1); - END_CALCS - // init_matrix GenGamData(1,n_bts,1,2); - - // ------------------------------------------------------ - // read in wt-age data by - !! ad_comm::change_datafile_name(Wtage_file);cout<<"Opening "<1) phase_d_scale = 3; else phase_d_scale = -1; - !! write_log(endyr_wt); - !! write_log(styr_wt); - !! write_log(log_sd_coh); - !! write_log(log_sd_yr); - !! write_log(wt_obs); - !! write_log(sd_obs); - !! endyr_wt = endyr_wt - int(ctrl_flag(28)); - int max_nyrs_data; - !! cout<<(endyr_wt)<0) - { - mnCPUE(iyr,ist) = mean(CPUE(iyr,ist)); - mntemp(iyr,ist) = mean(temp(iyr,ist)); - } - } - } - END_CALCS - vector FW_fsh(1,4); - number FW_bts; - number FW_ats; - - -// read in data for doing temp-recuitment model, and spatial predation - !! ad_comm::change_datafile_name(Temp_Cons_Dist_file); - init_vector SST(styr-1,endyr-1); - number SST_fut - !! SST_fut = mean(SST(endyr-7,endyr-3)); - init_number n_pred_grp_nonpoll // the number of non-pollock predator groups - init_number n_pred_grp_poll // the number of pollock predator groups - number n_pred_grp // the total number of predator groups - !! n_pred_grp = n_pred_grp_nonpoll + n_pred_grp_poll; // the number of predator groups - init_matrix N_pred(1,n_pred_grp,styr,endyr) // the abindance of the predator groups - init_int nstrata_pred // the number of strata for computing spatial predation - init_vector strata(1,nstrata_pred) // vector of strata names - init_number n_pred_ages // the number of ages which are preyed upon - init_vector pred_ages(1,n_pred_ages) // the ages which are preyed upon - init_3darray poll_dist(1,n_pred_ages,styr,endyr,1,nstrata_pred) // the distribution of age 1 pollock across strata, by year - init_3darray pred_dist_nonpoll(1,n_pred_grp,styr,endyr,1,nstrata_pred) // the distribution of predators by strata, by year - 3darray Npred_bystrata_nonpoll(styr,endyr,1,n_pred_grp_nonpoll,1,nstrata_pred) // the number of non-pollock predators by strata for a given year - init_vector area_pred(1,nstrata_pred) // the strata area for getting predatuion by area (sq km) - // end test imput for ragged array - init_ivector nyrs_cons_nonpoll(1,n_pred_grp_nonpoll) // the number of years for which we have estimates of nonpollock predator consumption - init_matrix yrs_cons_nonpoll(1,n_pred_grp_nonpoll,1,nyrs_cons_nonpoll) // the years for which we have predator consumption estimates - init_matrix obs_cons_nonpoll(1,n_pred_grp_nonpoll,1,nyrs_cons_nonpoll) // the time series of consumption estiates for each predator (kilotons) - init_3darray oac_cons_nonpoll(1,n_pred_grp_nonpoll,1,nyrs_cons_nonpoll,1,n_pred_ages) // the observed age comps for predation for each predator, by year - init_matrix sam_oac_cons_nonpoll_raw(1,n_pred_grp_nonpoll,1,nyrs_cons_nonpoll) // the sample size for the consumption age comps - 3darray obs_cons_wgt_atage_nonpoll(1,n_pred_grp_nonpoll,1,nyrs_cons_nonpoll,1,n_pred_ages) // the observed weight consumed at age (by predator and year) - 3darray obs_cons_natage_nonpoll(1,n_pred_grp_nonpoll,1,nyrs_cons_nonpoll,1,n_pred_ages) // the observed consumed number at age (by predator and year) - 3darray obs_cpup_nonpoll(1,n_pred_grp_nonpoll,1,nyrs_cons_nonpoll,1,n_pred_ages) // the observed number consumed per unit predator - - - init_vector C_a_nonpoll(1,n_pred_grp_nonpoll) // C_a parameter for max consumption - init_vector C_b_nonpoll(1,n_pred_grp_nonpoll) // C_b parameter for max consumption - init_vector TCM_nonpoll(1,n_pred_grp_nonpoll) // TCM parameter for max consumption temperature function - init_vector TC0_nonpoll(1,n_pred_grp_nonpoll) // TC0 parameter for max consumption temperature funcion - init_vector CQ_nonpoll(1,n_pred_grp_nonpoll) // CQ parameter for max consumption temperature function - init_matrix temp_bystrata(styr,endyr,1,nstrata_pred) // the temperature by strata (year before 1982 are an average) - init_matrix mn_wgt_nonpoll(1,n_pred_grp_nonpoll,styr,endyr) // the mean weight of the predator groups - vector Y_nonpoll(1,n_pred_grp_nonpoll) // Y number for max consumption temperature function - vector Z_nonpoll(1,n_pred_grp_nonpoll) // Z number for max consumption temperature function - vector X_nonpoll(1,n_pred_grp_nonpoll) // X number for max consumption temperature function - 3darray V_nonpoll(styr,endyr,1,n_pred_grp_nonpoll,1,nstrata_pred) // V matrix for max consumption temperature function - 3darray F_t_nonpoll(styr,endyr,1,n_pred_grp_nonpoll,1,nstrata_pred) - 3darray Cmax_nonpoll(styr,endyr,1,n_pred_grp_nonpoll,1,nstrata_pred) // Cmax by year, predator group, and region - init_vector Cmax_avg(1,n_pred_grp); // Cmax for avg weight and temperature, for functional response - vector atf_wgts(1,n_pred_grp); // mean atf weights for computing functional response - vector poll_wgts(1,n_pred_ages); // mean poll weights for computing functional response - ivector comp_nr_ub(1,n_pred_grp_nonpoll); // define upper bound for ragged matrix cons_nr as integer vector - !! comp_nr_ub = ivector(nyrs_cons_nonpoll*n_pred_ages); - init_int phase_cope; - init_int n_cope; - init_ivector yrs_cope(1,n_cope); - init_vector obs_cope(1,n_cope); - init_vector obs_cope_std(1,n_cope); - vector lse_cope(1,n_cope); - vector lvar_cope(1,n_cope); - - - int k // added by Paul (k,m,yr_ind,z) - int m - int yr_ind - int z - - LOCAL_CALCS - lse_cope = elem_div(obs_cope_std,obs_cope); - lse_cope = sqrt(log(square(lse_cope) + 1.)); - lvar_cope = square(lse_cope); - write_log(SST); - write_log( n_pred_grp_nonpoll); - write_log( n_pred_grp_poll); - write_log( n_pred_grp); - write_log( N_pred); - write_log( nstrata_pred); - write_log( strata); - write_log( n_pred_ages); - write_log( pred_ages); - write_log( poll_dist); - write_log( pred_dist_nonpoll); - write_log( area_pred); - write_log( nyrs_cons_nonpoll); - write_log( yrs_cons_nonpoll); - write_log( obs_cons_nonpoll); - write_log(sam_oac_cons_nonpoll_raw); - write_log(temp_phase); - write_log(C_a_nonpoll); - write_log(C_b_nonpoll); - write_log(TCM_nonpoll); - write_log(TC0_nonpoll); - write_log(CQ_nonpoll); - write_log(temp_bystrata); - write_log(mn_wgt_nonpoll); - write_log(Cmax_avg); - write_log(phase_cope); - write_log(n_cope); - write_log(yrs_cope); - write_log(obs_cope); - write_log(obs_cope_std); - - // Assign values to temperature and predation phases (if estimating predation mortality or climate enhanced recruitment) - if(do_pred==1 && do_mult_func_resp==1) - { - do_pred_phase_ms = pred_phase; - do_pred_phase_ss = -1; - do_pred_phase = pred_phase; - } - else if (do_pred==1 && do_mult_func_resp!=1) - { - do_pred_phase_ms = -1; - do_pred_phase_ss = pred_phase; - do_pred_phase = pred_phase; - } - else - { - do_pred_phase_ms = -1; - do_pred_phase_ss = -1; - do_pred_phase = -1; - } - if(do_temp==1) - do_temp_phase = temp_phase; - else - do_temp_phase = -1; - // If estimating predation mortality, do a bunch of preliminary calculations - if (do_pred==1) - { - //Rescale spatial distribution matrices to add to one, and compute predator by year and strata - for (i=1;i<=n_pred_ages;i++) - { - for (ii=styr;ii<=endyr;ii++) - { - poll_dist(i,ii) = poll_dist(i,ii)/sum(poll_dist(i,ii)); - } - } - - for (i=1;i<=n_pred_grp;i++) - { - for (ii=styr;ii<=endyr;ii++) - { - pred_dist_nonpoll(i,ii) = pred_dist_nonpoll(i,ii)/sum(pred_dist_nonpoll(i,ii)); - Npred_bystrata_nonpoll(ii,i) = N_pred(i,ii)*pred_dist_nonpoll(i,ii); - } - } - - // Compute the catch per unit predator (CPUP) in numbers - for (m=1;m<=n_pred_grp;m++) - { - for (i=1;i<=nyrs_cons_nonpoll(m);i++) - { - yr_ind = yrs_cons_nonpoll(m,i) - 1981; // for getting the index for the wt_bts - if(yr_ind<1) - yr_ind = 1; - if(yr_ind>38) - yr_ind = 38; - obs_cons_wgt_atage_nonpoll(m,i) = oac_cons_nonpoll(m,i)*obs_cons_nonpoll(m,i); // kilotons - obs_cons_natage_nonpoll(m,i) = elem_div(obs_cons_wgt_atage_nonpoll(m,i),wt_bts(yr_ind)(1,n_pred_ages)); - obs_cpup_nonpoll(m,i) = obs_cons_natage_nonpoll(m,i)/N_pred(m,yrs_cons_nonpoll(m,i)); - } - } - - // Compute things for Cmax - for (i=1;i<=n_pred_grp_nonpoll;i++) - { - Y_nonpoll(i) = log(CQ_nonpoll(i))*(TCM_nonpoll(i)-TC0_nonpoll(i)+2); - Z_nonpoll(i) = log(CQ_nonpoll(i))*(TCM_nonpoll(i)-TC0_nonpoll(i)); - X_nonpoll(i) = Z_nonpoll(i)*Z_nonpoll(i)*pow(1+sqrt(1+40.0/Y_nonpoll(i)),2)/400.0; - for (j=styr;j<=endyr;j++) - { - V_nonpoll(j,i) = (TCM_nonpoll(i) - temp_bystrata(j))/(TCM_nonpoll(i) - TC0_nonpoll(i)); - F_t_nonpoll(j,i) = elem_prod(pow(V_nonpoll(j,i),X_nonpoll(i)),mfexp(X_nonpoll(i)*(1.0- V_nonpoll(j,i)))); - Cmax_nonpoll(j,i) = 365*C_a_nonpoll(i)*pow(mn_wgt_nonpoll(i,j),C_b_nonpoll(i))*F_t_nonpoll(j,i); - } - } - - // Compute average of atf and pollock weights for computing functional response - for (i=1;i<=n_pred_grp;i++) - atf_wgts(i) = mean(mn_wgt_nonpoll(i)); - for(i=1;i<=n_pred_ages;i++) - poll_wgts(i) = 1000.0*mean(column(wt_bts,i)); // convert from kg to g - } - write_log(do_temp_phase); // copy the temperture phase here - write_log(do_pred_phase_ss); // phase for estimating the predation parameters, single species function response - write_log(do_pred_phase_ms); // phase for estimating the predation parameters, multi-species function response - write_log(do_pred_phase); // phase for estimating the residual M - - END_CALCS - init_number test_2 - !! write_log(test_2); - !! if(test_2!=1234567){ cout<<"Failed on data read "<0) ph_log_fsh2 = phase_logist_fsh+1;else ph_log_fsh2 = phase_logist_fsh; - init_number sel_dif2_fsh(ph_log_fsh2) - - init_bounded_dev_vector sel_dif1_fsh_dev(styr,endyr_r,-5,5,phase_logist_fsh_dev) // allow for variability in survey selectivity inflection - init_bounded_dev_vector sel_a501_fsh_dev(styr,endyr_r,-5,5,phase_logist_fsh_dev) // allow for variability in survey selectivity inflection - init_bounded_dev_vector sel_trm2_fsh_dev(styr,endyr_r,-.5,.5,-2) // allow for variability in survey selectivity inflection - - // Parameters for computing SPR rates - number SPR_ABC // (0.05,2.,phase_F40) - sdreport_vector endyr_N(1,nages); - sdreport_number B_Bnofsh; - sdreport_number Bzero; - sdreport_number Percent_Bzero; - sdreport_number Percent_Bzero_1; - sdreport_number Percent_Bzero_2; - sdreport_number Percent_B100; - sdreport_number Percent_B100_1; - sdreport_number Percent_B100_2; - // Added to test for significance of "mean temperature" from earlier assessments (effect on survey q) - // NOT presented as a sensitivity in 2008 nor 2009 - // sdreport_vector q_temp(1,5); - vector q_temp(1,5); - sdreport_number SPR_OFL // (0.05,2.,phase_F40) - sdreport_number F40 // (0.05,2.,phase_F40) - sdreport_number F35 // (0.05,2.,phase_F40) - sdreport_vector SSB(styr,endyr_r) - - // Stuff for SPR and yield projections - number sigmarsq_out - number ftmp - number SB0 - number SBF40 - number SBF35 - number sprpen - number F_pen - number meanrec - vector SR_resids(styr_est,endyr_est); - matrix Nspr(1,4,1,nages) - sdreport_vector sel_fut(1,nages) - - 3darray natage_future(1,nscen,styr_fut,endyr_fut,1,nages) - init_vector rec_dev_future(styr_fut,endyr_fut,phase_F40); - - 3darray F_future(1,nscen,styr_fut,endyr_fut,1,nages); - matrix Z_future(styr_fut,endyr_fut,1,nages); - matrix S_future(styr_fut,endyr_fut,1,nages); - matrix catage_future(styr_fut,endyr_fut,1,nages); - number avg_rec_dev_future - - matrix eac_fsh(1,n_fsh_r,1,nbins) //--Expected proportion at age in Fishery - vector elc_fsh(1,nlbins) //--Expected proportion at length in Fishery - matrix eac_bts(1,n_bts_r,1,nbins) //--Expected proportion at age in trawl survey - matrix eac_cmb(1,n_bts_r,1,nbins) //--Expected proportion at age in combined surveys - matrix oac_cmb(1,n_bts_r,1,nbins) //--observed proportion at age in combined surveys - matrix eac_ats(1,n_ats_ac_r,1,nbins)//--Expected proportion at age in hydro survey - vector ea1_ats(1,n_ats_ac_r) //--Expected age 1 index from hydro survey - vector et_fsh(1,n_fsh_r) //--Expected total numbers in Fishery - vector et_bts(1,n_bts_r) //--Expected total numbers in Survey - vector et_cmb(1,n_bts_r) //--Expected total numbers in HydroSurvey - vector avail_bts(1,n_bts_r) //--Availability estimates in BTS - vector avail_ats(1,n_bts_r) //--Availability estimates in HydroSurvey - vector sigma_cmb(1,n_bts_r) //--Std errors of combined surveys (by year) - vector var_cmb(1,n_bts_r) //--Std errors of combined surveys (by year) - vector ot_cmb(1,n_bts_r) //--Observed total numbers in combined Surveys - vector eb_bts(1,n_bts_r) //--Expected total biomass in Survey - vector eb_ats(1,n_ats_r) //--Expected total biomass in Survey - vector et_ats(1,n_ats_r) //--Expected total numbers in HydroSurvey - vector lse_ats(1,n_ats_r) - vector lvar_ats(1,n_ats_r) - vector et_avo(1,n_avo_r) //--Expected total numbers in HydroSurvey - vector et_cpue(1,n_cpue) //--Expected total numbers in CPUE - vector Fmort(styr,endyr_r) - - matrix catage(styr,endyr_r,1,nages) - vector pred_catch(styr,endyr_r) - vector Pred_N_bts(styr,endyr_r) - vector Pred_N_ats(styr,endyr_r) - vector pred_cpue(1,n_cpue) - vector pred_cope(1,n_cope) - sdreport_vector Nage_3(styr,endyr_r) - vector pred_avo(1,n_avo) - // vector SSB(styr,endyr_r) - matrix natage(styr,endyr_r,1,nages); - vector srmod_rec(styr_est,endyr_est); - matrix Z(styr,endyr_r,1,nages); - matrix F(styr,endyr_r,1,nages); - matrix S(styr,endyr_r,1,nages); - matrix M(styr,endyr_r,1,nages); - init_bounded_vector M_dev(styr+1,endyr_r,-1.,1.,-8); - matrix log_sel_fsh(styr,endyr_r,1,nages); - matrix sel_fsh(styr,endyr_r,1,nages); - matrix log_sel_bts(styr,endyr_r,1,nages); - matrix log_sel_ats(styr,endyr_r,1,nages); - number ff; - number catch_like; - number avgsel_fsh; - number avgsel_bts; - number avgsel_ats; - number bzero; - number surv; - number nthisage; - vector surv_like(1,3); - number cpue_like; - number cope_like; - number avo_like; - vector sel_like(1,3); - vector sel_like_dev(1,3); - vector age_like(1,ngears); - matrix age_like_yr(1,ngears,1,nagecomp) - number len_like; - number wt_like; - vector age_like_offset(1,ngears); - number len_like_offset; - number MN_const // Multinomial constant - !! MN_const = 1e-3; // Multinomial constant - vector Priors(1,4); - vector rec_like(1,7); - vector all_like(1,26); - number sumtmp; - number tmpsp; - vector log_initage(2,nages); - vector pred_biom(styr,endyr_r); - vector SRR_SSB(1,20) - vector fake_SST(1,40) // added by Paul - vector fake_dens(1,40) // added by Paul - - // Average weight stuff - init_bounded_number L1(10,50,2); - init_bounded_number L2(30,90,3); - init_number log_alpha(-1); - init_number log_K(4); - vector wt_inc(age_st,age_end-1); - init_matrix d_scale(1,nscale_parm,age_st,age_end,phase_d_scale); - number K; - // init_number log_t0(3); - // Predicted weight matrix - 3darray wt_hat(1,ndat_wt,1,max_nyrs_data,age_st,age_end); - number alphawt; - matrix wt_pre(styr_wt,endyr_wt,age_st,age_end) - vector mnwt(age_st,age_end); - init_bounded_vector coh_eff(styr_wt-nages_wt-age_st+1,endyr_wt-age_st+3,-15,15,phase_coheff); - init_bounded_vector yr_eff(styr_wt,endyr_wt+3,-15,15,phase_yreff); - - sdreport_vector wt_last(age_st,age_end); - sdreport_vector wt_cur(age_st,age_end); - sdreport_vector wt_next(age_st,age_end); - sdreport_vector wt_yraf(age_st,age_end); - - sdreport_number avg_age_msy; - // sdreport_number avgln_msy; - sdreport_number avgwt_msy; - sdreport_number MSY; - sdreport_number MSY_wt; - sdreport_number Fmsy; - sdreport_number Fmsy_wt; // includes vector on value...default is equal weight - sdreport_number Fmsy2; - sdreport_number Fmsy2_wt; - sdreport_vector Fmsy2_dec(1,10); // uncertain next year weight previous decade selectivity estimates - sdreport_vector Fmsy2_decwt(1,10); // uncertain next year selectivity but previous decade weight estimates - sdreport_number Bmsy2; - sdreport_number Bmsy2_wt; - sdreport_number lnFmsy; - sdreport_number lnFmsy2; - sdreport_number SER_Fmsy; - sdreport_number Fendyr_Fmsy; - sdreport_number Rmsy; - sdreport_number Bmsy; - - // Decision table variables............................... - sdreport_vector Fcur_Fmsy(1,nscen); - sdreport_vector Bcur_Bmsy(1,nscen); - sdreport_vector Bcur_Bmean(1,nscen); - sdreport_vector Bcur3_Bcur(1,nscen); - sdreport_vector Bcur3_Bmean(1,nscen); - sdreport_vector Bcur2_Bmsy(1,nscen); - sdreport_vector Bcur2_B20(1,nscen); - sdreport_vector LTA1_5R(1,nscen); // long term average age 1_5 Ratio - sdreport_vector LTA1_5(1,nscen); // long term average age 1_5 - sdreport_vector MatAgeDiv1(1,nscen); // Diversity of Age structure in mature population - sdreport_vector MatAgeDiv2(1,nscen); // Diversity of Age structure in mature population - vector H(styr,endyr_r); - vector avg_age_mature(styr,endyr_r); - sdreport_vector RelEffort(1,nscen); // Effort relative to endyr - - sdreport_number F40_spb; - // sdreport_number F40_catch; - sdreport_number begbiom; - sdreport_number DepletionSpawners; - sdreport_number SB100; - sdreport_number Current_Spawners; - sdreport_vector pred_rec(styr,endyr_r); - sdreport_vector age_3_plus_biom(styr,endyr_r+2); - sdreport_vector ABC_biom(1,10); - sdreport_vector ABC_biom2(1,10); - sdreport_vector rechat(1,20); - vector SRresidhat(1,40); // added by Paul, the predicted residual - sdreport_vector SER(styr,endyr_r); - matrix future_SER(1,nscen,styr_fut,endyr_fut); - matrix future_catch(1,nscen,styr_fut,endyr_fut); - sdreport_matrix future_SSB(1,nscen,endyr_r,endyr_fut) - sdreport_vector Age3_Abund(styr,endyr_r) - vector age_1_7_biomass(styr,endyr_r); - vector NLL(1,20) - objective_function_value fff; - - // Things for estimating the relation between temperature and SR residuals - init_number resid_temp_x1(do_temp_phase); - init_number resid_temp_x2(do_temp_phase); - vector SR_resids_temp(styr_est,endyr_est); - number SR_resids_like; - - // Things for estimating the predation mortality rates (added by Paul) - init_bounded_matrix log_a_II(1,n_pred_grp,1,n_pred_ages,-12,0,do_pred_phase_ss) // the "a" parameter for Holling type 2(log scale, by predator and age) - init_bounded_matrix log_b_II(1,n_pred_grp,1,n_pred_ages,0,12,do_pred_phase_ss) // the "b" parameter for Holling type 2(log scale, by predator and age) - init_bounded_vector log_a_II_vec(1,n_pred_grp,-12,0,do_pred_phase_ms) // the "a" parameter for Holling type 2(log scale, by predator) - init_bounded_vector log_b_II_vec(1,n_pred_grp,0,12,do_pred_phase_ms) // the "b" parameter for Holling type 2(log scale, by predator) - init_bounded_matrix log_rho(1,n_pred_grp,1,n_pred_ages,-12,0,do_pred_phase_ms) // the proportion of proportion for MS Holling - //matrix log_rho2(1,n_pred_grp,1,n_pred_ages); // after rescaling - - init_bounded_vector log_resid_M(1,n_pred_ages,-3,0.1,do_pred_phase) // the residual natural mortality for the ages that are preyed upon (on log scale) - vector resid_M(1,n_pred_ages) - vector resid_M_like(1,n_pred_ages) // the likelihood for fitting the residual natural mortality rates - matrix a_II(1,n_pred_grp,1,n_pred_ages) // unlogged versions of the Holling type II parameters - matrix b_II(1,n_pred_grp,1,n_pred_ages) - matrix rho(1,n_pred_grp,1,n_pred_ages) // the proportion of proportion for MS Holling - vector a_II_vec(1,n_pred_grp) // unlogged versions of the Holling type II parameters - vector b_II_vec(1,n_pred_grp) - - - 3darray natage_strat(styr,endyr_r,1,n_pred_ages,1,nstrata_pred) // the number of poll by strata for a given year and age - 3darray natage_strat_dens(styr,endyr_r,1,n_pred_ages,1,nstrata_pred) // the density of poll by strata for a given year and age - - matrix meannatage(styr,endyr_r,1,nages) // the mean N summed across the strata - 3darray meannatage_bystrata(styr,endyr_r,1,nages,1,nstrata_pred) // the mean N by strata - 3darray mean_dens_bystrata(styr,endyr_r,1,nages,1,nstrata_pred) // mean density at age by strata - matrix mean_dens(styr,endyr_r,1,n_pred_ages) // the mean density over all the area, within year and strata - 4darray cons(1,n_pred_grp,styr,endyr_r,1,n_pred_ages,1,nstrata_pred); // the consumption of pollock, by predator, year, pollock age, and area - 3darray natmort_pred(1,n_pred_ages,styr,endyr_r,1,nstrata_pred); // the natural mortality with spatial predation - 4darray M_pred(styr,endyr_r,1,n_pred_ages,1,nstrata_pred,1,n_pred_grp_nonpoll); // the inst predation rate for pollock, of pollock, by predator, year, pollock age, and area - matrix M_pred_tmp(1,n_pred_grp_nonpoll,1,n_pred_ages) // temporary place for output of function get_Mpred2 - 3darray M_pred_sum(styr,endyr_r,1,n_pred_ages,1,nstrata_pred) // the total M across predators within a year, pollock age, and strata - 3darray Z_pred(styr,endyr_r,1,n_pred_ages,1,nstrata_pred); // the Z for a1 pollock, by year and area - 3darray S_pred(styr,endyr_r,1,n_pred_ages,1,nstrata_pred); // the S for a1 pollock, by year and area - matrix M_pred_avg(1,n_pred_ages,styr,endyr_r); // average of predation mortality across the strata, weighted by poll abundance at start of year - 3darray cons_atage(1,n_pred_grp,styr,endyr_r,1,n_pred_ages); // total consumption by predator, year and age (sums over strata) - 3darray cons_atage_wt(1,n_pred_grp,styr,endyr_r,1,n_pred_ages); // total weight of consumption by predator, year and age (sums over strata) - matrix pred_cons(1,n_pred_grp,styr,endyr_r); // the predicted total consumption bny predator within a year (across ages) - 3darray eac_cons(1,n_pred_grp,styr,endyr_r,1,n_pred_ages); // the estimated age comp for consumption(by predator and year, sums over strata) - vector ssq_cons(1,n_pred_grp); // the sum of squares for total consumption - vector oac_cons_like_offset(1,n_pred_grp); // offset for multinomial - vector age_like_cons(1,n_pred_grp); // likelihood for the age comps - 3darray pred_cpup(1,n_pred_grp,styr,endyr_r,1,n_pred_ages) // the predicted CPUP by predator, year, and age - 4darray implied_cpuppa(1,n_pred_grp,styr,endyr_r,1,n_pred_ages,1,nstrata_pred); // predicted catch per unit predator per area - 4darray implied_obs_cons_bystrata(1,n_pred_grp,styr,endyr_r,1,n_pred_ages,1,nstrata_pred) // implied observed consumption by strata - 4darray implied_prop_Cmax(1,n_pred_grp,styr,endyr_r,1,n_pred_ages,1,nstrata_pred) // implied consumption rate as proportion of Cmax - - // things added by Paul for the yield curve - number max_F_yldcrv; // maximum F for yield curve, based on spr calc - vector F_yldcrv(1,40); // vector of values for yield curve - vector yield_curve(1,40); // equilibrium yield curve - //vector yield_curve_fut(1,40); // Eq. yield curve with future sst - vector natmort_fut(1,nages); // for getting yield curve and Fmsy; reflects recent predation morality - - // updated compweights -- added by Paul, uses McAllister-Ianelli method - vector compweightsnew(1,n_pred_grp_nonpoll) - matrix comp_mcian_wgt_inv(1,n_pred_grp_nonpoll,1,nyrs_cons_nonpoll) // inverse of McAllister-Ianelli weights, consumption estimates - matrix comp_mcian_wgt(1,n_pred_grp_nonpoll,1,nyrs_cons_nonpoll) - - // consumption estimates normalized resids -- added by Paul to use sdnr to iterate the consumtion weights - vector consweightsnew(1,n_pred_grp_nonpoll) - matrix cons_nr(1,n_pred_grp_nonpoll,1,nyrs_cons_nonpoll) - matrix comp_nr(1,n_pred_grp_nonpoll,1,comp_nr_ub) // for holding the composition standardized resids - - // alpha terms that scale the environment effect to log_avg_rec and modeled recruitment (following Maunder and Watters 2002 paper) - number pred_rec_alpha - number srmod_rec_alpha - - - - - -PRELIMINARY_CALCS_SECTION - // fixed_catch_fut1 = fixed_catch_fut1 + 0.1 ; - wt_fut = wt_fsh(endyr_r); // initializes estimates to correct values...Eq. 21 - // base_natmort(1)=.9; base_natmort(2)=.45; for (j=3 ;j<=nages;j++) base_natmort(j)=natmortprior; - base_natmort = natmort_in; - - natmort = base_natmort; - // cout <<"M input= "< 0.4 ) ignore_last_ats_age1 = 1; else ignore_last_ats_age1=0; - // cout <<" Last age comp in BTS: " << endl << oac_bts_data(n_bts) << endl; - // cout <<" Age_like_offset: " << endl << age_like_offset << endl; - Cat_Fut(1) = next_yrs_catch; // catch guess - // Simple decrement of future cathes to capture relationship between adjustments (below Bmsy) w/in same year - for (i=2;i<=10;i++) - Cat_Fut(i) = Cat_Fut(i-1)*.90; - write_log(Cat_Fut); - - // cout << "Next year's catch and decrements"<0)) - Profile_F(); - if (mceval_phase()) - write_eval(); - update_compweights(); // added by Paul - -FUNCTION update_compweights // added by Paul, update if the comp is estmated, otherwise carry over previous comp weight - if(active(log_resid_M)) - { - for (i=1;i<=n_pred_grp_nonpoll;i++) - { - compweightsnew(i) = 1.0/mean(comp_mcian_wgt_inv(i)); // weights for consumption age comps - consweightsnew(i) = std_dev(cons_nr(i)); // weights for consumption estimates - } - } - else - { - compweightsnew = compweights; - consweightsnew = consweights; - } - -FUNCTION Get_Selectivity - avgsel_fsh.initialize(); - avgsel_bts.initialize(); - avgsel_ats.initialize(); - //cout<<"InSel"<0) // For logistic fishery selectivity (not default, NOT USED) - { - dvariable dif; - dvariable inf; - if (active(sel_a501_fsh_dev)) - { - dvar_matrix seldevs_tmp(1,3,styr,endyr_r); - seldevs_tmp(1) = sel_dif1_fsh_dev; - seldevs_tmp(2) = sel_a501_fsh_dev; - seldevs_tmp(3) = sel_trm2_fsh_dev; - log_sel_fsh = compute_selectivity1(styr,sel_dif1_fsh,sel_a501_fsh,sel_trm2_fsh,seldevs_tmp); - } - else - log_sel_fsh = compute_selectivity1(styr,sel_dif1_fsh,sel_a501_fsh,sel_trm2_fsh); - } - - //cout<<"InSel"<38) - yr_ind = 38; - - for (k=1;k<=n_pred_ages;k++) // loop over ages of pollock that are preyed upon - { - natage_strat(i,k) = natage(i,pred_ages(k))*poll_dist(k,i); // distribute the age k pollock in each area - natage_strat_dens(i,k) = elem_div(natage_strat(i,k),area_pred); // the density of age k pollock in each area - } - - // do the multispecies funcction response thing - if (do_mult_func_resp==1) - { - for (j=1;j<=nstrata_pred;j++) // loop over strata, and get the mean abundance accounting for all mortality - { - M_pred_tmp = get_Mpred2(column(natage_strat_dens(i),j),resid_M,F(i)(1,n_pred_ages), - column(Npred_bystrata_nonpoll(i),j), - column(mn_wgt_nonpoll,i), - wt_bts(yr_ind)(1,n_pred_ages)*1000, - column(Cmax_nonpoll(i),j),rho, - a_II_vec,b_II_vec, j); - - // cout << "M_pred_tmp is " << M_pred_tmp << endl; - // loop to get the results in the right structure - for (ii=1;ii<=n_pred_ages;ii++) - { - for (jj=1;jj<=n_pred_grp;jj++) - { - M_pred(i,ii,j,jj) = M_pred_tmp(jj,ii); - //cout <<" prey age is " << ii<< " pred_grp is "<< jj << " strata is " << j << " M_pred is " << M_pred(i,ii,j,jj) << endl; - } - } - } - } // loop if doing the multispecies functional response - - // loop over ages of pollock that are preyed upon - for (k=1;k<=n_pred_ages;k++) - { - // loop over strata, and get the mean abundance accounting for all mortality - for (j=1;j<=nstrata_pred;j++) - { - if (do_mult_func_resp != 1) // loop is doing the single species functional response - { - // natage_strat(styr,endyr_r,1,n_pred_ages,1,nstrata_pred) // the number of poll by strata for a given year and age - // Jim thinks this may be easier implemented by passing the indices (i,j,k)... - M_pred(i,k,j) = get_Mpred(natage_strat_dens(i,k,j), - F(i,pred_ages(k)), - resid_M(k), - column(natage_strat_dens(i),j), - column(a_II,k), - column(b_II,k), - column(Npred_bystrata_nonpoll(i),j), - column(mn_wgt_nonpoll,i)/(wt_bts(yr_ind,k)*1000), - column(Cmax_nonpoll(i),j),rho); - // e.g.,: M_pred(i,k,j) = get_Mpred(i,k,j); - } - M_pred_sum(i,k,j) = sum(M_pred(i,k,j)); // sum across the different predators - natmort_pred(k,i,j) = M_pred_sum(i,k,j) + resid_M(k); - Z_pred(i,k,j) = F(i,pred_ages(k)) + natmort_pred(k,i,j); // get the total Z by year, pollock age, and strata - S_pred(i,k,j) = mfexp(-1.0*Z_pred(i,k,j)); - } // close strata loop - // cout << "M_pred is " << M_pred<< endl; - meannatage_bystrata(i,k) = elem_div(elem_prod(natage_strat(i,k),(1.-S_pred(i,k))),Z_pred(i,k)); // the mean number at age by strata - meannatage(i,k) = sum(meannatage_bystrata(i,k)); // the mean N summed across the strata - mean_dens_bystrata(i,k) = elem_div(meannatage_bystrata(i,k),area_pred); - mean_dens(i,k) = (meannatage(i,k))/sum(area_pred); // the mean density summed over the strata within a year and age - - for (j=1;j<=nstrata_pred;j++) - { // loop over strata, and get the mean abundance accounting for all mortality - for (m=1;m<=n_pred_grp;m++) - { // loop over predators, and get the consumption and M by predator, pollock age, year, and strata - cons(m,i,k,j) = M_pred(i,k,j,m)*meannatage_bystrata(i,k,j); - } - } - M_pred_avg(k,i) = (M_pred_sum(i,k)*meannatage_bystrata(i,k))/sum(meannatage_bystrata(i,k)); // get an average M_pred across the strata, weighted by the beginning abundance in each strata - // get the total survival by year, pollock age, and strata - if(i < endyr_r) - natage(i+1,k+1) = natage_strat(i,k)*S_pred(i,k); // the natage after predation (summed over strata) - SSB(i) += natage_strat(i,k)*pow(S_pred(i,k),yrfrac)*p_mature(pred_ages(k))*wt_ssb(i,pred_ages(k)); - } // end age loop - SSB(i) += elem_prod(elem_prod(natage(i)(n_pred_ages+1,nages), - pow(S(i)(n_pred_ages+1,nages),yrfrac)),p_mature(n_pred_ages+1,nages))*wt_ssb(i)(n_pred_ages+1,nages); // Eq. 1 - - meannatage(i)(n_pred_ages+1,nages) = elem_prod(elem_div(1.-S(i)(n_pred_ages+1,nages),Z(i)(n_pred_ages+1,nages)), - natage(i)(n_pred_ages+1,nages)); // the mean n at age for the ages not preyed upon - if (i < endyr_r) - { - natage(i+1)(n_pred_ages+2,nages) = ++elem_prod(natage(i)(n_pred_ages+1,nages-1), S(i)(n_pred_ages+1,nages-1)); - natage(i+1,nages) += natage(i,nages)*S(i,nages); // Eq. 1 - } - } // end year loop - // added by Paul to get M values the reflect the recent predation - for (i=1;i<=n_pred_ages;i++) - { - natmort_fut(i) = mean(M_pred_avg(i)(endyr_r-4,endyr_r) + resid_M(i)); - } - } // ********* end of thing by Paul ******************* - else - { - SSB(styr) = elem_prod(elem_prod(natage(styr),pow(S(styr),yrfrac)),p_mature)*wt_ssb(styr); // Eq. 1 - natage(styr+1)(2,nages) = ++elem_prod(natage(styr)(1,nages-1), S(styr)(1,nages-1)); // Eq. 1 - natage(styr+1,nages) += natage(styr,nages)*S(styr,nages); // Eq. 1 - for (i=styr+1;i0) - { - tmp_abun_end = tmp_abun_bg*mfexp(-tmp_F -tmp_M); - avg_N = tmp_abun_bg*(1-mfexp(-tmp_F -tmp_M))/(tmp_F + tmp_M); - M_pred = elem_prod(elem_prod(wt_ratio,elem_prod(tmp_Cmax,elem_prod(a,Npred))),(1/(b+avg_N))); - M_pred_sum = sum(M_pred); - //M_pred = elem_prod(a,Npred)*(1/(b+avg_N)); - cout << "the weight ratio is " << endl; - cout << wt_ratio << endl; - cout << "the tmp_Cmax is " << endl; - cout << tmp_Cmax << endl; - cout << "the a_II " << endl; - cout << a << endl; - cout << "the Npred " << endl; - cout << Npred << endl; - cout << "The b_II " << endl; - cout << b << endl; - cout << "The avg_N " << endl; - cout << avg_N << endl; - cout << "the M_pred " << endl; - cout << elem_prod(wt_ratio,elem_prod(tmp_Cmax,elem_prod(a,Npred)))*(1/(b+avg_N)) << endl; - - // Check this should be double not dvariable... - double dd = 10.; - int iter = 0; - // Check differentiability here... - while (iter < 10) - { - iter++; - prev_tmp_abun_end = tmp_abun_end; - tmp_abun_end = tmp_abun_bg*mfexp(-tmp_F -tmp_M -M_pred_sum); - avg_N = tmp_abun_bg*(1-mfexp(-tmp_F -tmp_M -M_pred_sum))/(tmp_F + tmp_M + M_pred_sum); - M_pred = elem_prod(elem_prod(wt_ratio,elem_prod(tmp_Cmax,elem_prod(a,Npred))),(1/(b+avg_N))); - M_pred_sum = sum(M_pred); - //M_pred = elem_prod(a,Npred)*(1/(b+avg_N)); - - dd = value(prev_tmp_abun_end) / value(tmp_abun_end) - 1.; - if (dd<0.) dd *= -1.; - // if(active(log_a_II)){ - // cout <<"in loop "<0) - { - tmp_abun_end = tmp_abun_bg*mfexp(-tmp_F -tmp_M); - avg_N = tmp_abun_bg*(1-mfexp(-tmp_F -tmp_M))/(tmp_F + tmp_M); - M_pred = elem_prod(elem_prod(wt_ratio,elem_prod(tmp_Cmax,elem_prod(a,Npred))),(1/(b+avg_N))); - M_pred_sum = sum(M_pred); - //M_pred = elem_prod(a,Npred)*(1/(b+avg_N)); - /*cout << "the weight ratio is " << endl; - cout << wt_ratio << endl; - cout << "the tmp_Cmax is " << endl; - cout << tmp_Cmax << endl; - cout << "the a_II " << endl; - cout << a << endl; - cout << "the Npred " << endl; - cout << Npred << endl; - cout << "The b_II " << endl; - cout << b << endl; - cout << "The avg_N " << endl; - cout << avg_N << endl; - cout << "the M_pred " << endl; - cout << elem_prod(wt_ratio,elem_prod(tmp_Cmax,elem_prod(a,Npred)))*(1/(b+avg_N)) << endl; - */ - - // Check this should be double not dvariable... - dvariable dd = 10.; - int iter = 0; - // Check differentiability here... - while (dd > 1e-6) - { - iter++; - prev_tmp_abun_end = tmp_abun_end; - tmp_abun_end = tmp_abun_bg*mfexp(-tmp_F -tmp_M -M_pred_sum); - avg_N = tmp_abun_bg*(1-mfexp(-tmp_F -tmp_M -M_pred_sum))/(tmp_F + tmp_M + M_pred_sum); - M_pred = elem_prod(elem_prod(wt_ratio,elem_prod(tmp_Cmax,elem_prod(a,Npred))),(1/(b+avg_N))); - M_pred_sum = sum(M_pred); - dd = prev_tmp_abun_end / tmp_abun_end - 1.; - //M_pred = elem_prod(a,Npred)*(1/(b+avg_N)); - - if (dd<0.) - dd *= -1.; - // if(active(log_a_II)){ - // cout <<"in loop "<0.0) - { - tmp_abun_end_vec = elem_prod(tmp_abun_vec,mfexp(-F_vec -tmp_M_vec)); - avg_N_vec = elem_prod(tmp_abun_vec,elem_div((1-mfexp(-F_vec -tmp_M_vec)),(F_vec + tmp_M_vec))); - - for (ii=1;ii<=n_pred_ages;ii++) // loop over prey ages - { - for (jj=1;jj<=n_pred_grp;jj++) // loop over number of predators - { - M_pred_mat(jj,ii) = (tmp_Cmax(jj)* (wts_pred(jj)/wts_prey(ii)) *a_vec(jj)*Npred(jj)*rho(jj,ii))/ - (b_vec(jj) + rho(jj)*avg_N_vec); - } - M_pred_sum_vec(ii) = sum(column(M_pred_mat,ii)); - } - - dvector dd_vec(1,n_pred_ages); - double dd_vec_sum = 10; - int iter = 0; - - while (dd_vec_sum > 1e-6) - { - iter++; - prev_tmp_abun_end_vec = tmp_abun_end_vec; - tmp_abun_end_vec = elem_prod(tmp_abun_vec,mfexp(-F_vec -tmp_M_vec -M_pred_sum_vec)); - avg_N_vec = elem_prod(tmp_abun_vec,elem_div((1-mfexp(-F_vec -tmp_M_vec -M_pred_sum_vec)),(F_vec + tmp_M_vec +M_pred_sum_vec))); - - for (ii=1;ii<=n_pred_ages;ii++) // loop over prey ages - { - for (jj=1;jj<=n_pred_grp;jj++) // loop over number of predators - { - M_pred_mat(jj,ii) = (tmp_Cmax(jj)*(wts_pred(jj)/wts_prey(ii))*a_vec(jj)*Npred(jj)*rho(jj,ii))/(b_vec(jj) + rho(jj)*avg_N_vec); - } - M_pred_sum_vec(ii) = sum(column(M_pred_mat,ii)); - } - - for (ii=1;ii<=n_pred_ages;ii++) - { - if (tmp_abun_end_vec(ii)>0.0) - { - dd_vec(ii) = value(prev_tmp_abun_end_vec(ii)) / value(tmp_abun_end_vec(ii)) - 1.; - if (dd_vec(ii)<0.) dd_vec(ii) *= -1.; - } - else dd_vec(ii) = 0.0; - } - dd_vec_sum = sum(dd_vec); - } - } - else - { - M_pred_mat = 0.0; - } - for (ii=1;ii<=n_pred_ages;ii++) - { - if (tmp_abun_end_vec(ii)==0.0) - { - for (jj=1;jj<=n_pred_grp;jj++) M_pred_mat(jj,ii) = 0.0; - } - } - RETURN_ARRAYS_DECREMENT(); - return M_pred_mat; - - -FUNCTION GetDependentVar - // For making some output for spiffy output - // Spiffy SR output - tmpsp=1.1*max(SSB); - if (tmpsp=2020) - { - regime(2) = mean(pred_rec(1978,2020)); - regime(5) = mean(pred_rec(1990,2020)); - regime(7) = mean(pred_rec(2000,2020)); - regime(8) = mean(pred_rec(1964,2020)); - regime(6) = mean(pred_rec(1990,1999)); - regime(3) = mean(pred_rec(1978,1999)); - } - else - { - if (endyr_r<=2000) - regime(7) = pred_rec(endyr_r); - else - regime(7) = mean(pred_rec(2000,endyr_r)); - - regime(2) = mean(pred_rec(1978,endyr_r)); - regime(5) = mean(pred_rec(1990,endyr_r)); - regime(8) = mean(pred_rec(1964,endyr_r)); - if (endyr_r<=1999) - { - regime(3) = mean(pred_rec(1978,endyr_r)); - regime(6) = mean(pred_rec(1990,endyr_r)); - } - else - { - regime(6) = mean(pred_rec(1990,1999)); - regime(3) = mean(pred_rec(1978,1999)); - } - } - if (!mceval_phase()) get_msy(); - - dvar_vector res(1,4); - res.initialize(); - res = get_msy_wt(); - Fmsy_wt = res(1); - MSY_wt = res(2); - Bmsy2_wt = res(3); - Fmsy2_wt = res(4); - // cout <<"Msy stuff: "<=1;iyr--) - { - res.initialize(); - sel_fut = sel_fsh(endyr_r-iyr+1); - // sel_fut /=sel_fut(6); // NORMALIZE TO AGE 6 - sel_fut /= mean(sel_fut); // NORMALIZE TO mean - if (!mceval_phase()) res = get_msy_wt(); - Fmsy2_dec(iyr) = res(4); - // cout <=1;iyr--) - { - res.initialize(); - wt_fut(3,nages) = wt_pre(endyr_r-iyr+1); - if (!mceval_phase()) res = get_msy_wt(); - Fmsy2_decwt(iyr) = res(4); - // cout <8) - ftmp= mean(F(endyr_r)) * ((double(k-1)-1.)*.05 + 0.5); // this takes endyr F and brackets it...for mean - else - ftmp = SolveF2(natage_future(k,styr_fut),dec_tab_catch(k)); - } - for (i=styr_fut;i<=endyr_fut;i++) - { - F_future(k,i) = sel_fut*ftmp; - Z_future(i) = F_future(k,i) + natmort; - S_future(i) = mfexp(-Z_future(i)); - dvariable criterion; - dvariable Bref ; - future_SSB(k,i) = elem_prod(elem_prod(natage_future(k,i),pow(S_future(i),yrfrac)), p_mature) * wt_ssb(endyr_r); - // Now compute the time of spawning SSB for everything else.... - // future_SSB(k,i) = elem_prod(elem_prod(natage_future(k,i),pow(S_future(i),yrfrac)), p_mature) * wt_ssb(endyr_r); - - if (phase_sr<0) //No Stock-recruitment curve for future projections-------- - natage_future(k,i, 1) = mfexp(log_avgrec + rec_dev_future(i)); - else //Use Stock-recruitment curve --------- - { - Xspawn =future_SSB(k,i-1); - natage_future(k,i,1) = SRecruit(Xspawn) * mfexp(rec_dev_future(i) ); - } - - if (i0 ) - { - for (i=endyr_r-(nyrs_sel_avg+1);i<=endyr_r;i++) - sel_fut = sel_fut + sel_fsh(i); - sel_fut/=nyrs_sel_avg; - } - else - sel_fut = sel_fsh(endyr_r+nyrs_sel_avg); // negative nyrs_sel_avg can be used to pick years for evaluation - - //sel_fut/=sel_fut(6); // NORMALIZE TO AGE 6 - sel_fut/=mean(sel_fut); // NORMALIZE TO mean - -FUNCTION compute_spr_rates - //Compute SPR Rates - F35 = get_spr_rates(.35); - F40 = get_spr_rates(.40); - -FUNCTION dvariable get_spr_rates(double spr_percent,dvar_vector sel) - RETURN_ARRAYS_INCREMENT(); - double df=1.e-3; - dvariable F1 ; - F1.initialize(); - F1 = .9*natmort(9); - dvariable F2; - dvariable F3; - dvariable yld1; - dvariable yld2; - dvariable yld3; - dvariable dyld; - dvariable dyldp; - // Newton Raphson stuff to go here - for (int ii=1;ii<=6;ii++) - { - F2 = F1 + df; - F3 = F1 - df; - yld1 = -100.*square(log(spr_percent/spr_ratio(F1,sel))); - yld2 = -100.*square(log(spr_percent/spr_ratio(F2,sel))); - yld3 = -100.*square(log(spr_percent/spr_ratio(F3,sel))); - dyld = (yld2 - yld3)/(2*df); // First derivative (to find the root of this) - dyldp = (yld3-(2*yld1)+yld2)/(df*df); // Newton-Raphson approximation for second derivitive - F1 -= dyld/dyldp; - } - RETURN_ARRAYS_DECREMENT(); - return(F1); -FUNCTION dvariable spr_ratio(dvariable trial_F,dvar_vector& sel) - RETURN_ARRAYS_INCREMENT(); - dvariable SBtmp; - dvar_vector Ntmp(1,nages); - dvar_vector srvtmp(1,nages); - SBtmp.initialize(); - Ntmp.initialize(); - srvtmp.initialize(); - dvar_vector Ftmp(1,nages); - Ftmp = sel*trial_F; - srvtmp = exp(-(Ftmp + natmort) ); - dvar_vector wttmp = wt_ssb(endyr_r); - Ntmp(1)=1.; - j=1; - SBtmp += Ntmp(j)*p_mature(j)*wttmp(j)*pow(srvtmp(j),yrfrac); - for (j=2;j3) get_combined_index(); - -FUNCTION Get_Cons_at_Age - // added by Paul - // if doing the spatial predation thing, get the consumption at age - for (i=styr;i<=endyr_r;i++) - { - for (k=1;k<=n_pred_ages;k++) - { - for (m=1;m<=n_pred_grp;m++) - { - cons_atage(m,i,k) = sum(cons(m,i,k)); // get ths consumption by predator, year, and age (summed over strata) - pred_cpup(m,i,k) = cons_atage(m,i,k)/N_pred(m,i); - } // predator loop - } // age loop - } // year loop - - for (i=styr;i<=endyr_r;i++) - { - yr_ind = i-1981; // for getting the index for the wt_bts - if(yr_ind<1) - yr_ind = 1; - if(yr_ind>38) - yr_ind = 38; - for (m=1;m<=n_pred_grp;m++) - { - cons_atage_wt(m,i) = elem_prod(cons_atage(m,i),wt_bts(yr_ind)(1,n_pred_ages)); - pred_cons(m,i) = sum(cons_atage_wt(m,i)); // the predicted consumption by predator and year (summed over age and strata) - eac_cons(m,i) = cons_atage_wt(m,i)/pred_cons(m,i); // the predicted consumption age comp (by predator and year) - } // predator loop - } // year loop - - for (i=styr;i<=endyr_r;i++) - { - yr_ind = i-1981; // for getting the index for the wt_bts - if(yr_ind<1) - yr_ind = 1; - if(yr_ind>38) - yr_ind = 38; - - for (m=1;m<=n_pred_grp;m++) - { - for (k=1;k<=n_pred_ages;k++) - { - // implied_obs_cons_bystrata(m,i,k) = obs_cons_natage_nonpoll(m,i,k)*(cons(m,iyr,k)/cons_atage(m,iyr,k)); - implied_obs_cons_bystrata(m,i,k) = cons(m,i,k); - for (z=1;z<=nstrata_pred;z++) - { - if (Npred_bystrata_nonpoll(i,m,z)>0) - { - implied_cpuppa(m,i,k,z) = implied_obs_cons_bystrata(m,i,k,z)/(Npred_bystrata_nonpoll(i,m,z)*area_pred(z)); - implied_prop_Cmax(m,i,k,z) = (implied_cpuppa(m,i,k,z)*area_pred(z))/ - ((mn_wgt_nonpoll(m,i)/(wt_bts(yr_ind,k)*1000))*Cmax_nonpoll(i,m,z)); - - // old code here -- used in model runs for Brazil paper - // implied_cpuppa(m,i,k,z) = - // (implied_obs_cons_bystrata(m,i,k,z)*((wt_bts(yr_ind,k)*1000)/mn_wgt_nonpoll(m,iyr)))/ - // (Cmax_nonpoll(iyr,m,z)*Npred_bystrata_nonpoll(iyr,m,z)*area_pred(z)); - // implied_prop_Cmax(m,i,k,z) = implied_cpuppa(m,i,k,z)*area_pred(z); - } - else - { - implied_cpuppa(m,i,k,z) = -9; - implied_prop_Cmax(m,i,k,z) = -9; - } - - /*cout <<"year is "<5||F1<0.01)) - { - ii=9; - count_Ffail++; - cout<5||F1<0.01)) - { - ii=5; - count_Ffail++; - cout< 0) - rec_like(5) += 10.*norm2(log_rec_devs(endyr_est,endyr_r))/(2.*sigmarsq_out+.001);// WILL BREAK ON RETROSPECTIVE - - /* Larval drift contribution to recruitment prediction (not used in recent years) Eq. 8 - if (active(larv_rec_devs)) - rec_like(3) += ctrl_flag(23)*norm2(larv_rec_devs); - if (ctrl_flag(27)>0) - { - larv_rec_trans=trans(larv_rec_devs); - // Do third differencing on spatial aspect... - for (i=1;i<=11;i++) - { - rec_like(6) += ctrl_flag(27)*norm2(first_difference( first_difference(first_difference( larv_rec_devs(i))))); - rec_like(6) += ctrl_flag(27)*norm2(first_difference( first_difference(first_difference(larv_rec_trans(i))))); - } - } - */ - - // +===+====+==+==+==+==+==+==+==+====+====+==+==+===+====+==+==+==+==+==+==+==+====+====+====+ -FUNCTION Evaluate_Objective_Function - // if (active(repl_F)) - Recruitment_Likelihood(); - Surv_Likelihood(); //-survey Deviations - Selectivity_Likelihood(); - - catch_like = norm2(log(obs_catch(styr,endyr_r)+1e-4)-log(pred_catch+1e-4)); - - if (current_phase() >= robust_phase) - Robust_Likelihood(); //-Robust AGE Likelihood part - else - Multinomial_Likelihood(); //-Multinomial AGE Likelihood part - - NLL.initialize(); - NLL(1) += ctrl_flag(1) * catch_like; - NLL(2) += ctrl_flag(2) * surv_like(1); - NLL(3) += ctrl_flag(2) * surv_like(2); - NLL(4) += ctrl_flag(2) * surv_like(3); - NLL(5) += ctrl_flag(12) * cpue_like; - NLL(6) += ctrl_flag(6) * avo_like; - NLL(7) += ctrl_flag(3) * sum(rec_like); - if (phase_cope>0 & current_phase()>=phase_cope) - NLL(8) += cope_like; - F_pen = norm2(log_F_devs); - NLL(9) += ctrl_flag(4) * F_pen; - - NLL(10) += ctrl_flag(7)*age_like(1); - NLL(11) += ctrl_flag(8)*age_like(2); - NLL(12) += ctrl_flag(9)*age_like(3); - - if (use_endyr_len>0) - NLL(13)+= ctrl_flag(7)*len_like; - - NLL(14)+= sum(sel_like); - NLL(15)+= sum(sel_like_dev); - - //COUT(sel_like_dev); - // COUT(age_like); - // COUT(avo_like); - // COUT(surv_like); - - // Condition model in early phases to stay reasonable - // Removed at the end - if (current_phase()<3) - { - fff += 10.*square(log(mean(Fmort)/.2)); - fff += 10.*square(log_avginit-log_avgrec) ; //This is to make the initial comp not stray too far - } - - Priors.initialize(); - - if (active(natmort_phi)) // Sensitivity approach for estimating natural mortality (as offset of input vector, NOT USED, NOT IN DOC) - Priors(3) = norm2( log(natmort(3,nages) / natmortprior) )/(2.*cvnatmortprior*cvnatmortprior); - - // Prior on combined-survey catchability, idea is to have sum of two surveys tend to the prior distribution - // q_all.initialize(); - dvariable q_bts_tmp; - dvariable q_ats_tmp; - q_bts_tmp.initialize(); - q_ats_tmp.initialize(); - - for (i=1;i<=n_bts_r;i++) - { - iyr = yrs_bts_data(i); - // Note this is to correct for reduced survey strata coverage pre 1985 and in 86 - if (!(iyr<1985||iyr==1986)) - { - q_bts_tmp += sum(mfexp(log_q_bts + log_sel_bts(iyr)(q_amin,q_amax))); - } - } - q_bts_tmp /= ((q_amax-q_amin+1)*(n_bts_r-4)) ; // 4 years not in main q calc...OjO will break if BTS series length changes between 1982 and 86 - for ( i=1;i<=n_ats_r;i++) - { - iyr = yrs_ats_data(i); - q_ats_tmp += sum(mfexp(log_q_ats + log_sel_ats(iyr)(q_amin,q_amax))); - } - q_ats_tmp /= ((q_amax-q_amin+1)*n_ats_r) ; - q_all= log(q_bts_tmp + q_ats_tmp) ; - - // Note: optional ability to put a prior on "combined" surveys overall catchability - if (q_all_sigma<1.) - Priors(2) = square( q_all- q_all_prior )/(2.*q_all_sigma*q_all_sigma); - q_all = exp(q_all); - - // Prior on BTS catchability - if (active(log_q_bts)&&q_bts_sigma<1.) - { - Priors(2) = square( log_q_bts - q_bts_prior )/(2.*q_bts_sigma*q_bts_sigma); - // cout<log_sel_fsh(i,j+1)) - sel_like(1)+=ctrl_flag(13)*square(log_sel_fsh(i,j)-log_sel_fsh(i,j+1)); - - - if (active(sel_coffs_bts)) - { - for (i=styr_bts;i<= endyr_r;i++) //--This is for controlling the selectivity shape BOTTOM TRAWL SURVEY - for (j=6;j<=n_selages_bts;j++) - if (log_sel_bts(i,j)>log_sel_bts(i,j+1)) - sel_like(2)+=ctrl_flag(14)*square(log_sel_bts(i,j)-log_sel_bts(i,j+1)); - } - - for (i=styr_ats;i<= endyr_r;i++) //--This is for selectivity shape HYDROA TRAWL SURVEY - for (j=mina_ats;j<=n_selages_ats;j++) - if (log_sel_ats(i,j)0.){ - dvar_matrix lnseltmp = trans(log_sel_bts); - for (j=q_amin;j3) - { - for (i=1;i<=n_bts_r;i++) - { - // Under development--not used (yet). Idea to combine surveys for directly accounting for differences - // of water column densities... - surv_like(1) += square(ot_cmb(i)-et_cmb(i))/(2.*var_cmb(i)); - // slight penalty here to get the q to scale correctly (note--development, not used) - surv_like(1) += .01*square(ot_bts(i)-et_bts(i))/(2.*var_ot_bts(i)); - } - // slight penalty here to get the q to scale correctly - for (i=1;i<=n_ats_r;i++) - surv_like(2) += .01*square(log(ot_ats(i)+.01)-log(et_ats(i)+.01))/ (2.*lvar_ats(i)) ; - } - else - { - // This is used (standard approach) Eq. 19, historically used normal distribution since year-by-year var_bts avail - // for (i=1;i<=n_bts_r;i// ++) - // surv_like(1) += square(ot_bts(i)-et_bts(i))/(2.*var_bts(i)); - - //dvar_vector srv_tmp = log(ot_bts + 1e-8)-log(et_bts + 1e-8); - // Note not logged... - - dvar_vector srv_tmp(1,n_bts_r); - // eb_bts *= q_bts; - // NOTE for VAST estimates the biomass is simply scaled - q_bts =mean(ob_bts)/mean(eb_bts); - eb_bts *= q_bts; - if (do_bts_bio) - srv_tmp = (ob_bts) - (eb_bts ); - else - srv_tmp = (ot_bts )-(et_bts ); - - switch (DoCovBTS) - { - case 0: // normal design-based estimates (no Covariance) - if (do_bts_bio) - { - srv_tmp = log(ob_bts) - log(eb_bts ); - for (i=1;i<=n_bts_r;i++) - surv_like(1) += square(srv_tmp(i))/(2.*var_ob_bts(i)); - } - - case 1: - surv_like(1) = .5 * srv_tmp * inv_bts_cov * srv_tmp; - break; - - case 2: // Test gen gamma - if (do_bts_bio) - { - for (i=1;i<=n_bts_r;i++) - if (sd_GenGam(i)>0) - surv_like(1) += dgengamma(ob_bts(i), eb_bts(i), sd_GenGam(i), q_GenGam(i)); - else{ - surv_like(1) += square(log(ob_bts(i))- log(eb_bts(i)))/ (2*var_ob_bts(i)); - } - } - break; - case 3: // lognormal - if (do_bts_bio) - { - // Need to fix up variance for this option... - srv_tmp = log(ob_bts) - log(eb_bts ); - } - break; - } - } - /* else - { - if (do_bts_bio) - { - for (i=1;i<=n_bts_r;i++) - surv_like(1) += square(srv_tmp(i))/(2.*var_ob_bts(i)); - } - else - { - for (i=1;i<=n_bts_r;i++) - surv_like(1) += square(srv_tmp(i))/(2.*var_ot_bts(i)); - } - } - surv_like(1) *= ctrl_flag(5); - */ - - // AT Biomass section - /* - */ - if (do_ats_bio) - { - for (i=1;i<=n_ats_r;i++) // Eq. 19 - surv_like(2) += square(log(ob_ats(i)+.01)-log(eb_ats(i)+.01))/ (2.*lvarb_ats(i)) ; - // #surv_like(2) += square(log(ob_ats(i)+.01)-log(eb_ats(i)+.01))/ (2.*var_ob_ats(i)) ; - } - else - { - for (i=1;i<=n_ats_r;i++) // Eq. 19 - surv_like(2) += square(log(ot_ats(i)+.01)-log(et_ats(i)+.01))/ (2.*lvar_ats(i)) ; - } - surv_like(2) *= ctrl_flag(2); - - if (use_age1_ats) - { - // Compute q for this age1 index... - dvariable qtmp = mfexp(mean(log(oa1_ats)-log(ea1_ats))); - if (ignore_last_ats_age1) - surv_like(3) = 0.5*norm2(log(oa1_ats(1,n_ats_r-1)+.01)-log(ea1_ats(1,n_ats_r-1)*qtmp +.01))/(age1_sigma_ats*age1_sigma_ats) ; - else - surv_like(3) = 0.5*norm2(log(oa1_ats+.01)-log(ea1_ats*qtmp +.01))/(age1_sigma_ats*age1_sigma_ats) ; - } - avo_like.initialize(); - cpue_like.initialize(); - cope_like.initialize(); - - dvar_vector cpue_dev = obs_cpue-pred_cpue; - for (i=1;i<=n_cpue;i++) - cpue_like += square(cpue_dev(i))/(2.*obs_cpue_var(i)); - - dvar_vector avo_dev = obs_avo-pred_avo; - for (i=1;i<=n_avo_r;i++) - avo_like += square(avo_dev(i))/(2.*obs_avo_var(i)); - - // if (phase_cope>0 & current_phase()>=phase_cope) - if (last_phase()) - { - if (phase_cope>0) - { - // Compute q for this age1 index... - int ntmp = n_cope - (yrs_cope(n_cope)+3-endyr_r); - dvariable qtmp = mfexp(mean(log(obs_cope(1,ntmp))-log(pred_cope(1,ntmp)))); - for (i=ntmp+1;i<=n_cope;i++) - pred_cope(i) = obs_cope(i)/qtmp; - pred_cope *= qtmp; - - // dvar_vector cope_dev = obs_cope-pred_cope; - for (i=1;i<=n_cope;i++) - cope_like += square(log(obs_cope(i))-log(pred_cope(i)))/ (2.*lvar_cope(i)) ; - // cope_like += square(cope_dev(i))/(2.*square(obs_cope_std(i))); - } - } - /* // This is for projecting Nage_3 but left out for now - else - { - if (sd_phase()) - { - for (i=1;i<=n_cope;i++) - { - if (yrs_cope(i)>endyr_r) - Nage_3(i) = natage_future(3,yrs_cope(i),3); - else - Nage_3(i) = natage(yrs_cope(i),3); - } - } - } - */ - if (sd_phase()) - Nage_3 = column(natage,3); - -FUNCTION Robust_Likelihood - age_like.initialize(); - len_like.initialize(); - // dvariable robust_p(const dmatrix& obs,const dvar_matrix& pred,const dvariable& a,const dvariable& b) - dvariable rf=.1/nages; - age_like(1) = robust_p(oac_fsh,eac_fsh,rf,sam_fsh); - age_like(2) = robust_p(oac_bts,eac_bts,rf,sam_bts); - - - if (current_phase() >= ats_robust_phase ) // ats robustness phase big number means do multinomial, not robust - age_like(3) = robust_p(oac_ats,eac_ats,rf,sam_ats,mina_ats,nages); - else // Multinomial for EIT - for (i=1; i <= nagecomp(3); i++) - age_like(3) -= sam_ats(i)*oac_ats(i)(mina_ats,nages) * log(eac_ats(i)(mina_ats,nages) + MN_const); - - len_like = robust_p(olc_fsh,elc_fsh,rf,50); - -FUNCTION Multinomial_Likelihood - age_like.initialize(); - len_like.initialize(); - //-Likelihood due to Age compositions-------------------------------- - for (int igear =1;igear<=ngears;igear++) - { - for (i=1; i <= nagecomp(igear); i++) - { - switch (igear) - { - case 1: - age_like(igear) -= sam_fsh(i)*oac_fsh(i)*log(eac_fsh(i) + MN_const); - break; - case 2: - age_like(igear) -= sam_bts(i)*oac_bts(i)*log(eac_bts(i) + MN_const); - break; - case 3: - age_like(igear) -= sam_ats(i)*oac_ats(i)(mina_ats,nages) * log(eac_ats(i)(mina_ats,nages) + MN_const); - break; - } - } - age_like(igear)-=age_like_offset(igear); - } - if (sd_phase() || mceval_phase()) // This to write out annual likelihood components - { - for (int igear =1;igear<=ngears;igear++) - { - for (i=1; i <= nagecomp(igear); i++) - { - switch (igear) - { - case 1: - age_like_yr(igear,i) -= sam_fsh(i)*oac_fsh(i)*log(eac_fsh(i) + MN_const); - break; - case 2: - age_like_yr(igear,i) -= sam_bts(i)*oac_bts(i)*log(eac_bts(i) + MN_const); - break; - default: - age_like_yr(igear,i) -= sam_ats(i)*oac_ats(i)(mina_ats,nages)*log(eac_ats(i)(mina_ats,nages) +MN_const); - break; - } - } - } - } - //len_like = sam_fsh(n_fsh_r)*olc_fsh*log(elc_fsh+MN_const); - len_like = -50*olc_fsh*log(elc_fsh+MN_const) - len_like_offset ; - - // this one allows a concentrated range of ages (last two args are min and max age range) -FUNCTION dvariable robust_p(dmatrix& obs,dvar_matrix& pred,const dvariable& a, const data_ivector& b, const int amin, const int amax) - if (obs.indexmin() != pred.indexmin() || obs.indexmax() != pred.indexmax() ) - cerr << "Index limits on observed vector are not equal to the Index\n" - "limits on the predicted vector in robust_p function\n"; - RETURN_ARRAYS_INCREMENT(); //Need this statement because the function returns a variable type - dvar_matrix v(obs.indexmin(),obs.indexmax(),amin,amax); - // dvar_matrix l = elem_div(square(pred - obs), v); - dvariable log_likelihood = 0.; - for (i=obs.indexmin();i<= obs.indexmax() ;i++) - { - v(i) = a + 2. * elem_prod(obs(i)(amin,amax) ,1. - obs(i)(amin,amax)); - dvar_vector l = elem_div(square(pred(i)(amin,amax) - obs(i)(amin,amax)), v(i)); - log_likelihood -= sum(log(mfexp(-1.* double(b(i)) * l) + .01)); - } - log_likelihood += 0.5 * sum(log(v)); - RETURN_ARRAYS_DECREMENT(); // Need this to decrement the stack increment - return(log_likelihood); - -FUNCTION dvariable robust_p(const dmatrix& obs,const dvar_matrix& pred,const dvariable& a, const data_ivector& b) - RETURN_ARRAYS_INCREMENT(); //Need this statement because the function returns a variable type - if (obs.indexmin() != pred.indexmin() || obs.indexmax() != pred.indexmax() ) - cerr << "Index limits on observed vector are not equal to the Index\n" - "limits on the predicted vector in robust_p function\n"; - // dvar_matrix v = a + 2. * elem_prod(pred ,1. - pred ); - dvar_matrix v = a + 2. * elem_prod(obs ,1. - obs ); - dvar_matrix l = elem_div(square(pred - obs), v); - dvariable log_likelihood = 0.; - for (i=obs.indexmin();i<= obs.indexmax() ;i++) - { - log_likelihood -= sum(log(mfexp(-1.* double(b(i)) * l(i)) + .01)); - } - log_likelihood += 0.5 * sum(log(v)); - RETURN_ARRAYS_DECREMENT(); // Need this to decrement the stack increment - return(log_likelihood); - -FUNCTION dvariable robust_p(const dmatrix& obs, const dvar_matrix& pred, const dvariable& a, const dvariable& b) - RETURN_ARRAYS_INCREMENT(); //Need this statement because the function - if (obs.indexmin() != pred.indexmin() || obs.indexmax() != pred.indexmax() ) - cerr << "Index limits on observed vector are not equal to the Index\n" - "limits on the predicted vector in robust_p function\n"; - //returns a variable type - - // dvar_matrix v = a + 2. * elem_prod(pred ,1. - pred ); - dvar_matrix v = a + 2. * elem_prod(obs ,1. - obs ); - dvar_matrix l = mfexp(- b * elem_div(square(pred - obs), v )); - dvariable log_likelihood = -1.0*sum(log(l + .01)); - log_likelihood += 0.5 * sum(log(v)); - RETURN_ARRAYS_DECREMENT(); // Need this to decrement the stack increment - return(log_likelihood); - -FUNCTION dvariable robust_p(const dvector& obs, const dvar_vector& pred, const dvariable& a, const dvariable& b) - RETURN_ARRAYS_INCREMENT(); //Need this statement because the function - if (obs.indexmin() != pred.indexmin() || obs.indexmax() != pred.indexmax() ) - cerr << "Index limits on observed vector are not equal to the Index\n" - "limits on the predicted vector in robust_p function\n"; - //returns a variable type - - // dvar_matrix v = a + 2. * elem_prod(pred ,1. - pred ); - dvar_vector v = a + 2. * elem_prod(obs ,1. - obs ); - dvar_vector l = mfexp(- b * elem_div(square(pred - obs), v )); - dvariable log_likelihood = -1.0*sum(log(l + .01)); - log_likelihood += 0.5 * sum(log(v)); - RETURN_ARRAYS_DECREMENT(); // Need this to decrement the stack increment - return(log_likelihood); - -FUNCTION write_srec - srecout << "Yearclass SSB Recruit Pred_Rec Residual"<0) - SimulateData1(); - else - { - // if (!pflag) - /* - for (int k=1;k<=nscen;k++) - { - write_mceval(future_SSB(k)); - write_mceval(future_catch(k)); - } - write_mceval(Fcur_Fmsy); - write_mceval(Bcur_Bmsy); - write_mceval(Bcur_Bmean); - write_mceval(Bcur2_Bmsy); - write_mceval(Bcur2_B20); - write_mceval(Bcur3_Bcur); - write_mceval(Bcur3_Bmean); - write_mceval(LTA1_5R); - write_mceval(LTA1_5); - write_mceval(MatAgeDiv1); - write_mceval(MatAgeDiv2); - write_mceval(RelEffort); - write_mceval <=1;iyr--) - { - sel_fut = sel_fsh(endyr_r-iyr+1); - sel_fut /= mean(sel_fut); // NORMALIZE TO mean - report<=1;iyr--) - { - wt_fut(3,nages) = wt_pre(endyr_r-iyr+1); - report< age1tmp) report << age1tmp-1 << " " < age1tmp) report << age1tmp-1 << " " < age1tmp) report << age1tmp-1 << " " < age1tmp) report << age1tmp-1 << " " < age1tmp) report << age1tmp-1 << " " < age1tmp) report << age1tmp-1 << " " <LENGTH-------------------------------------------------- -FUNCTION Get_Age2length - // Linf=Linfprior;// Asymptotic length - // k_coeff=kprior; - // Lo=Loprior;// first length (corresponds to first age-group) - // sdage=sdageprior;// coefficient of variation of length-at-age - // if some of these are estimated. - // L1 25.23382299 k 0.139339914 Linf 68.43045172 - double Linf = 68.43045172 ; - double k_coeff = 0.139339914 ; - double Lo = 25.23382299 ; - double sdage = .06; - dvar_vector mu_age(1,nages); - dvar_vector sigma_age(1,nages); - int i, j; - mu_age(1)=Lo; // first length (modal) - for (i=2;i<=nages;i++) - mu_age(i) = Linf*(1.-exp(-k_coeff))+exp(-k_coeff)*mu_age(i-1); // the mean length by age group - // wt_age_vb(r) = lw_a * pow(mu_age, lw_b); - // maturity_vb(r) = 1/(1 + exp(32.93 - 1.45*mu_age(r))); - sigma_age = sdage*mu_age; // standard deviation of length-at-age - age_len = value(Age_Len_Conversion( mu_age, sigma_age, lens)); - write_log(sigma_age); - -FUNCTION dvar_matrix Age_Len_Conversion(dvar_vector& mu, dvar_vector& sig, dvector& x) - //RETURN_ARRAYS_INCREMENT(); - int i, j; - dvariable z1; - dvariable z2; - int si,ni; si=mu.indexmin(); ni=mu.indexmax(); - int sj,nj; sj=x.indexmin(); nj=x.indexmax(); - dvar_matrix pdf(si,ni,sj,nj); - double xs; - pdf.initialize(); - for(i=si;i<=ni;i++) //loop over ages - { - for(j=sj;j<=nj;j++) //loop over length bins - { - if (jint(1) ) - wt_hat(h,i) = elem_prod(d_scale(h-1) , wt_pre(iyr) ); - else - wt_hat(h,i) = wt_pre(iyr); - - for (int j=age_st;j<=age_end;j++) - { - wt_like += square(wt_obs(h,i,j) - mnwt(j)) /(2.*square(sd_obs(h,i,j))); - wt_like += square(wt_obs(h,i,j) - wt_hat(h,i,j))/(2.*square(sd_obs(h,i,j))); - } - } - } - wt_like += 0.5*norm2(coh_eff); - wt_like += 0.5*norm2( yr_eff); - fff += wt_like; - wt_last = wt_pre(endyr_r-1); //*exp(sigma_coh*sigma_coh/2. + sigma_yr*sigma_yr/2.);; - wt_cur = wt_pre(endyr_r ); //*exp(sigma_coh*sigma_coh/2. + sigma_yr*sigma_yr/2.);; - wt_next = wt_pre(endyr_r+1); //*exp(sigma_coh*sigma_coh/2. + sigma_yr*sigma_yr/2.);; - wt_yraf = wt_pre(endyr_r+2); //*exp(sigma_coh*sigma_coh/2. + sigma_yr*sigma_yr/2.);; - - // Condition on using this based on wt flag - if (wt_fut_phase>0) - { - // Use cohort and year effects fits for current year catch - pred_catch(endyr_r) = catage(endyr_r)(3,nages) * wt_cur; - - // Set future catch equal to estimate from model - // Only model wts-at-age from 3+ so this is the 1's and 2's - pred_catch(endyr_r) = catage(endyr_r)(1,2) * wt_fsh(endyr_r)(1,2); - // wt_fut = wt_fsh(endyr_r); // initializes estimates to correct values...Eq. 21 - wt_fut(3,nages) = wt_next; // initializes estimates to correct values...Eq. 21 - } -FUNCTION dvariable dgengamma(const double& x, dvariable mean2, const double& sigma, const double& Q) - // returns negative log-likelihood of generalized gamma distribution - RETURN_ARRAYS_INCREMENT(); - double k = pow( Q, -2 ); - double Beta = pow( sigma, -1 ) * Q; - dvariable log_theta = log(mean2) - lgamma( (k*Beta+1)/Beta ) + lgamma( k ); - dvariable mu = log_theta + log(k) / Beta; - dvariable w = (log(x) - mu) / sigma; - double abs_q = sqrt(Q*Q); // = abs(Q); not differentiable! - double qi = 1/square(Q); - dvariable qw = Q*w; - dvariable logres = -log(sigma*x) + log(abs_q) * (1. - 2. * qi) + qi * - (qw - exp(qw)) - lgamma(qi); - RETURN_ARRAYS_DECREMENT(); - return(-logres); - -FUNCTION double sdnr(const dvector& obs, const dvar_vector& pred, const dvector& sig) - RETURN_ARRAYS_INCREMENT(); - double sdnr; - sdnr = std_dev(elem_div((obs-value(pred)),sig)); - RETURN_ARRAYS_DECREMENT(); - return sdnr; - - /* ****************************************************** - FUNCTION dvariable Check_Parm(const double& Pmin, const double& Pmax, const double& jitter, const prevariable& Pval) - { - dvariable NewVal; - dvariable temp; - NewVal=Pval; - if(PvalPmax) - {N_warn++; warning<<" parameter init value is greater than parameter max "< "<0.0) - { - temp=log((Pmax-Pmin+0.0000002)/(NewVal-Pmin+0.0000001)-1.)/(-2.); // transform the parameter - temp += randn(radm) * jitter; - NewVal=Pmin+(Pmax-Pmin)/(1.+mfexp(-2.*temp)); - } - return NewVal; - } - */ - - /** - * @brief Calculate sdnr and MAR - **/ - /** - FUNCTION void get_all_sdnr_MAR() - { - for ( int k = 1; k <= nSurveys; k++ ) - { - dvector stdtmp = cpue_sd(k) * 1.0 / cpue_lambda(k); - dvar_vector restmp = elem_div(log(elem_div(obs_cpue(k), pre_cpue(k))), stdtmp) + 0.5 * stdtmp; - sdnr_MAR_cpue(k) = calc_sdnr_MAR(value(restmp)); - } - for ( int k = 1; k <= nSizeComps; k++ ) - { - //dvector sdtmp = cpue_sd(k) * 1.0 / cpue_lambda(i); - sdnr_MAR_lf(k) = calc_sdnr_MAR(value(d3_res_size_comps(k))); - } - Francis_weights = calc_Francis_weights(); - } - - - FUNCTION dvector calc_sdnr_MAR(dvector tmpVec) - { - dvector out(1,2); - dvector tmp = fabs(tmpVec); - dvector w = sort(tmp); - out(1) = std_dev(tmpVec); // sdnr - out(2) = w(floor(0.5*(size_count(w)+1))); // MAR - return out; - } - - - FUNCTION dvector calc_sdnr_MAR(dmatrix tmpMat) - { - dvector tmpVec(1,size_count(tmpMat)); - dvector out(1,2); - int dmin,dmax; - dmin = 1; - for ( int ii = tmpMat.indexmin(); ii <= tmpMat.indexmax(); ii++ ) - { - dmax = dmin + size_count(tmpMat(ii)) - 1; - tmpVec(dmin,dmax) = tmpMat(ii).shift(dmin); - dmin = dmax + 1; - } - dvector tmp = fabs(tmpVec); - dvector w = sort(tmp); - out(1) = std_dev(tmpVec); // sdnr - out(2) = w(floor(0.5*(size_count(w)+1))); // MAR - return out; - } - **/ - -FUNCTION dvector rmultinomial(const dvector ptmp, const int n_sam) - int i1=ptmp.indexmin(); - int i2=ptmp.indexmax(); - dvector freq(i1,i2); - dvector p(i1,i2); - ivector bin(1,n_sam); - p = ptmp; - p /= sum(p); - bin.fill_multinomial(rng,p); // fill a vector v - for (int j=1;j<=n_sam;j++) - freq(bin(j))++; - // Apply ageing error to samples.............. - return( freq/sum(freq) ); - /** - * @brief Calculate Francis weights - * @details this code based on equation TA1.8 in Francis(2011) should be changed so separate weights if by sex - * - * Produces the new weight that should be used. - **/ -FUNCTION double calc_Francis_weights(const dmatrix oac, const dvar_matrix eac, const ivector sam ) - { - int nobs; - int i1=oac.rowmin(); - int i2=oac.rowmax(); - double lfwt,Var,Pre,Obs; - dvector ages(oac.colmin(),nages); - for (int i=oac.colmin();i<=nages;i++) - ages(i) = double(i)+.5; - nobs = oac.rowsize(); - dvector resid(i1,i2); - resid.initialize(); - for ( int i = i1; i <= i2; i++ ) - { - // Obs = sum(elem_prod(oac(i), ages+.5)); - Obs = oac(i) * (ages+.5); - // Pre = sum(elem_prod(value(eac(i)), ages+.5)); - Pre = value(eac(i)) * (ages+.5); - Var = value(eac(i)) * square(ages+.5); - Var -= square(Pre); - resid(i) = (Obs - Pre) / sqrt(Var * 1.0 / (sam(i) )); - // cout<<"FW "<0) - SimulateData1(); - //cout<< "Estimated and SR-predicted recruits"<0) - report<< -dgengamma(ob_bts(i), eb_bts(i), sd_GenGam(i), q_GenGam(i))< - #include - - #undef COUT - #define COUT(object) cout << #object "," << object << endl; - - ofstream misc_out("misc_out.rep"); - #undef misc_out - #define misc_out(object) misc_out << #object "\n" << object << endl; - - ofstream for_sd("extra_sd.rep"); - #undef for_sd - #define for_sd(object) for_sd << #object "\n" << object << endl; - - ofstream write_mcFmort("mcFmort.rep"); - #undef write_mcFmort - #define write_mcFmort(object) write_mcFmort << " " << object ; - - ofstream write_mcSRR("mcSRR.rep"); - #undef write_mcSRR - #define write_mcSRR(object) write_mcSRR << " " << object ; - - ofstream write_mceval("mceval.rep"); - #undef write_mceval - #define write_mceval(object) write_mceval << " " << object ; - - ofstream mceval_ac_ppl("mceval_ac_ppl.csv"); - #undef write_mceval_ac_ppl - #define write_mceval_ac_ppl(object) mceval_ac_ppl << "," << object ; - - - ofstream mceval_ppl("mceval_ppl.csv"); - #undef write_mceval_ppl - #define write_mceval_ppl(object) mceval_ppl << "," << object ; - - ofstream write_mceval_eb_bts("mceval_eb_bts.rep"); - #undef write_mceval_eb_bts - #define write_mceval_eb_bts(object) write_mceval_eb_bts << " " << object ; - - ofstream write_mceval_eb_ats("mceval_eb_ats.rep"); - #undef write_mceval_eb_ats - #define write_mceval_eb_ats(object) write_mceval_eb_ats << " " << object ; - - ofstream write_mceval_ea1_ats("mceval_ea1_ats.rep"); - #undef write_mceval_ea1_ats - #define write_mceval_ea1_ats(object) write_mceval_ea1_ats << " " << object ; - - ofstream write_mceval_pred_catch("mceval_pred_catch.rep"); - #undef write_mceval_pred_catch - #define write_mceval_pred_catch(object) write_mceval_pred_catch << " " << object ; - - ofstream write_mceval_eac_fsh("mceval_eac_fsh.rep"); - #undef write_mceval_eac_fsh - #define write_mceval_eac_fsh(object) write_mceval_eac_fsh << " " << object ; - - ofstream write_mceval_eac_bts("mceval_eac_bts.rep"); - #undef write_mceval_eac_bts - #define write_mceval_eac_bts(object) write_mceval_eac_bts << " " << object ; - - ofstream write_mceval_eac_ats("mceval_eac_ats.rep"); - #undef write_mceval_eac_ats - #define write_mceval_eac_ats(object) write_mceval_eac_ats << " " << object ; - - ofstream write_mceval_pred_cpue("mceval_pred_cpue.rep"); - #undef write_mceval_pred_cpue - #define write_mceval_pred_cpue(object) write_mceval_pred_cpue << " " << object ; - - ofstream write_mceval_pred_avo("mceval_pred_avo.rep"); - #undef write_mceval_pred_avo - #define write_mceval_pred_avo(object) write_mceval_pred_avo << " " << object ; - - ofstream write_mceval_mnwt("mceval_mnwt.rep"); - #undef write_mceval_mnwt - #define write_mceval_mnwt(object) write_mceval_mnwt << " " << object ; - - ofstream write_mceval_fsh_wt("mceval_fsh_wt.rep"); - #undef write_mceval_fsh_wt - #define write_mceval_fsh_wt(object) write_mceval_fsh_wt << " " << object ; - - ofstream write_mceval_srv_wt("mceval_srv_wt.rep"); - #undef write_mceval_srv_wt - #define write_mceval_srv_wt(object) write_mceval_srv_wt << " " << object ; - - ofstream write_log("Input_Log.rep"); - #undef write_log - #define write_log(object) write_log << #object "\n" << object << endl; - - ofstream retro_out("retro.rep",ios::app); - #undef write_retro - #define write_retro(object) retro_out << #object " " <% ggplot(aes(x=Age,y=sel,color=Model)) + geom_line(linewidth=1.5) + ggthemes::theme_few() + ylab("Selectivity") + scale_x_continuous(breaks=1:15);p1 - p1 <- plot_sel(sel=modtune[[2]]$sel_bts,styr=1982,fill="darkblue") ;p1 + library(viridis) + library(scales) + library(ggridges) - plot_agefit(modlst[[1]],case_label=af_title,gear="bts",type="survey",styr=1982,ageplus=10) + p1 <- plot_sel(sel=modlst[[7]]$sel_bts,styr=1982,fill="darkblue") ;p1 + p1 <- plot_sel(sel=modlst[[7]]$sel_ats,styr=1994,fill="darkblue") ;p1 - plot_agefit(modlst[[2]],case_label=af_title,gear="bts",type="survey",styr=1982,ageplus=10) +af_title<-"2023 Assessment" +plot_agefit(modlst[[7]],case_label=af_title,gear="bts",type="survey",styr=1984,ageplus=10) +plot_agefit(modlst[[7]],case_label=af_title,gear="fsh",type="fishery",styr=1984,ageplus=10) +plot_agefit(modlst[[6]],case_label=af_title,gear="bts",type="survey",styr=1984,ageplus=10) #ggsave("figs/sel_comp_vast.pdf",plot=p1,width=8,height=4.0,units="in") #---Age diversity df <- data.frame(Year=M$Yr,Age=M$H,Measure="Population Age\n diversity") - p1 <- plot_avo(modtune[c(1,2)], ylim=NULL) + xlim(c(2005,2023)) + theme_few(base_size = 18);p1 + p1 <- plot_avo(modlst[c(2)], ylim=NULL) + xlim(c(2005,2024)) + theme_few(base_size = 18);p1 p1 <- plot_avo(modtune[c(1,2)], ylim=NULL) + xlim(c(2005,2023)) + facet_grid(Model~.,scales='free')+theme_few(base_size = 18);p1 #--tuned figures @@ -45,9 +57,12 @@ if (doplots) { # "SSB Emp. wt-age", #10 # "SSB RE wt-age") #11 - p1 <- plot_ssb(modlst[c(1,2)],xlim=c(2008.5,2023.5));p1 - p2 <- plot_recruitment(modlst[c(1,2)],xlim=c(1990.5,2023.5));p2 - p3 <- plot_srr(modlst[c(1,2)],alpha=.2,xlim=c(0,5200),ylim=c(0,70000));p3 + .OVERLAY=TRUE + p1 <- plot_ssb(modlst[c(2,7)],xlim=c(2008.5,2023.5));p1 + p2 <- plot_recruitment(modlst[c(2,7)],xlim=c(2000.5,2023.5)) + + theme_few(base_size = 15);p2 + p3 <- ebswp::plot_srr(modlst[c(1,7)],alpha=.2,xlim=c(0,5200),ylim=c(0,70000)) + + theme_few(base_size = 15);p3 p1 <- plot_bts(modlst[c(8,5)],xlim=c(1981.5,2023.5),ylim=c(0,15500)) ;p1 ggsave("doc/figs/mod_bts_gengam.pdf",plot=p1,width=8,height=5.0,units="in") @@ -107,7 +122,9 @@ p1 <- plot_ats(modlst[c(8,5)],xlim=c(1993.5,2023.5),ylim=c(0,10000)) ;p1 #---Age diversity df <- data.frame(Year=M$Yr,Age=M$H,Measure="Population Age\n diversity") df <- rbind(df,data.frame(Year=M$Yr,Age=M$avg_age_mature,Measure="SSB Age\n diversity")) - p1 <- df %>% filter(Year>1979) %>% ggplot(aes(x=Year,y=Age,color=Measure)) + geom_line(size=2) + theme_few() + scale_x_continuous(limits=c(1980,2022), breaks=seq(1980,2022,5));p1 + p1 <- df %>% filter(Year>1979) %>% ggplot(aes(x=Year,y=Age,color=Measure)) + geom_line(size=2) + + theme_few() + scale_x_continuous(limits=c(1980,2023), breaks=seq(1980,2023,5));p1 + p1 ggsave("figs/age_diversity.pdf",plot=p1,width=8,height=3.0,units="in") #---Recruit #p1 <- plot_recruitment(modlst[c(1,2,3)],xlim=c(2008.5,2019.5));p1 @@ -166,7 +183,7 @@ p1 <- plot_ats(modlst[c(8,5)],xlim=c(1993.5,2023.5),ylim=c(0,10000)) ;p1 #sel #M<- modlst[[thismod]] #---Selectivity----------------- - yr=c(M$Yr,2023);sel<-rbind(M$sel_fsh,M$sel_fut) + yr=c(M$Yr,2024);sel<-rbind(M$sel_fsh,M$sel_fut) p1 <- plot_sel(Year=yr,sel=sel,scale=3); p1 #p1 <- plot_sel();p1 #dtmp <- p1$data %>% filter(Year==2022) @@ -367,6 +384,12 @@ modlst #ggsave("figs/mod_ats_eval2.pdf",plot=p5,height=7,width=6.0,units="in"); #---SSB figure---------------------------- + M <- modlst[[7]] + P <- modlst[[1]] + M$future_catch[4,2:6] + M$future_catch + M$future_SSB[5,2:6 ] + M$future_SSB df <- rbind( rbind(data.frame(Model="This year",Year= M$SSB[,1], SSB=M$SSB[,2], lb=M$SSB[,4], ub=M$SSB[,5]), data.frame(Model="This year",Year= nextyr:(nextyr+4), SSB=M$future_SSB[4,2:6], lb = M$future_SSB[4,2:6] -2*M$future_SSB.sd[4,2:6], ub = M$future_SSB[4,2:6] +2*M$future_SSB.sd[4,2:6])), @@ -376,10 +399,11 @@ data.frame(Model="Last year",Year= (nextyr-1):(nextyr+3), SSB=P$future_SSB[4,2:6 p1 <- ggplot(df,aes(x=Year,y=SSB,ymax=ub,ymin=lb,fill=Model)) + geom_ribbon(alpha=.6) + geom_line() + theme_few() + scale_x_continuous(limits=c(2002,2028),breaks=seq(2002,2028,2)) + ylab("Female spawning biomass (kt)") + - geom_vline(xintercept=2022,col="grey",size=1.2); p1 + geom_vline(xintercept=2023,col="grey",size=1.2); p1 ggsave("figs/proj_ssb.pdf",plot=p1,width=7.4,height=4.5,units="in") #---R/S------------------ + M <- modlst[[7]] nyrs=length(M$SSB[,1]) dt <- data.table(yr=M$SSB[1:nyrs,1],ssb= M$SSB[1:nyrs,2], r=M$R[2:(nyrs-1),2] ) dt <- dt %>% mutate(Year=substr(as.character(yr),3,4),rs= log(dt$r/dt$ssb)) @@ -391,14 +415,18 @@ p1 <- ggplot(df,aes(x=Year,y=SSB,ymax=ub,ymin=lb,fill=Model)) + geom_ribbon(alph p3 <- p1 / p2;p3 ggsave("figs/mod_rs.pdf",plot=p3,width=5.2,height=7.5,units="in") + thismod=7 p1 <- plot_ser(modlst[thismod],xlim=c(1964,thisyr+1),alpha=.7) + scale_x_continuous(breaks=seq(1965,thisyr+1,5)) + p1 ggsave("figs/mod_ser.pdf",plot=p1,width=9.2,height=7.0,units="in"); p1 #---fishing mortality mod_F.pdf----------------------------------------------------------------- df <-data.frame(Year=M$Yr,M$F); names(df) <- c("Year",1:15); df.g <- gather(df,age,F,2:16,-Year) + df p1 <- df.g %>% mutate(age=as.numeric(age)) %>% filter(age<11)%>% ggplot(aes(y=age,x=Year,fill=F)) + geom_tile() + .THEME + ylab("Age")+ geom_contour(aes(z=F),color="darkgrey",size=.5,alpha=.4) + scale_fill_gradient(low = "white", high = "red") + scale_x_continuous(breaks=seq(1965,thisyr,5)) + geom_line(data=df.g[df.g$age=="6",],aes(x=Year,y=F*10)) + annotate("text", label = "Age 6 F (x10)" , x = 2015, y = 1.2, size = 5, colour = "black") + scale_y_continuous(breaks=seq(0,10,1)) + p1 ggsave("figs/mod_F.pdf",plot=p1,width=9.2,height=6.0,units="in") #---Historical assessment retrospectives-------------------------------------------------------- @@ -1381,15 +1409,16 @@ print(cd_tab) df.g %>% filter(Measure=="catch") %>% select(year,scen,val) %>% ggplot(aes(x=year,y=scen)) + .THEME + geom_density_2d() ## VAST figure - idxOut <- read.csv("data/vast.csv",header=TRUE) + idxOut <- read.csv("doc/data/vast.csv",header=TRUE) idxOut - df <- mutate(idxOut, CV = sd/Biomass, sd=sd/1e6, Biomass=Biomass / 1e6, ll=Biomass-2*sd,ul=2*sd+Biomass) %>% + df <- filter(Year!=2020) |> mutate(idxOut, CV = sd/Biomass, sd=sd/1e6, Biomass=Biomass / 1e6, ll=Biomass-2*sd,ul=2*sd+Biomass) %>% select(Model, Area, Year,Biomass, CV,ul,ll) - p1 <- ggplot(df,aes(x=Year,y=Biomass)) + geom_point(position=position_dodge(width=0.5)) + geom_smooth(span=.1,se=F) + + p1 <- ggplot(df |> filter(Year!=2020) ,aes(x=Year,y=Biomass)) + geom_point(position=position_dodge(width=0.5)) + geom_smooth(span=.1,se=F) + geom_errorbar(aes(ymin=ll,ymax=ul),width=.5,position=position_dodge(width=0.5)) + theme_few(base_size=19) + scale_x_continuous(breaks=seq(1982,2022,by=4)) + geom_hline(yintercept=0) + facet_grid(Area~.,scales="free") + ylab("Biomass (millions of t)"); p1 - ggsave("figs/vast_idx.pdf",plot=p1,width=9.4,height=5,units="in") + p1 + ggsave("doc/figs/vast_idx.pdf",plot=p1,width=9.4,height=5,units="in") p1 <- ggplot(df%>%filter(Area=="Both"),aes(x=Year,y=Biomass)) + geom_point(size=3,color="red",position=position_dodge(width=0.5)) + geom_errorbar(aes(ymin=ll,ymax=ul),width=.5,position=position_dodge(width=0.5)) + theme_few(base_size=19) + @@ -1572,16 +1601,17 @@ ggsave("figs/catch.pdf",plot=p1,width=7.5,height=4.5,units="in") # flextable(tsum) #---Roe data------------------ # Table comes from AKFIN: https://akfinbi.psmfc.org/analytics/saw.dll?Dashboard&PortalPath=%2Fshared%2FStock%20Assessment%2F_portal%2FStock%20Assessment&Page=Ianelli%20-%20Monthly%20EBS%20pollock%20production&Action=Navigate - rd <- read_csv("../doc/data/roe.csv") + rd <- read_csv("doc/data/roe.csv") glimpse(rd) - names(rd) <- c("Year","Month","TYPE","SECTOR","SPECIES_GROUP_CODE","FMP_SUBAREA","Prod","PERMITS","t","CONFIDENTIAL_FLAG") + summary(rd) + names(rd) <- c("Year","Month","TYPE","SECTOR","SPECIES_GROUP_CODE","FMP_SUBAREA","Prod","PERMITS","t","CONFIDENTIAL_FLAG") rd <- rd %>% mutate(Season=ifelse(Month<6,"A","B")) rd <- select(rd,Year, Season, Prod,t) %>% group_by(Year,Season,Prod) %>% summarize(t=sum(t)) p1 <- rd %>% filter(Prod=="ROE") %>% ggplot(aes(x=Year,y=t,shape=Season,col=Season)) + geom_point(size=3) + geom_line(size=1.0) + - theme_few(base_size=14) + ylab("Tons of roe produced") + expand_limits(y=0) + - scale_x_continuous(breaks=seq(2000,2022,2)); #p1 + ggthemes::theme_few(base_size=14) + ylab("Tons of roe produced") + expand_limits(y=0) + + scale_x_continuous(breaks=seq(2000,2023,2)); #p1 # + geom_hline(yintercept=mst,linetype="dashed") + geom_hline(yintercept=mbt,linetype="dashed") p1 ggsave("figs/roe.pdf",plot=p1,width=7.5,height=4.5,units="in") @@ -1590,25 +1620,25 @@ ggsave("figs/catch.pdf",plot=p1,width=7.5,height=4.5,units="in") #--- Catch-age estimates for sex catch age -------------------------- # Updated edt <- NULL - for (i in 2021:1991){ - edt <- rbind(edt,read_table(paste0("../data/fishery/sampler/results/Est_",i,".dat"))) + for (i in 2022:1991){ + edt <- rbind(edt,read_table(paste0("~/_mymods/ebs_main/data/fishery/sampler/results/Est_",i,".dat"))) } tot <- edt %>% filter(type=="N",sex<3,stratum<5) %>% mutate(season=ifelse(stratum==1,"A","B"),sex=ifelse(sex==1,"Male",ifelse(sex==2,"Female","Total"))) %>% group_by(year,sex) %>% summarize(Catch=sum(value)) tot$season <-"all" p1 <- edt %>% filter(type=="N",sex<3,stratum<5) %>% mutate(season=ifelse(stratum==1,"A","B"),sex=ifelse(sex==1,"Male",ifelse(sex==2,"Female","Total"))) %>% group_by(year,sex,season) %>% summarize(Catch=sum(value)) %>% rbind(tot) %>% - ggplot(aes(x=year,y=Catch,color=season,shape=sex,linetype=season)) + theme_few(base_size=12) + geom_point(size=3.5) + geom_line(size=1.2) + expand_limits(y=0) + - scale_x_continuous(breaks=seq(1990,2020,2)) + ylab("Catch (thousands)") + xlab("Year") - p1 - ggsave("figs/catch_sex.pdf",plot=p1,width=7.5,height=4.5,units="in") + ggplot(aes(x=year,y=Catch,color=season,shape=sex,linetype=season)) + theme_few(base_size=11) + + geom_point(size=3.5) + geom_line(size=1.2) + expand_limits(y=0) + + scale_y_continuous(label=comma) + scale_x_continuous(breaks=seq(1990,2022,2)) + + ylab("Catch (thousands)") + xlab("Year"); p1 + ggsave("doc/figs/catch_sex.pdf",plot=p1,width=7.5,height=4.5,units="in") # Look at catch in Number vs catch in weight edt_df <- edt |> mutate(Year=year) |> group_by(Year) |> summarise(Number=sum(value)) cdf_df <- df |> group_by(Year) |> summarise(Biomass=sum(Catch)) left_join(edt_df, cdf_df) |> ggplot(aes(x=Number,y=Biomass,label=Year)) + theme_few() + geom_text(size=3) + geom_smooth(method=lm) #left_join(edt_df, cdf_df) |> mutate(Number=Number/mean(Number),Biomass=Biomass/mean(Biomass) ) |> pivot_longer(!Year,names_to="Type",values_to="Catch") |> ggplot(aes(x=Year,y=Catch,color=Type)) + - library(scales) left_join(edt_df, cdf_df) |> pivot_longer(!Year,names_to="Type",values_to="Catch") |> ggplot(aes(x=Year,y=Catch,color=Type)) + scale_y_continuous(label=comma) + geom_line(size=2) + theme_few(base_size = 18) + facet_grid(Type~.,scales="free") + theme(legend.position="none")