diff --git a/.Rbuildignore b/.Rbuildignore index 3d49c9a..086246f 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -7,3 +7,4 @@ src\pm.cpp ^docs$ ^pkgdown$ ^\.github$ +^src$ diff --git a/2023_runs/usehulsoniss/Makefile b/2023_runs/usehulsoniss/Makefile index a9d37fa..931e7af 100644 --- a/2023_runs/usehulsoniss/Makefile +++ b/2023_runs/usehulsoniss/Makefile @@ -2,9 +2,9 @@ EXEC = pm PROJ = spm PROJDIST = ~/_mymods/afsc-assessments/spmR/src/ ifeq ($(SAFE),TRUE) - DIST = ../../source/ + DIST = ../../src/ else - DIST = ../../source/ + DIST = ../../src/ endif ARGS = -nox -iprint 150 ARGS2 = -nox -iprint 150 -binp pm.bar -phase 22 -sdonly diff --git a/2023_runs/usehulsoniss/pm b/2023_runs/usehulsoniss/pm index b22c3ba..5134594 120000 --- a/2023_runs/usehulsoniss/pm +++ b/2023_runs/usehulsoniss/pm @@ -1 +1 @@ -../../source/pm \ No newline at end of file +../../src/pm \ No newline at end of file diff --git a/2023_runs/usehulsoniss/pm.tpl b/2023_runs/usehulsoniss/pm.tpl index 35c02d5..3dc208b 120000 --- a/2023_runs/usehulsoniss/pm.tpl +++ b/2023_runs/usehulsoniss/pm.tpl @@ -1 +1 @@ -../../source/pm.tpl \ No newline at end of file +../../src/pm.tpl \ No newline at end of file diff --git a/src/2022/Makefile b/src/2022/Makefile new file mode 100644 index 0000000..e154650 --- /dev/null +++ b/src/2022/Makefile @@ -0,0 +1,18 @@ + +ifdef ComSpec + RM=del /F /Q + COPY=copy +else + RM=rm -rf + COPY=cp +endif + +all: pm + +pm: pm.tpl + @admb -f pm.tpl + +clean: + @$(RM) pm.cpp + @$(RM) pm.htp + @$(RM) pm.obj diff --git a/src/2022/pm b/src/2022/pm new file mode 100755 index 0000000..0c8501b Binary files /dev/null and b/src/2022/pm differ diff --git a/src/2022/pm.cpp b/src/2022/pm.cpp new file mode 100644 index 0000000..2bd3297 --- /dev/null +++ b/src/2022/pm.cpp @@ -0,0 +1,7181 @@ +#ifdef DEBUG + #ifndef __SUNPRO_C + #include + #include + #endif +#endif + #include + #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 " " < +#endif +#include +#ifdef USE_ADMB_CONTRIBS +#include + +#endif + extern "C" { + void ad_boundf(int i); + } +#include + +model_data::model_data(int argc,char * argv[]) : ad_comm(argc,argv) +{ + adstring tmpstring; + tmpstring=adprogram_name + adstring(".dat"); + if (argc > 1) + { + int on=0; + if ( (on=option_match(argc,argv,"-ind"))>-1) + { + if (on>argc-2 || argv[on+1][0] == '-') + { + cerr << "Invalid input data command line option" + " -- ignored" << endl; + } + else + { + tmpstring = adstring(argv[on+1]); + } + } + } + global_datafile = new cifstream(tmpstring); + if (!global_datafile) + { + cerr << "Error: Unable to allocate global_datafile in model_data constructor."; + ad_exit(1); + } + if (!(*global_datafile)) + { + delete global_datafile; + global_datafile=NULL; + } + *(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) >> Temp_Cons_Dist_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(Temp_Cons_Dist_file); + count_mcmc=0; + count_mcsave=0; + q_amin = 3; q_amax= 15; // age range overwhich q applies (for prior specifications) + selages.allocate(1,15); + selages=1.0;selages(1)=0;selages(2)=0; + avo_sel.allocate(1,15); + avo_sel(1)=0.0; avo_sel(2)=1; avo_sel(3)=1; avo_sel(4)=0.85; avo_sel(5)=0.7; avo_sel(6)=0.55; avo_sel(7)=0.3; avo_sel(8)=0.15; avo_sel(9)=0.05; avo_sel(10)=0.01; avo_sel(11)=0.01; avo_sel(12)=0.01; avo_sel(13)=0.01; avo_sel(14)=0.01; avo_sel(15)=0.01; + Cat_Fut.allocate(1,10); + do_EIT1=1; // flag to carry EIT out in future year (for simulations only) + pflag=0; + pad_srecout = new ofstream("srec_Ass_out.rep"); + pad_projout = new ofstream("pm.prj"); + pad_nofish = new ofstream("nofish.rep"); + pad_projout2 = new ofstream("pmsr.prj"); + pad_eval = new ofstream("pm_eval.rep") ; + 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; + } + phase_selcoffs_ats.allocate("phase_selcoffs_ats"); + phase_selcoffs_ats_dev.allocate("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)); + 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; + 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;} + yrs_ch_fsh.allocate(1,nch_fsh); + sel_ch_sig_fsh.allocate(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;} + 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); + oac_fsh.allocate(1,n_fsh_r,1,nbins); + oac_bts.allocate(1,n_bts_r,1,nbins); + 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; + oac_ats.allocate(1,n_ats_ac_r,1,nbins); + oa1_ats.allocate(1,n_ats_ac_r); + ot_fsh.allocate(1,n_fsh_r); + ot_bts.allocate(1,n_bts_r); + std_ob_bts.allocate(1,n_bts_r); + ob_bts.allocate(1,n_bts_r); + ot_ats.allocate(1,n_ats_r); + ob_ats.allocate(1,n_ats_r); + std_ob_ats.allocate(1,n_ats_r); + lseb_ats.allocate(1,n_ats_r); + lvarb_ats.allocate(1,n_ats_r); + iseed=0; + do_fmort=0; + adj_1.allocate(1,10); + adj_2.allocate(1,10); + SSB_1.allocate(1,10); + SSB_2.allocate(1,10); + 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; +long int lseed=iseed; + pad_rng = new random_number_generator(iseed);; + ad_comm::change_datafile_name(Wtage_file); + log_sd_coh.allocate("log_sd_coh"); + log_sd_yr.allocate("log_sd_yr"); + cur_yr.allocate("cur_yr"); + styr_wt.allocate("styr_wt"); + endyr_wt.allocate("endyr_wt"); + ndat_wt.allocate("ndat_wt"); + nyrs_data.allocate(1,ndat_wt,"nyrs_data"); + yrs_data.allocate(1,ndat_wt,1,nyrs_data,"yrs_data"); + if (endyr_r < endyr) + for (int h=1;h<=ndat_wt;h++) + { + int itmp = 1; + for (i=styr;i<=endyr_r;i++) + { + if (i == yrs_data(h,itmp) ) + itmp++; + } + nyrs_data(h)= itmp-1; + } + age_st.allocate("age_st"); + age_end.allocate("age_end"); + nages_wt = age_end - age_st + 1; + nscale_parm = ndat_wt-1; + wt_obs.allocate(1,ndat_wt,1,nyrs_data,age_st,age_end,"wt_obs"); + sd_obs.allocate(1,ndat_wt,1,nyrs_data,age_st,age_end,"sd_obs"); + if (ndat_wt>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)); + cout<<(endyr_wt)<0) + { + mnCPUE(iyr,ist) = mean(CPUE(iyr,ist)); + mntemp(iyr,ist) = mean(temp(iyr,ist)); + } + } + } + FW_fsh.allocate(1,4); + ad_comm::change_datafile_name(Temp_Cons_Dist_file); + SST.allocate(styr-1,endyr-1,"SST"); + SST_fut = mean(SST(endyr-7,endyr-3)); + n_pred_grp_nonpoll.allocate("n_pred_grp_nonpoll"); + n_pred_grp_poll.allocate("n_pred_grp_poll"); + n_pred_grp = n_pred_grp_nonpoll + n_pred_grp_poll; // the number of predator groups + N_pred.allocate(1,n_pred_grp,styr,endyr,"N_pred"); + nstrata_pred.allocate("nstrata_pred"); + strata.allocate(1,nstrata_pred,"strata"); + n_pred_ages.allocate("n_pred_ages"); + pred_ages.allocate(1,n_pred_ages,"pred_ages"); + poll_dist.allocate(1,n_pred_ages,styr,endyr,1,nstrata_pred,"poll_dist"); + pred_dist_nonpoll.allocate(1,n_pred_grp,styr,endyr,1,nstrata_pred,"pred_dist_nonpoll"); + Npred_bystrata_nonpoll.allocate(styr,endyr,1,n_pred_grp_nonpoll,1,nstrata_pred); + area_pred.allocate(1,nstrata_pred,"area_pred"); + nyrs_cons_nonpoll.allocate(1,n_pred_grp_nonpoll,"nyrs_cons_nonpoll"); + yrs_cons_nonpoll.allocate(1,n_pred_grp_nonpoll,1,nyrs_cons_nonpoll,"yrs_cons_nonpoll"); + obs_cons_nonpoll.allocate(1,n_pred_grp_nonpoll,1,nyrs_cons_nonpoll,"obs_cons_nonpoll"); + oac_cons_nonpoll.allocate(1,n_pred_grp_nonpoll,1,nyrs_cons_nonpoll,1,n_pred_ages,"oac_cons_nonpoll"); + sam_oac_cons_nonpoll_raw.allocate(1,n_pred_grp_nonpoll,1,nyrs_cons_nonpoll,"sam_oac_cons_nonpoll_raw"); + obs_cons_wgt_atage_nonpoll.allocate(1,n_pred_grp_nonpoll,1,nyrs_cons_nonpoll,1,n_pred_ages); + obs_cons_natage_nonpoll.allocate(1,n_pred_grp_nonpoll,1,nyrs_cons_nonpoll,1,n_pred_ages); + obs_cpup_nonpoll.allocate(1,n_pred_grp_nonpoll,1,nyrs_cons_nonpoll,1,n_pred_ages); + C_a_nonpoll.allocate(1,n_pred_grp_nonpoll,"C_a_nonpoll"); + C_b_nonpoll.allocate(1,n_pred_grp_nonpoll,"C_b_nonpoll"); + TCM_nonpoll.allocate(1,n_pred_grp_nonpoll,"TCM_nonpoll"); + TC0_nonpoll.allocate(1,n_pred_grp_nonpoll,"TC0_nonpoll"); + CQ_nonpoll.allocate(1,n_pred_grp_nonpoll,"CQ_nonpoll"); + temp_bystrata.allocate(styr,endyr,1,nstrata_pred,"temp_bystrata"); + mn_wgt_nonpoll.allocate(1,n_pred_grp_nonpoll,styr,endyr,"mn_wgt_nonpoll"); + Y_nonpoll.allocate(1,n_pred_grp_nonpoll); + Z_nonpoll.allocate(1,n_pred_grp_nonpoll); + X_nonpoll.allocate(1,n_pred_grp_nonpoll); + V_nonpoll.allocate(styr,endyr,1,n_pred_grp_nonpoll,1,nstrata_pred); + F_t_nonpoll.allocate(styr,endyr,1,n_pred_grp_nonpoll,1,nstrata_pred); + Cmax_nonpoll.allocate(styr,endyr,1,n_pred_grp_nonpoll,1,nstrata_pred); + Cmax_avg.allocate(1,n_pred_grp,"Cmax_avg"); + atf_wgts.allocate(1,n_pred_grp); + poll_wgts.allocate(1,n_pred_ages); + comp_nr_ub.allocate(1,n_pred_grp_nonpoll); + comp_nr_ub = ivector(nyrs_cons_nonpoll*n_pred_ages); + phase_cope.allocate("phase_cope"); + n_cope.allocate("n_cope"); + yrs_cope.allocate(1,n_cope,"yrs_cope"); + obs_cope.allocate(1,n_cope,"obs_cope"); + obs_cope_std.allocate(1,n_cope,"obs_cope_std"); + lse_cope.allocate(1,n_cope); + lvar_cope.allocate(1,n_cope); + 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 + test_2.allocate("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; + sel_dif2_fsh.allocate(ph_log_fsh2,"sel_dif2_fsh"); + sel_dif1_fsh_dev.allocate(styr,endyr_r,-5,5,phase_logist_fsh_dev,"sel_dif1_fsh_dev"); + sel_a501_fsh_dev.allocate(styr,endyr_r,-5,5,phase_logist_fsh_dev,"sel_a501_fsh_dev"); + sel_trm2_fsh_dev.allocate(styr,endyr_r,-.5,.5,-2,"sel_trm2_fsh_dev"); + SPR_ABC.allocate("SPR_ABC"); + #ifndef NO_AD_INITIALIZE + SPR_ABC.initialize(); + #endif + endyr_N.allocate(1,nages,"endyr_N"); + B_Bnofsh.allocate("B_Bnofsh"); + Bzero.allocate("Bzero"); + Percent_Bzero.allocate("Percent_Bzero"); + Percent_Bzero_1.allocate("Percent_Bzero_1"); + Percent_Bzero_2.allocate("Percent_Bzero_2"); + Percent_B100.allocate("Percent_B100"); + Percent_B100_1.allocate("Percent_B100_1"); + Percent_B100_2.allocate("Percent_B100_2"); + q_temp.allocate(1,5,"q_temp"); + #ifndef NO_AD_INITIALIZE + q_temp.initialize(); + #endif + SPR_OFL.allocate("SPR_OFL"); + F40.allocate("F40"); + F35.allocate("F35"); + SSB.allocate(styr,endyr_r,"SSB"); + sigmarsq_out.allocate("sigmarsq_out"); + #ifndef NO_AD_INITIALIZE + sigmarsq_out.initialize(); + #endif + ftmp.allocate("ftmp"); + #ifndef NO_AD_INITIALIZE + ftmp.initialize(); + #endif + SB0.allocate("SB0"); + #ifndef NO_AD_INITIALIZE + SB0.initialize(); + #endif + SBF40.allocate("SBF40"); + #ifndef NO_AD_INITIALIZE + SBF40.initialize(); + #endif + SBF35.allocate("SBF35"); + #ifndef NO_AD_INITIALIZE + SBF35.initialize(); + #endif + sprpen.allocate("sprpen"); + #ifndef NO_AD_INITIALIZE + sprpen.initialize(); + #endif + F_pen.allocate("F_pen"); + #ifndef NO_AD_INITIALIZE + F_pen.initialize(); + #endif + meanrec.allocate("meanrec"); + #ifndef NO_AD_INITIALIZE + meanrec.initialize(); + #endif + SR_resids.allocate(styr_est,endyr_est,"SR_resids"); + #ifndef NO_AD_INITIALIZE + SR_resids.initialize(); + #endif + Nspr.allocate(1,4,1,nages,"Nspr"); + #ifndef NO_AD_INITIALIZE + Nspr.initialize(); + #endif + sel_fut.allocate(1,nages,"sel_fut"); + natage_future.allocate(1,nscen,styr_fut,endyr_fut,1,nages,"natage_future"); + #ifndef NO_AD_INITIALIZE + natage_future.initialize(); + #endif + rec_dev_future.allocate(styr_fut,endyr_fut,phase_F40,"rec_dev_future"); + F_future.allocate(1,nscen,styr_fut,endyr_fut,1,nages,"F_future"); + #ifndef NO_AD_INITIALIZE + F_future.initialize(); + #endif + Z_future.allocate(styr_fut,endyr_fut,1,nages,"Z_future"); + #ifndef NO_AD_INITIALIZE + Z_future.initialize(); + #endif + S_future.allocate(styr_fut,endyr_fut,1,nages,"S_future"); + #ifndef NO_AD_INITIALIZE + S_future.initialize(); + #endif + catage_future.allocate(styr_fut,endyr_fut,1,nages,"catage_future"); + #ifndef NO_AD_INITIALIZE + catage_future.initialize(); + #endif + avg_rec_dev_future.allocate("avg_rec_dev_future"); + #ifndef NO_AD_INITIALIZE + avg_rec_dev_future.initialize(); + #endif + eac_fsh.allocate(1,n_fsh_r,1,nbins,"eac_fsh"); + #ifndef NO_AD_INITIALIZE + eac_fsh.initialize(); + #endif + elc_fsh.allocate(1,nlbins,"elc_fsh"); + #ifndef NO_AD_INITIALIZE + elc_fsh.initialize(); + #endif + eac_bts.allocate(1,n_bts_r,1,nbins,"eac_bts"); + #ifndef NO_AD_INITIALIZE + eac_bts.initialize(); + #endif + eac_cmb.allocate(1,n_bts_r,1,nbins,"eac_cmb"); + #ifndef NO_AD_INITIALIZE + eac_cmb.initialize(); + #endif + oac_cmb.allocate(1,n_bts_r,1,nbins,"oac_cmb"); + #ifndef NO_AD_INITIALIZE + oac_cmb.initialize(); + #endif + eac_ats.allocate(1,n_ats_ac_r,1,nbins,"eac_ats"); + #ifndef NO_AD_INITIALIZE + eac_ats.initialize(); + #endif + ea1_ats.allocate(1,n_ats_ac_r,"ea1_ats"); + #ifndef NO_AD_INITIALIZE + ea1_ats.initialize(); + #endif + et_fsh.allocate(1,n_fsh_r,"et_fsh"); + #ifndef NO_AD_INITIALIZE + et_fsh.initialize(); + #endif + et_bts.allocate(1,n_bts_r,"et_bts"); + #ifndef NO_AD_INITIALIZE + et_bts.initialize(); + #endif + et_cmb.allocate(1,n_bts_r,"et_cmb"); + #ifndef NO_AD_INITIALIZE + et_cmb.initialize(); + #endif + avail_bts.allocate(1,n_bts_r,"avail_bts"); + #ifndef NO_AD_INITIALIZE + avail_bts.initialize(); + #endif + avail_ats.allocate(1,n_bts_r,"avail_ats"); + #ifndef NO_AD_INITIALIZE + avail_ats.initialize(); + #endif + sigma_cmb.allocate(1,n_bts_r,"sigma_cmb"); + #ifndef NO_AD_INITIALIZE + sigma_cmb.initialize(); + #endif + var_cmb.allocate(1,n_bts_r,"var_cmb"); + #ifndef NO_AD_INITIALIZE + var_cmb.initialize(); + #endif + ot_cmb.allocate(1,n_bts_r,"ot_cmb"); + #ifndef NO_AD_INITIALIZE + ot_cmb.initialize(); + #endif + eb_bts.allocate(1,n_bts_r,"eb_bts"); + #ifndef NO_AD_INITIALIZE + eb_bts.initialize(); + #endif + eb_ats.allocate(1,n_ats_r,"eb_ats"); + #ifndef NO_AD_INITIALIZE + eb_ats.initialize(); + #endif + et_ats.allocate(1,n_ats_r,"et_ats"); + #ifndef NO_AD_INITIALIZE + et_ats.initialize(); + #endif + lse_ats.allocate(1,n_ats_r,"lse_ats"); + #ifndef NO_AD_INITIALIZE + lse_ats.initialize(); + #endif + lvar_ats.allocate(1,n_ats_r,"lvar_ats"); + #ifndef NO_AD_INITIALIZE + lvar_ats.initialize(); + #endif + et_avo.allocate(1,n_avo_r,"et_avo"); + #ifndef NO_AD_INITIALIZE + et_avo.initialize(); + #endif + et_cpue.allocate(1,n_cpue,"et_cpue"); + #ifndef NO_AD_INITIALIZE + et_cpue.initialize(); + #endif + Fmort.allocate(styr,endyr_r,"Fmort"); + #ifndef NO_AD_INITIALIZE + Fmort.initialize(); + #endif + catage.allocate(styr,endyr_r,1,nages,"catage"); + #ifndef NO_AD_INITIALIZE + catage.initialize(); + #endif + pred_catch.allocate(styr,endyr_r,"pred_catch"); + #ifndef NO_AD_INITIALIZE + pred_catch.initialize(); + #endif + Pred_N_bts.allocate(styr,endyr_r,"Pred_N_bts"); + #ifndef NO_AD_INITIALIZE + Pred_N_bts.initialize(); + #endif + Pred_N_ats.allocate(styr,endyr_r,"Pred_N_ats"); + #ifndef NO_AD_INITIALIZE + Pred_N_ats.initialize(); + #endif + pred_cpue.allocate(1,n_cpue,"pred_cpue"); + #ifndef NO_AD_INITIALIZE + pred_cpue.initialize(); + #endif + pred_cope.allocate(1,n_cope,"pred_cope"); + #ifndef NO_AD_INITIALIZE + pred_cope.initialize(); + #endif + Nage_3.allocate(styr,endyr_r,"Nage_3"); + pred_avo.allocate(1,n_avo,"pred_avo"); + #ifndef NO_AD_INITIALIZE + pred_avo.initialize(); + #endif + natage.allocate(styr,endyr_r,1,nages,"natage"); + #ifndef NO_AD_INITIALIZE + natage.initialize(); + #endif + srmod_rec.allocate(styr_est,endyr_est,"srmod_rec"); + #ifndef NO_AD_INITIALIZE + srmod_rec.initialize(); + #endif + Z.allocate(styr,endyr_r,1,nages,"Z"); + #ifndef NO_AD_INITIALIZE + Z.initialize(); + #endif + F.allocate(styr,endyr_r,1,nages,"F"); + #ifndef NO_AD_INITIALIZE + F.initialize(); + #endif + S.allocate(styr,endyr_r,1,nages,"S"); + #ifndef NO_AD_INITIALIZE + S.initialize(); + #endif + M.allocate(styr,endyr_r,1,nages,"M"); + #ifndef NO_AD_INITIALIZE + M.initialize(); + #endif + M_dev.allocate(styr+1,endyr_r,-1.,1.,-8,"M_dev"); + log_sel_fsh.allocate(styr,endyr_r,1,nages,"log_sel_fsh"); + #ifndef NO_AD_INITIALIZE + log_sel_fsh.initialize(); + #endif + sel_fsh.allocate(styr,endyr_r,1,nages,"sel_fsh"); + #ifndef NO_AD_INITIALIZE + sel_fsh.initialize(); + #endif + log_sel_bts.allocate(styr,endyr_r,1,nages,"log_sel_bts"); + #ifndef NO_AD_INITIALIZE + log_sel_bts.initialize(); + #endif + log_sel_ats.allocate(styr,endyr_r,1,nages,"log_sel_ats"); + #ifndef NO_AD_INITIALIZE + log_sel_ats.initialize(); + #endif + ff.allocate("ff"); + #ifndef NO_AD_INITIALIZE + ff.initialize(); + #endif + catch_like.allocate("catch_like"); + #ifndef NO_AD_INITIALIZE + catch_like.initialize(); + #endif + avgsel_fsh.allocate("avgsel_fsh"); + #ifndef NO_AD_INITIALIZE + avgsel_fsh.initialize(); + #endif + avgsel_bts.allocate("avgsel_bts"); + #ifndef NO_AD_INITIALIZE + avgsel_bts.initialize(); + #endif + avgsel_ats.allocate("avgsel_ats"); + #ifndef NO_AD_INITIALIZE + avgsel_ats.initialize(); + #endif + bzero.allocate("bzero"); + #ifndef NO_AD_INITIALIZE + bzero.initialize(); + #endif + surv.allocate("surv"); + #ifndef NO_AD_INITIALIZE + surv.initialize(); + #endif + nthisage.allocate("nthisage"); + #ifndef NO_AD_INITIALIZE + nthisage.initialize(); + #endif + surv_like.allocate(1,3,"surv_like"); + #ifndef NO_AD_INITIALIZE + surv_like.initialize(); + #endif + cpue_like.allocate("cpue_like"); + #ifndef NO_AD_INITIALIZE + cpue_like.initialize(); + #endif + cope_like.allocate("cope_like"); + #ifndef NO_AD_INITIALIZE + cope_like.initialize(); + #endif + avo_like.allocate("avo_like"); + #ifndef NO_AD_INITIALIZE + avo_like.initialize(); + #endif + sel_like.allocate(1,3,"sel_like"); + #ifndef NO_AD_INITIALIZE + sel_like.initialize(); + #endif + sel_like_dev.allocate(1,3,"sel_like_dev"); + #ifndef NO_AD_INITIALIZE + sel_like_dev.initialize(); + #endif + age_like.allocate(1,ngears,"age_like"); + #ifndef NO_AD_INITIALIZE + age_like.initialize(); + #endif + len_like.allocate("len_like"); + #ifndef NO_AD_INITIALIZE + len_like.initialize(); + #endif + wt_like.allocate("wt_like"); + #ifndef NO_AD_INITIALIZE + wt_like.initialize(); + #endif + age_like_offset.allocate(1,ngears,"age_like_offset"); + #ifndef NO_AD_INITIALIZE + age_like_offset.initialize(); + #endif + len_like_offset.allocate("len_like_offset"); + #ifndef NO_AD_INITIALIZE + len_like_offset.initialize(); + #endif + MN_const.allocate("MN_const"); + #ifndef NO_AD_INITIALIZE + MN_const.initialize(); + #endif + MN_const = 1e-3; // Multinomial constant + Priors.allocate(1,4,"Priors"); + #ifndef NO_AD_INITIALIZE + Priors.initialize(); + #endif + rec_like.allocate(1,7,"rec_like"); + #ifndef NO_AD_INITIALIZE + rec_like.initialize(); + #endif + all_like.allocate(1,26,"all_like"); + #ifndef NO_AD_INITIALIZE + all_like.initialize(); + #endif + sumtmp.allocate("sumtmp"); + #ifndef NO_AD_INITIALIZE + sumtmp.initialize(); + #endif + tmpsp.allocate("tmpsp"); + #ifndef NO_AD_INITIALIZE + tmpsp.initialize(); + #endif + log_initage.allocate(2,nages,"log_initage"); + #ifndef NO_AD_INITIALIZE + log_initage.initialize(); + #endif + pred_biom.allocate(styr,endyr_r,"pred_biom"); + #ifndef NO_AD_INITIALIZE + pred_biom.initialize(); + #endif + SRR_SSB.allocate(1,20,"SRR_SSB"); + #ifndef NO_AD_INITIALIZE + SRR_SSB.initialize(); + #endif + fake_SST.allocate(1,40,"fake_SST"); + #ifndef NO_AD_INITIALIZE + fake_SST.initialize(); + #endif + fake_dens.allocate(1,40,"fake_dens"); + #ifndef NO_AD_INITIALIZE + fake_dens.initialize(); + #endif + L1.allocate(10,50,2,"L1"); + L2.allocate(30,90,3,"L2"); + log_alpha.allocate(-1,"log_alpha"); + log_K.allocate(4,"log_K"); + wt_inc.allocate(age_st,age_end-1,"wt_inc"); + #ifndef NO_AD_INITIALIZE + wt_inc.initialize(); + #endif + d_scale.allocate(1,nscale_parm,age_st,age_end,phase_d_scale,"d_scale"); + K.allocate("K"); + #ifndef NO_AD_INITIALIZE + K.initialize(); + #endif + wt_hat.allocate(1,ndat_wt,1,max_nyrs_data,age_st,age_end,"wt_hat"); + #ifndef NO_AD_INITIALIZE + wt_hat.initialize(); + #endif + alphawt.allocate("alphawt"); + #ifndef NO_AD_INITIALIZE + alphawt.initialize(); + #endif + wt_pre.allocate(styr_wt,endyr_wt,age_st,age_end,"wt_pre"); + #ifndef NO_AD_INITIALIZE + wt_pre.initialize(); + #endif + mnwt.allocate(age_st,age_end,"mnwt"); + #ifndef NO_AD_INITIALIZE + mnwt.initialize(); + #endif + coh_eff.allocate(styr_wt-nages_wt-age_st+1,endyr_wt-age_st+3,-15,15,phase_coheff,"coh_eff"); + yr_eff.allocate(styr_wt,endyr_wt+3,-15,15,phase_yreff,"yr_eff"); + wt_last.allocate(age_st,age_end,"wt_last"); + wt_cur.allocate(age_st,age_end,"wt_cur"); + wt_next.allocate(age_st,age_end,"wt_next"); + wt_yraf.allocate(age_st,age_end,"wt_yraf"); + avg_age_msy.allocate("avg_age_msy"); + avgwt_msy.allocate("avgwt_msy"); + MSY.allocate("MSY"); + MSY_wt.allocate("MSY_wt"); + Fmsy.allocate("Fmsy"); + Fmsy_wt.allocate("Fmsy_wt"); + Fmsy2.allocate("Fmsy2"); + Fmsy2_wt.allocate("Fmsy2_wt"); + Fmsy2_dec.allocate(1,10,"Fmsy2_dec"); + Fmsy2_decwt.allocate(1,10,"Fmsy2_decwt"); + Bmsy2.allocate("Bmsy2"); + Bmsy2_wt.allocate("Bmsy2_wt"); + lnFmsy.allocate("lnFmsy"); + lnFmsy2.allocate("lnFmsy2"); + SER_Fmsy.allocate("SER_Fmsy"); + Fendyr_Fmsy.allocate("Fendyr_Fmsy"); + Rmsy.allocate("Rmsy"); + Bmsy.allocate("Bmsy"); + Fcur_Fmsy.allocate(1,nscen,"Fcur_Fmsy"); + Bcur_Bmsy.allocate(1,nscen,"Bcur_Bmsy"); + Bcur_Bmean.allocate(1,nscen,"Bcur_Bmean"); + Bcur3_Bcur.allocate(1,nscen,"Bcur3_Bcur"); + Bcur3_Bmean.allocate(1,nscen,"Bcur3_Bmean"); + Bcur2_Bmsy.allocate(1,nscen,"Bcur2_Bmsy"); + Bcur2_B20.allocate(1,nscen,"Bcur2_B20"); + LTA1_5R.allocate(1,nscen,"LTA1_5R"); + LTA1_5.allocate(1,nscen,"LTA1_5"); + MatAgeDiv1.allocate(1,nscen,"MatAgeDiv1"); + MatAgeDiv2.allocate(1,nscen,"MatAgeDiv2"); + H.allocate(styr,endyr_r,"H"); + #ifndef NO_AD_INITIALIZE + H.initialize(); + #endif + avg_age_mature.allocate(styr,endyr_r,"avg_age_mature"); + #ifndef NO_AD_INITIALIZE + avg_age_mature.initialize(); + #endif + RelEffort.allocate(1,nscen,"RelEffort"); + F40_spb.allocate("F40_spb"); + begbiom.allocate("begbiom"); + DepletionSpawners.allocate("DepletionSpawners"); + SB100.allocate("SB100"); + Current_Spawners.allocate("Current_Spawners"); + pred_rec.allocate(styr,endyr_r,"pred_rec"); + age_3_plus_biom.allocate(styr,endyr_r+2,"age_3_plus_biom"); + ABC_biom.allocate(1,10,"ABC_biom"); + ABC_biom2.allocate(1,10,"ABC_biom2"); + rechat.allocate(1,20,"rechat"); + SRresidhat.allocate(1,40,"SRresidhat"); + #ifndef NO_AD_INITIALIZE + SRresidhat.initialize(); + #endif + SER.allocate(styr,endyr_r,"SER"); + future_SER.allocate(1,nscen,styr_fut,endyr_fut,"future_SER"); + #ifndef NO_AD_INITIALIZE + future_SER.initialize(); + #endif + future_catch.allocate(1,nscen,styr_fut,endyr_fut,"future_catch"); + #ifndef NO_AD_INITIALIZE + future_catch.initialize(); + #endif + future_SSB.allocate(1,nscen,endyr_r,endyr_fut,"future_SSB"); + Age3_Abund.allocate(styr,endyr_r,"Age3_Abund"); + age_1_7_biomass.allocate(styr,endyr_r,"age_1_7_biomass"); + #ifndef NO_AD_INITIALIZE + age_1_7_biomass.initialize(); + #endif + NLL.allocate(1,20,"NLL"); + #ifndef NO_AD_INITIALIZE + NLL.initialize(); + #endif + fff.allocate("fff"); + prior_function_value.allocate("prior_function_value"); + likelihood_function_value.allocate("likelihood_function_value"); + resid_temp_x1.allocate(do_temp_phase,"resid_temp_x1"); + resid_temp_x2.allocate(do_temp_phase,"resid_temp_x2"); + SR_resids_temp.allocate(styr_est,endyr_est,"SR_resids_temp"); + #ifndef NO_AD_INITIALIZE + SR_resids_temp.initialize(); + #endif + SR_resids_like.allocate("SR_resids_like"); + #ifndef NO_AD_INITIALIZE + SR_resids_like.initialize(); + #endif + log_a_II.allocate(1,n_pred_grp,1,n_pred_ages,-12,0,do_pred_phase_ss,"log_a_II"); + log_b_II.allocate(1,n_pred_grp,1,n_pred_ages,0,12,do_pred_phase_ss,"log_b_II"); + log_a_II_vec.allocate(1,n_pred_grp,-12,0,do_pred_phase_ms,"log_a_II_vec"); + log_b_II_vec.allocate(1,n_pred_grp,0,12,do_pred_phase_ms,"log_b_II_vec"); + log_rho.allocate(1,n_pred_grp,1,n_pred_ages,-12,0,do_pred_phase_ms,"log_rho"); + log_resid_M.allocate(1,n_pred_ages,-3,0.1,do_pred_phase,"log_resid_M"); + resid_M.allocate(1,n_pred_ages,"resid_M"); + #ifndef NO_AD_INITIALIZE + resid_M.initialize(); + #endif + resid_M_like.allocate(1,n_pred_ages,"resid_M_like"); + #ifndef NO_AD_INITIALIZE + resid_M_like.initialize(); + #endif + a_II.allocate(1,n_pred_grp,1,n_pred_ages,"a_II"); + #ifndef NO_AD_INITIALIZE + a_II.initialize(); + #endif + b_II.allocate(1,n_pred_grp,1,n_pred_ages,"b_II"); + #ifndef NO_AD_INITIALIZE + b_II.initialize(); + #endif + rho.allocate(1,n_pred_grp,1,n_pred_ages,"rho"); + #ifndef NO_AD_INITIALIZE + rho.initialize(); + #endif + a_II_vec.allocate(1,n_pred_grp,"a_II_vec"); + #ifndef NO_AD_INITIALIZE + a_II_vec.initialize(); + #endif + b_II_vec.allocate(1,n_pred_grp,"b_II_vec"); + #ifndef NO_AD_INITIALIZE + b_II_vec.initialize(); + #endif + natage_strat.allocate(styr,endyr_r,1,n_pred_ages,1,nstrata_pred,"natage_strat"); + #ifndef NO_AD_INITIALIZE + natage_strat.initialize(); + #endif + natage_strat_dens.allocate(styr,endyr_r,1,n_pred_ages,1,nstrata_pred,"natage_strat_dens"); + #ifndef NO_AD_INITIALIZE + natage_strat_dens.initialize(); + #endif + meannatage.allocate(styr,endyr_r,1,nages,"meannatage"); + #ifndef NO_AD_INITIALIZE + meannatage.initialize(); + #endif + meannatage_bystrata.allocate(styr,endyr_r,1,nages,1,nstrata_pred,"meannatage_bystrata"); + #ifndef NO_AD_INITIALIZE + meannatage_bystrata.initialize(); + #endif + mean_dens_bystrata.allocate(styr,endyr_r,1,nages,1,nstrata_pred,"mean_dens_bystrata"); + #ifndef NO_AD_INITIALIZE + mean_dens_bystrata.initialize(); + #endif + mean_dens.allocate(styr,endyr_r,1,n_pred_ages,"mean_dens"); + #ifndef NO_AD_INITIALIZE + mean_dens.initialize(); + #endif + cons.allocate(1,n_pred_grp,styr,endyr_r,1,n_pred_ages,1,nstrata_pred,"cons"); + #ifndef NO_AD_INITIALIZE + cons.initialize(); + #endif + natmort_pred.allocate(1,n_pred_ages,styr,endyr_r,1,nstrata_pred,"natmort_pred"); + #ifndef NO_AD_INITIALIZE + natmort_pred.initialize(); + #endif + M_pred.allocate(styr,endyr_r,1,n_pred_ages,1,nstrata_pred,1,n_pred_grp_nonpoll,"M_pred"); + #ifndef NO_AD_INITIALIZE + M_pred.initialize(); + #endif + M_pred_tmp.allocate(1,n_pred_grp_nonpoll,1,n_pred_ages,"M_pred_tmp"); + #ifndef NO_AD_INITIALIZE + M_pred_tmp.initialize(); + #endif + M_pred_sum.allocate(styr,endyr_r,1,n_pred_ages,1,nstrata_pred,"M_pred_sum"); + #ifndef NO_AD_INITIALIZE + M_pred_sum.initialize(); + #endif + Z_pred.allocate(styr,endyr_r,1,n_pred_ages,1,nstrata_pred,"Z_pred"); + #ifndef NO_AD_INITIALIZE + Z_pred.initialize(); + #endif + S_pred.allocate(styr,endyr_r,1,n_pred_ages,1,nstrata_pred,"S_pred"); + #ifndef NO_AD_INITIALIZE + S_pred.initialize(); + #endif + M_pred_avg.allocate(1,n_pred_ages,styr,endyr_r,"M_pred_avg"); + #ifndef NO_AD_INITIALIZE + M_pred_avg.initialize(); + #endif + cons_atage.allocate(1,n_pred_grp,styr,endyr_r,1,n_pred_ages,"cons_atage"); + #ifndef NO_AD_INITIALIZE + cons_atage.initialize(); + #endif + cons_atage_wt.allocate(1,n_pred_grp,styr,endyr_r,1,n_pred_ages,"cons_atage_wt"); + #ifndef NO_AD_INITIALIZE + cons_atage_wt.initialize(); + #endif + pred_cons.allocate(1,n_pred_grp,styr,endyr_r,"pred_cons"); + #ifndef NO_AD_INITIALIZE + pred_cons.initialize(); + #endif + eac_cons.allocate(1,n_pred_grp,styr,endyr_r,1,n_pred_ages,"eac_cons"); + #ifndef NO_AD_INITIALIZE + eac_cons.initialize(); + #endif + ssq_cons.allocate(1,n_pred_grp,"ssq_cons"); + #ifndef NO_AD_INITIALIZE + ssq_cons.initialize(); + #endif + oac_cons_like_offset.allocate(1,n_pred_grp,"oac_cons_like_offset"); + #ifndef NO_AD_INITIALIZE + oac_cons_like_offset.initialize(); + #endif + age_like_cons.allocate(1,n_pred_grp,"age_like_cons"); + #ifndef NO_AD_INITIALIZE + age_like_cons.initialize(); + #endif + pred_cpup.allocate(1,n_pred_grp,styr,endyr_r,1,n_pred_ages,"pred_cpup"); + #ifndef NO_AD_INITIALIZE + pred_cpup.initialize(); + #endif + implied_cpuppa.allocate(1,n_pred_grp,styr,endyr_r,1,n_pred_ages,1,nstrata_pred,"implied_cpuppa"); + #ifndef NO_AD_INITIALIZE + implied_cpuppa.initialize(); + #endif + implied_obs_cons_bystrata.allocate(1,n_pred_grp,styr,endyr_r,1,n_pred_ages,1,nstrata_pred,"implied_obs_cons_bystrata"); + #ifndef NO_AD_INITIALIZE + implied_obs_cons_bystrata.initialize(); + #endif + implied_prop_Cmax.allocate(1,n_pred_grp,styr,endyr_r,1,n_pred_ages,1,nstrata_pred,"implied_prop_Cmax"); + #ifndef NO_AD_INITIALIZE + implied_prop_Cmax.initialize(); + #endif + max_F_yldcrv.allocate("max_F_yldcrv"); + #ifndef NO_AD_INITIALIZE + max_F_yldcrv.initialize(); + #endif + F_yldcrv.allocate(1,40,"F_yldcrv"); + #ifndef NO_AD_INITIALIZE + F_yldcrv.initialize(); + #endif + yield_curve.allocate(1,40,"yield_curve"); + #ifndef NO_AD_INITIALIZE + yield_curve.initialize(); + #endif + natmort_fut.allocate(1,nages,"natmort_fut"); + #ifndef NO_AD_INITIALIZE + natmort_fut.initialize(); + #endif + compweightsnew.allocate(1,n_pred_grp_nonpoll,"compweightsnew"); + #ifndef NO_AD_INITIALIZE + compweightsnew.initialize(); + #endif + comp_mcian_wgt_inv.allocate(1,n_pred_grp_nonpoll,1,nyrs_cons_nonpoll,"comp_mcian_wgt_inv"); + #ifndef NO_AD_INITIALIZE + comp_mcian_wgt_inv.initialize(); + #endif + comp_mcian_wgt.allocate(1,n_pred_grp_nonpoll,1,nyrs_cons_nonpoll,"comp_mcian_wgt"); + #ifndef NO_AD_INITIALIZE + comp_mcian_wgt.initialize(); + #endif + consweightsnew.allocate(1,n_pred_grp_nonpoll,"consweightsnew"); + #ifndef NO_AD_INITIALIZE + consweightsnew.initialize(); + #endif + cons_nr.allocate(1,n_pred_grp_nonpoll,1,nyrs_cons_nonpoll,"cons_nr"); + #ifndef NO_AD_INITIALIZE + cons_nr.initialize(); + #endif + comp_nr.allocate(1,n_pred_grp_nonpoll,1,comp_nr_ub,"comp_nr"); + #ifndef NO_AD_INITIALIZE + comp_nr.initialize(); + #endif + pred_rec_alpha.allocate("pred_rec_alpha"); + #ifndef NO_AD_INITIALIZE + pred_rec_alpha.initialize(); + #endif + srmod_rec_alpha.allocate("srmod_rec_alpha"); + #ifndef NO_AD_INITIALIZE + srmod_rec_alpha.initialize(); + #endif +} + +void model_parameters::preliminary_calculations(void) +{ + +#if defined(USE_ADPVM) + + admaster_slave_variable_interface(*this); + +#endif + // 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 +#ifdef DEBUG + std::cout << "DEBUG: Total gradient stack used is " << gradient_structure::get()->GRAD_STACK1->total() << " out of " << gradient_structure::get_GRADSTACK_BUFFER_SIZE() << std::endl;; + std::cout << "DEBUG: Total dvariable address used is " << gradient_structure::get()->GRAD_LIST->total_addresses() << " out of " << gradient_structure::get_MAX_DLINKS() << std::endl;; + std::cout << "DEBUG: Total dvariable address used is " << gradient_structure::get()->ARR_LIST1->get_max_last_offset() << " out of " << gradient_structure::get_ARRAY_MEMBLOCK_SIZE() << std::endl;; +#endif +} + +void model_parameters::update_compweights(void) +{ + ofstream& srecout= *pad_srecout; + ofstream& projout= *pad_projout; + ofstream& nofish= *pad_nofish; + ofstream& projout2= *pad_projout2; + ofstream& eval= *pad_eval; + random_number_generator& rng= *pad_rng; + 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; + } +} + +void model_parameters::Get_Selectivity(void) +{ + ofstream& srecout= *pad_srecout; + ofstream& projout= *pad_projout; + ofstream& nofish= *pad_nofish; + ofstream& projout2= *pad_projout2; + ofstream& eval= *pad_eval; + random_number_generator& rng= *pad_rng; + 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; +} + +void model_parameters::GetDependentVar(void) +{ + ofstream& srecout= *pad_srecout; + ofstream& projout= *pad_projout; + ofstream& nofish= *pad_nofish; + ofstream& projout2= *pad_projout2; + ofstream& eval= *pad_eval; + random_number_generator& rng= *pad_rng; + // 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 +} + +void model_parameters::compute_spr_rates(void) +{ + ofstream& srecout= *pad_srecout; + ofstream& projout= *pad_projout; + ofstream& nofish= *pad_nofish; + ofstream& projout2= *pad_projout2; + ofstream& eval= *pad_eval; + random_number_generator& rng= *pad_rng; + //Compute SPR Rates + F35 = get_spr_rates(.35); + F40 = get_spr_rates(.40); +} + +dvariable model_parameters::get_spr_rates(double spr_percent,dvar_vector sel) +{ + ofstream& srecout= *pad_srecout; + ofstream& projout= *pad_projout; + ofstream& nofish= *pad_nofish; + ofstream& projout2= *pad_projout2; + ofstream& eval= *pad_eval; + random_number_generator& rng= *pad_rng; + 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); +} + +dvariable model_parameters::spr_ratio(dvariable trial_F,dvar_vector& sel) +{ + ofstream& srecout= *pad_srecout; + ofstream& projout= *pad_projout; + ofstream& nofish= *pad_nofish; + ofstream& projout2= *pad_projout2; + ofstream& eval= *pad_eval; + random_number_generator& rng= *pad_rng; + 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(); +} + +void model_parameters::Get_Cons_at_Age(void) +{ + ofstream& srecout= *pad_srecout; + ofstream& projout= *pad_projout; + ofstream& nofish= *pad_nofish; + ofstream& projout2= *pad_projout2; + ofstream& eval= *pad_eval; + random_number_generator& rng= *pad_rng; + // 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))))); + } + } + */ + // +===+====+==+==+==+==+==+==+==+====+====+==+==+===+====+==+==+==+==+==+==+==+====+====+====+ +} + +void model_parameters::Evaluate_Objective_Function(void) +{ + ofstream& srecout= *pad_srecout; + ofstream& projout= *pad_projout; + ofstream& nofish= *pad_nofish; + ofstream& projout2= *pad_projout2; + ofstream& eval= *pad_eval; + random_number_generator& rng= *pad_rng; + // 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); + // 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 ); + // Covariance on observed population (numbers) switch + if (DoCovBTS && current_phase()>4) + { + surv_like(1) = .5 * srv_tmp * inv_bts_cov * srv_tmp; + // cout <<"Survey likelihood: " << surv_like(1) << endl; + } + 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); +} + +void model_parameters::Robust_Likelihood(void) +{ + ofstream& srecout= *pad_srecout; + ofstream& projout= *pad_projout; + ofstream& nofish= *pad_nofish; + ofstream& projout2= *pad_projout2; + ofstream& eval= *pad_eval; + random_number_generator& rng= *pad_rng; + 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); +} + +void model_parameters::Multinomial_Likelihood(void) +{ + ofstream& srecout= *pad_srecout; + ofstream& projout= *pad_projout; + ofstream& nofish= *pad_nofish; + ofstream& projout2= *pad_projout2; + ofstream& eval= *pad_eval; + random_number_generator& rng= *pad_rng; + 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; + default: + 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); + } + //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) +} + +dvariable model_parameters::robust_p(dmatrix& obs,dvar_matrix& pred,const dvariable& a, const data_ivector& b, const int amin, const int amax) +{ + ofstream& srecout= *pad_srecout; + ofstream& projout= *pad_projout; + ofstream& nofish= *pad_nofish; + ofstream& projout2= *pad_projout2; + ofstream& eval= *pad_eval; + random_number_generator& rng= *pad_rng; + 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); +} + +dvariable model_parameters::robust_p(const dmatrix& obs,const dvar_matrix& pred,const dvariable& a, const data_ivector& b) +{ + ofstream& srecout= *pad_srecout; + ofstream& projout= *pad_projout; + ofstream& nofish= *pad_nofish; + ofstream& projout2= *pad_projout2; + ofstream& eval= *pad_eval; + random_number_generator& rng= *pad_rng; + 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); +} + +dvariable model_parameters::robust_p(const dmatrix& obs, const dvar_matrix& pred, const dvariable& a, const dvariable& b) +{ + ofstream& srecout= *pad_srecout; + ofstream& projout= *pad_projout; + ofstream& nofish= *pad_nofish; + ofstream& projout2= *pad_projout2; + ofstream& eval= *pad_eval; + random_number_generator& rng= *pad_rng; + 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); +} + +dvariable model_parameters::robust_p(const dvector& obs, const dvar_vector& pred, const dvariable& a, const dvariable& b) +{ + ofstream& srecout= *pad_srecout; + ofstream& projout= *pad_projout; + ofstream& nofish= *pad_nofish; + ofstream& projout2= *pad_projout2; + ofstream& eval= *pad_eval; + random_number_generator& rng= *pad_rng; + 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); +} + +void model_parameters::write_srec(void) +{ + ofstream& srecout= *pad_srecout; + ofstream& projout= *pad_projout; + ofstream& nofish= *pad_nofish; + ofstream& projout2= *pad_projout2; + ofstream& eval= *pad_eval; + random_number_generator& rng= *pad_rng; + 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 << " " <int(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 + } +} + +double model_parameters::sdnr(const dvector& obs, const dvar_vector& pred, const dvector& sig) +{ + ofstream& srecout= *pad_srecout; + ofstream& projout= *pad_projout; + ofstream& nofish= *pad_nofish; + ofstream& projout2= *pad_projout2; + ofstream& eval= *pad_eval; + random_number_generator& rng= *pad_rng; + 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; + } + **/ +} + +dvector model_parameters::rmultinomial(const dvector ptmp, const int n_sam) +{ + ofstream& srecout= *pad_srecout; + ofstream& projout= *pad_projout; + ofstream& nofish= *pad_nofish; + ofstream& projout2= *pad_projout2; + ofstream& eval= *pad_eval; + random_number_generator& rng= *pad_rng; + 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. + **/ +} + +double model_parameters::calc_Francis_weights(const dmatrix oac, const dvar_matrix eac, const ivector sam ) +{ + ofstream& srecout= *pad_srecout; + ofstream& projout= *pad_projout; + ofstream& nofish= *pad_nofish; + ofstream& projout2= *pad_projout2; + ofstream& eval= *pad_eval; + random_number_generator& rng= *pad_rng; + { + 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"<(std::chrono::high_resolution_clock::now() - start).count() << " microseconds." << endl; + #ifndef __SUNPRO_C +bool failedtest = false; +if (std::fetestexcept(FE_DIVBYZERO)) + { cerr << "Error: Detected division by zero." << endl; failedtest = true; } +if (std::fetestexcept(FE_INVALID)) + { cerr << "Error: Detected invalid argument." << endl; failedtest = true; } +if (std::fetestexcept(FE_OVERFLOW)) + { cerr << "Error: Detected overflow." << endl; failedtest = true; } +if (std::fetestexcept(FE_UNDERFLOW)) + { cerr << "Error: Detected underflow." << endl; } +if (failedtest) { std::abort(); } + #endif +#endif + return 0; +} + +extern "C" { + void ad_boundf(int i) + { + /* so we can stop here */ + exit(i); + } +} diff --git a/src/2022/pm.tpl b/src/2022/pm.tpl new file mode 100644 index 0000000..ec38c28 --- /dev/null +++ b/src/2022/pm.tpl @@ -0,0 +1,6514 @@ +///////////////////////////////////////////////////////////// +// 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) >> 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(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 avo_sel(1,15) + !! avo_sel(1)=0.0; avo_sel(2)=1; avo_sel(3)=1; avo_sel(4)=0.85; avo_sel(5)=0.7; avo_sel(6)=0.55; avo_sel(7)=0.3; avo_sel(8)=0.15; avo_sel(9)=0.05; avo_sel(10)=0.01; avo_sel(11)=0.01; avo_sel(12)=0.01; avo_sel(13)=0.01; avo_sel(14)=0.01; avo_sel(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)); + 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); + + // ------------------------------------------------------ + // read in wt-age data by + !! ad_comm::change_datafile_name(Wtage_file); + // These were estimated in an RE model... + init_number log_sd_coh + init_number log_sd_yr + init_number cur_yr + init_int styr_wt + init_int endyr_wt + init_int ndat_wt + // Number of years of observations in each data set + init_ivector nyrs_data(1,ndat_wt); + // Actual years of observations in each data set (ragged array) + init_imatrix yrs_data(1,ndat_wt,1,nyrs_data); + LOCAL_CALCS + if (endyr_r < endyr) + for (int h=1;h<=ndat_wt;h++) + { + int itmp = 1; + for (i=styr;i<=endyr_r;i++) + { + if (i == yrs_data(h,itmp) ) + itmp++; + } + nyrs_data(h)= itmp-1; + } + END_CALCS + init_int age_st + init_int age_end + int nages_wt; + !! nages_wt = age_end - age_st + 1; + int nscale_parm; + !! nscale_parm = ndat_wt-1; + // init_matrix wt_obs(styr_wt,endyr_wt,age_st,age_end) + // init_matrix sd_obs(styr_wt,endyr_wt,age_st,age_end) + init_3darray wt_obs(1,ndat_wt,1,nyrs_data,age_st,age_end); + init_3darray sd_obs(1,ndat_wt,1,nyrs_data,age_st,age_end); + int phase_d_scale; + !! if (ndat_wt>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); + 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); + // 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 ); + + // Covariance on observed population (numbers) switch + if (DoCovBTS && current_phase()>4) + { + surv_like(1) = .5 * srv_tmp * inv_bts_cov * srv_tmp; + // cout <<"Survey likelihood: " << surv_like(1) << endl; + } + 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; + default: + 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); + } + //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 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"< + #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 " " <