#here I am going to make sure I am fishing at F40% #I will remove stochasticity #and fix the NE calcs library(batch) parseCommandArgs() library(Rcpp) library(inline) require(RcppArmadillo) runIBD = ' using namespace Rcpp; using namespace sugar; RNGScope scope; // Variables (and constants) section NumericVector inputs(INPUTS); int SSBinit; SSBinit=inputs[0]; int Nyrs; Nyrs=inputs[1]; int nsamp; nsamp=inputs[2]; int ngen; ngen=inputs[3]; int nmig; nmig=inputs[4]; double fishmort; fishmort=inputs[5]; int nsims; nsims=inputs[6]; int mmt; mmt=inputs[7]; int N; N=inputs[8]; int Yrs; Yrs=inputs[9]; double dpow; dpow=inputs[10]; double nages; nages=inputs[11]; double CVa; CVa=inputs[12]; double CVs; CVs=inputs[13]; double sigma_jr2=inputs[14]; double h=inputs[15]; Rcout<<"test"<(rnorm(Nyrs,0,pow(sigma_jr2,0.5))); // initial numbers section (if there was one recruit) N_til_init_a[0]=1; for (int i=1; i < 20; i++){ //fill in N_til_init; N_til_init_a[i]=N_til_init_a[i-1]*exp(-M);} N_til_init_a[20]=N_til_init_a[19]*exp(-M)/(1-exp(-M)); //establish the SBPR with 1 recruit and F40% N_til_init_a40[0]=1; for (int i=1; i < 20; i++){ //fill in N_til_init; N_til_init_a40[i]=N_til_init_a40[i-1]*exp(-(M+fishsel[i-1]*fishmort));} N_til_init_a40[20]=N_til_init_a40[19]*exp(-(M+fishsel[19]*fishmort))/(1-exp(-(M+fishsel[19]*fishmort))); // establish Q (Maturity at age) [from 1 to 21 - 0 to 20 actually] for (int i=0; i < 21; i++){ Q[i]=1.0/(1+exp(-(A+age[i]*B)));} Q[0]=0;Q[1]=0; // establish length, weight of males and females for (int i=0; i < 21; i++){ L_fem[i]=Linf_fem*(1-exp(-K_fem*(age[i]-A0_fem))); Wt_fem[i]=psi*pow((Linf_fem*(1-exp(-K_fem*(age[i]-A0_fem)))),theta); //0.001 is for kg to mt conversion Wt_mal[i]=psi*pow((Linf_mal*(1-exp(-K_mal*(age[i]-A0_mal)))),theta);} NumericVector tempvec1 = Wt_fem*N_til_init_a*Q; double den=std::accumulate(tempvec1.begin(), tempvec1.end(), 0.0); NumericVector Wt=(Wt_fem+Wt_mal)/2; // determine initial number of recruits and the SSB in each population (pops 1 and 2) S_l0=SSBinit*0.10; //S_l0 initial spawning stock biomass for pop1 R_l0A[0]=2.0*S_l0/den; //initial number recruits pop1 NOTE: if only one pop, all are in popA R_l0B[0]=R_l0A[0]; R_l0C[0]=R_l0A[0];R_l0D[0]=R_l0A[0]; R_l0E[0]=R_l0A[0]; R_l0F[0]=R_l0A[0];R_l0G[0]=R_l0A[0]; R_l0H[0]=R_l0A[0]; R_l0I[0]=R_l0A[0];R_l0J[0]=R_l0A[0]; double R_l0A_init=R_l0A[0]; //record the first recruits for beginning of mgmt_counter loop double R_l0B_init=R_l0B[0]; double R_l0C_init=R_l0C[0]; double R_l0D_init=R_l0D[0]; double R_l0E_init=R_l0E[0]; double R_l0F_init=R_l0F[0]; double R_l0G_init=R_l0G[0]; double R_l0H_init=R_l0H[0]; double R_l0I_init=R_l0I[0]; double R_l0J_init=R_l0J[0]; //n_init_a1 is the number of fish in each age class //this loop will also round to an even number in each age class so you can split into males and females for (int i=0; i<21; i++){ N_init_a1[i]=N_til_init_a[i]*R_l0A[0]; temp=int(N_init_a1[i]); if((N_init_a1[i]+0.5)>=(int(N_init_a1[i])+1)){ temp=int(N_init_a1[i])+1;} N_init_a1[i]=temp; if (temp%2!=0){ temp2=as(runif(1,-1,1)); if(temp2>=0){N_init_a1[i]=N_init_a1[i]+1;} else{N_init_a1[i]=N_init_a1[i]-1;} } N_init_a2[i]=N_til_init_a[i]*R_l0B[0]; temp=int(N_init_a2[i]); if((N_init_a2[i]+0.5)>=(int(N_init_a2[i])+1)){ temp=int(N_init_a2[i])+1;} N_init_a2[i]=temp; if (temp%2!=0){ temp2=as(runif(1,-1,1)); if(temp2>=0){N_init_a2[i]=N_init_a2[i]+1;} else{N_init_a2[i]=N_init_a2[i]-1;} } N_init_a3[i]=N_til_init_a[i]*R_l0C[0]; temp=int(N_init_a3[i]); if((N_init_a3[i]+0.5)>=(int(N_init_a3[i])+1)){ temp=int(N_init_a3[i])+1;} N_init_a3[i]=temp; if (temp%2!=0){ temp2=as(runif(1,-1,1)); if(temp2>=0){N_init_a3[i]=N_init_a3[i]+1;} else{N_init_a3[i]=N_init_a3[i]-1;} } N_init_a4[i]=N_til_init_a[i]*R_l0D[0]; temp=int(N_init_a4[i]); if((N_init_a4[i]+0.5)>=(int(N_init_a4[i])+1)){ temp=int(N_init_a4[i])+1;} N_init_a4[i]=temp; if (temp%2!=0){ temp2=as(runif(1,-1,1)); if(temp2>=0){N_init_a4[i]=N_init_a4[i]+1;} else{N_init_a4[i]=N_init_a4[i]-1;} } N_init_a5[i]=N_til_init_a[i]*R_l0E[0]; temp=int(N_init_a5[i]); if((N_init_a5[i]+0.5)>=(int(N_init_a5[i])+1)){ temp=int(N_init_a5[i])+1;} N_init_a5[i]=temp; if (temp%2!=0){ temp2=as(runif(1,-1,1)); if(temp2>=0){N_init_a5[i]=N_init_a5[i]+1;} else{N_init_a5[i]=N_init_a5[i]-1;} } N_init_a6[i]=N_til_init_a[i]*R_l0F[0]; temp=int(N_init_a6[i]); if((N_init_a6[i]+0.5)>=(int(N_init_a6[i])+1)){ temp=int(N_init_a6[i])+1;} N_init_a6[i]=temp; if (temp%2!=0){ temp2=as(runif(1,-1,1)); if(temp2>=0){N_init_a6[i]=N_init_a6[i]+1;} else{N_init_a6[i]=N_init_a6[i]-1;} } N_init_a7[i]=N_til_init_a[i]*R_l0G[0]; temp=int(N_init_a7[i]); if((N_init_a7[i]+0.5)>=(int(N_init_a7[i])+1)){ temp=int(N_init_a7[i])+1;} N_init_a7[i]=temp; if (temp%2!=0){ temp2=as(runif(1,-1,1)); if(temp2>=0){N_init_a7[i]=N_init_a7[i]+1;} else{N_init_a7[i]=N_init_a7[i]-1;} } N_init_a8[i]=N_til_init_a[i]*R_l0H[0]; temp=int(N_init_a8[i]); if((N_init_a8[i]+0.5)>=(int(N_init_a8[i])+1)){ temp=int(N_init_a8[i])+1;} N_init_a8[i]=temp; if (temp%2!=0){ temp2=as(runif(1,-1,1)); if(temp2>=0){N_init_a8[i]=N_init_a8[i]+1;} else{N_init_a8[i]=N_init_a8[i]-1;} } N_init_a9[i]=N_til_init_a[i]*R_l0I[0]; temp=int(N_init_a7[i]); if((N_init_a9[i]+0.5)>=(int(N_init_a9[i])+1)){ temp=int(N_init_a9[i])+1;} N_init_a9[i]=temp; if (temp%2!=0){ temp2=as(runif(1,-1,1)); if(temp2>=0){N_init_a9[i]=N_init_a9[i]+1;} else{N_init_a9[i]=N_init_a9[i]-1;} } N_init_a10[i]=N_til_init_a[i]*R_l0J[0]; temp=int(N_init_a10[i]); if((N_init_a10[i]+0.5)>=(int(N_init_a10[i])+1)){ temp=int(N_init_a10[i])+1;} N_init_a10[i]=temp; if (temp%2!=0){ temp2=as(runif(1,-1,1)); if(temp2>=0){N_init_a10[i]=N_init_a10[i]+1;} else{N_init_a10[i]=N_init_a10[i]-1;} } }//end this for loop N_init_a1_FIRST = clone(N_init_a1); N_init_a2_FIRST = clone(N_init_a2); N_init_a3_FIRST = clone(N_init_a3); N_init_a4_FIRST = clone(N_init_a4); N_init_a5_FIRST = clone(N_init_a5); N_init_a6_FIRST = clone(N_init_a6); N_init_a7_FIRST = clone(N_init_a7); N_init_a8_FIRST = clone(N_init_a8); N_init_a9_FIRST = clone(N_init_a9); N_init_a10_FIRST = clone(N_init_a10); N_init_a = N_init_a1+N_init_a2+N_init_a3+N_init_a4+N_init_a5+N_init_a6+N_init_a7+N_init_a8+N_init_a9+N_init_a10; nums=N_init_a; NumericVector tempvec = Wt_fem*N_init_a*Q; S_lyINIT=0.5*std::accumulate(tempvec.begin(), tempvec.end(), 0.0);//Initial Spawning bimass SPR_init=S_lyINIT/N_init_a[0]; //initial SPBR tempvec = Wt_fem*N_til_init_a40*Q; SPR_init_F40=0.5*std::accumulate(tempvec.begin(), tempvec.end(), 0.0);//Initial Spawning bimass per 1 recruit pop1_size = std::accumulate(N_init_a1.begin(),N_init_a1.end(),0.0); pop2_size = std::accumulate(N_init_a2.begin(),N_init_a2.end(),0.0); pop3_size = std::accumulate(N_init_a3.begin(),N_init_a3.end(),0.0); pop4_size = std::accumulate(N_init_a4.begin(),N_init_a4.end(),0.0); pop5_size = std::accumulate(N_init_a5.begin(),N_init_a5.end(),0.0); pop6_size = std::accumulate(N_init_a6.begin(),N_init_a6.end(),0.0); pop7_size = std::accumulate(N_init_a7.begin(),N_init_a7.end(),0.0); pop8_size = std::accumulate(N_init_a8.begin(),N_init_a8.end(),0.0); pop9_size = std::accumulate(N_init_a9.begin(),N_init_a9.end(),0.0); pop10_size = std::accumulate(N_init_a10.begin(),N_init_a10.end(),0.0); int N_init=pop1_size+pop2_size+pop3_size+pop4_size+pop5_size+pop6_size+pop7_size+pop8_size+pop9_size+pop10_size; NumericVector temp6=(N_init_a1*Wt_fem*Q); double S_lyINITA=0.5*std::accumulate(temp6.begin(),temp6.end(),0.0); NumericVector temp7=(N_init_a2*Wt_fem*Q); double S_lyINITB=0.5*std::accumulate(temp7.begin(),temp7.end(),0.0); NumericVector temp8=(N_init_a3*Wt_fem*Q); double S_lyINITC=0.5*std::accumulate(temp8.begin(),temp8.end(),0.0); NumericVector temp9=(N_init_a4*Wt_fem*Q); double S_lyINITD=0.5*std::accumulate(temp9.begin(),temp9.end(),0.0); NumericVector temp10=(N_init_a5*Wt_fem*Q); double S_lyINITE=0.5*std::accumulate(temp10.begin(),temp10.end(),0.0); NumericVector temp11=(N_init_a6*Wt_fem*Q); double S_lyINITF=0.5*std::accumulate(temp11.begin(),temp11.end(),0.0); NumericVector temp12=(N_init_a7*Wt_fem*Q); double S_lyINITG=0.5*std::accumulate(temp12.begin(),temp12.end(),0.0); NumericVector temp13=(N_init_a8*Wt_fem*Q); double S_lyINITH=0.5*std::accumulate(temp13.begin(),temp13.end(),0.0); NumericVector temp14=(N_init_a9*Wt_fem*Q); double S_lyINITI=0.5*std::accumulate(temp14.begin(),temp14.end(),0.0); NumericVector temp15=(N_init_a10*Wt_fem*Q); double S_lyINITJ=0.5*std::accumulate(temp15.begin(),temp15.end(),0.0); double SPR_initA; //initial SBPR popA double SPR_initB; //initial SBPR popB double SPR_initC; //initial SBPR popC double SPR_initD; //initial SBPR popD double SPR_initE; //initial SBPR popE double SPR_initF; //initial SBPR popF double SPR_initG; //initial SBPR popG double SPR_initH; //initial SBPR popH double SPR_initI; //initial SBPR popI double SPR_initJ; //initial SBPR popJ if(N_init_a1[0]==0){SPR_initA=0;} //have not yet reversed the order so zero is the young recruits else{SPR_initA=S_lyINITA/N_init_a1[0];} if(N_init_a2[0]==0){SPR_initB=0;} else{SPR_initB=S_lyINITB/N_init_a2[0];} if(N_init_a3[0]==0){SPR_initC=0;} //have not yet reversed the order so zero is the young recruits else{SPR_initC=S_lyINITC/N_init_a3[0];} if(N_init_a4[0]==0){SPR_initD=0;} else{SPR_initD=S_lyINITD/N_init_a4[0];} if(N_init_a5[0]==0){SPR_initE=0;} //have not yet reversed the order so zero is the young recruits else{SPR_initE=S_lyINITE/N_init_a5[0];} if(N_init_a6[0]==0){SPR_initF=0;} else{SPR_initF=S_lyINITF/N_init_a6[0];} if(N_init_a7[0]==0){SPR_initG=0;} //have not yet reversed the order so zero is the young recruits else{SPR_initG=S_lyINITG/N_init_a7[0];} if(N_init_a8[0]==0){SPR_initH=0;} else{SPR_initH=S_lyINITH/N_init_a8[0];} if(N_init_a9[0]==0){SPR_initI=0;} //have not yet reversed the order so zero is the young recruits else{SPR_initI=S_lyINITI/N_init_a9[0];} if(N_init_a10[0]==0){SPR_initJ=0;} else{SPR_initJ=S_lyINITJ/N_init_a10[0];} NumericVector TACA(Nyrs); NumericVector TACB(Nyrs); NumericVector TACC(Nyrs); NumericVector TACD(Nyrs); NumericVector TACE(Nyrs); NumericVector TACF(Nyrs); NumericVector TACG(Nyrs); NumericVector TACH(Nyrs); NumericVector TACI(Nyrs); NumericVector TACJ(Nyrs); NumericVector TAC(Nyrs); NumericVector TAC1(Nyrs); //for area 1 NumericVector TAC2(Nyrs); //for area 2 NumericVector SSB_hat(Nyrs); NumericVector SSB_hat_S(Nyrs); NumericVector Fishvec(Nyrs); //vector of fishing mortality for combined management NumericVector SSB_hatA(Nyrs); NumericVector SSB_hatB(Nyrs); NumericVector SSB_hatC(Nyrs); NumericVector SSB_hatD(Nyrs); NumericVector SSB_hatE(Nyrs); NumericVector SSB_hatF(Nyrs); NumericVector SSB_hatG(Nyrs); NumericVector SSB_hatH(Nyrs); NumericVector SSB_hatI(Nyrs); NumericVector SSB_hatJ(Nyrs); NumericVector SSB_hatA_S(Nyrs); NumericVector SSB_hatB_S(Nyrs); NumericVector SSB_hatC_S(Nyrs); NumericVector SSB_hatD_S(Nyrs); NumericVector SSB_hatE_S(Nyrs); NumericVector SSB_hatF_S(Nyrs); NumericVector SSB_hatG_S(Nyrs); NumericVector SSB_hatH_S(Nyrs); NumericVector SSB_hatI_S(Nyrs); NumericVector SSB_hatJ_S(Nyrs); NumericVector SPRAtrue(Nyrs); NumericVector SPRBtrue(Nyrs); NumericVector SPRCtrue(Nyrs); NumericVector SPRDtrue(Nyrs); NumericVector SPREtrue(Nyrs); NumericVector SPRFtrue(Nyrs); NumericVector SPRGtrue(Nyrs); NumericVector SPRHtrue(Nyrs); NumericVector SPRItrue(Nyrs); NumericVector SPRJtrue(Nyrs); NumericVector epsilonA(Nyrs); NumericVector epsilonB(Nyrs); NumericVector epsilonC(Nyrs); NumericVector epsilonD(Nyrs); NumericVector epsilonE(Nyrs); NumericVector epsilonF(Nyrs); NumericVector epsilonG(Nyrs); NumericVector epsilonH(Nyrs); NumericVector epsilonI(Nyrs); NumericVector epsilonJ(Nyrs); NumericVector epsilon(Nyrs); epsilonA[0]=0; epsilonB[0]=0; epsilonC[0]=0;epsilonD[0]=0;epsilonE[0]=0;epsilonF[0]=0;epsilonG[0]=0;epsilonH[0]=0;epsilonI[0]=0;epsilonJ[0]=0; epsilon[0]=0; NumericVector N_hat_aA(21); NumericVector N_hat_aB(21); NumericVector N_hat_aC(21); NumericVector N_hat_aD(21); NumericVector N_hat_aE(21); NumericVector N_hat_aF(21); NumericVector N_hat_aG(21); NumericVector N_hat_aH(21); NumericVector N_hat_aI(21); NumericVector N_hat_aJ(21); NumericVector N_hat_aA_S(21); NumericVector N_hat_aB_S(21); NumericVector N_hat_aC_S(21); NumericVector N_hat_aD_S(21); NumericVector N_hat_aE_S(21); NumericVector N_hat_aF_S(21); NumericVector N_hat_aG_S(21); NumericVector N_hat_aH_S(21); NumericVector N_hat_aI_S(21); NumericVector N_hat_aJ_S(21); NumericVector N_hat_a(21); double sigma_ja2A; double sigma_ja2B; double sigma_ja2C; double sigma_ja2D; double sigma_ja2E; double sigma_ja2F; double sigma_ja2G; double sigma_ja2H; double sigma_ja2I; double sigma_ja2J; double sigma_ja2_tot; double sigma_js2A; double sigma_js2B; double sigma_js2C; double sigma_js2D; double sigma_js2E; double sigma_js2F; double sigma_js2G; double sigma_js2H; double sigma_js2I; double sigma_js2J; double sigma_js2_tot; double sigma_jr2_tot; double psiA; double psiB; double psiC;double psiD;double psiE;double psiF;double psiG;double psiH;double psiI;double psiJ; double psi_star; double rho=0.707; double Fish=0; //set here so initially it is not a problem NumericMatrix Effort_mat = NumericMatrix(Nyrs,10); double S_lyA; double S_lyB; double S_lyC;double S_lyD;double S_lyE;double S_lyF;double S_lyG;double S_lyH;double S_lyI;double S_lyJ; double BioA; double BioB; double BioC; double BioD; double BioE; double BioF; double BioG; double BioH; double BioI; double BioJ; double eta_jyA; double eta_jyB; double eta_jyC;double eta_jyD;double eta_jyE;double eta_jyF;double eta_jyG;double eta_jyH;double eta_jyI;double eta_jyJ; double eta_jy; NumericVector temp16; NumericVector temp16a;NumericVector temp16b;NumericVector temp16c;NumericVector temp16d;NumericVector temp16e; NumericVector temp16f;NumericVector temp16g;NumericVector temp16h;NumericVector temp16i;NumericVector temp16j; NumericVector temp17; int temp18; int temp19; double Fish1; double Fish2; double FishA; double FishB; double FishC;double FishD;double FishE;double FishF;double FishG;double FishH;double FishI;double FishJ; double sepruleA; double sepruleB; double sepruleC;double sepruleD;double sepruleE;double sepruleF;double sepruleG;double sepruleH;double sepruleI;double sepruleJ; double seprule; NumericVector seprule_comb_vec(Nyrs); NumericVector d(10); //10 is right because it is for the 10 areas for (int i=0; i<10;i++){d[i]=pow((i+1)*100,dpow);} NumericVector Effort(10); NumericVector Effort_norm(10); NumericVector vec(10); NumericVector vec_CC(10);//for catch cascading NumericVector spawnersA(21);NumericVector spawnersB(21);NumericVector spawnersC(21);NumericVector spawnersD(21);NumericVector spawnersE(21); NumericVector spawnersF(21);NumericVector spawnersG(21);NumericVector spawnersH(21);NumericVector spawnersI(21);NumericVector spawnersJ(21); NumericVector spawnersA_rounded(21);NumericVector spawnersB_rounded(21);NumericVector spawnersC_rounded(21);NumericVector spawnersD_rounded(21);NumericVector spawnersE_rounded(21); NumericVector spawnersF_rounded(21);NumericVector spawnersG_rounded(21);NumericVector spawnersH_rounded(21);NumericVector spawnersI_rounded(21);NumericVector spawnersJ_rounded(21); NumericVector spawnersA_rounded_star(21);NumericVector spawnersB_rounded_star(21);NumericVector spawnersC_rounded_star(21);NumericVector spawnersD_rounded_star(21);NumericVector spawnersE_rounded_star(21); NumericVector spawnersF_rounded_star(21);NumericVector spawnersG_rounded_star(21);NumericVector spawnersH_rounded_star(21);NumericVector spawnersI_rounded_star(21);NumericVector spawnersJ_rounded_star(21); double pick_spawner_mal; double pick_spawner_fem; //int pick2; //this might not be necessary NumericVector which_spawnersA; int some_integerA; int some_integerB; int some_integerC;int some_integerD;int some_integerE;int some_integerF;int some_integerG;int some_integerH;int some_integerI;int some_integerJ; List tmpListA; List tmpListB; List tmpListC;List tmpListD;List tmpListE;List tmpListF;List tmpListG;List tmpListH;List tmpListI;List tmpListJ; int pickM; int pickF; double pick_moverA; double pick_moverB; int mov_pick1; int mov_pick2; NumericVector from_vec; NumericMatrix moverA2Brow1; NumericMatrix moverA2Brow2; NumericMatrix moverB2Arow1; NumericMatrix moverB2Arow2; NumericMatrix Replace1; NumericMatrix Replace2; double replace_moverA; double replace_moverB; int rep_pick1; int rep_pick2; NumericVector Nvec_temp; NumericMatrix RecA; NumericMatrix RecB; NumericMatrix RecC;NumericMatrix RecD;NumericMatrix RecE;NumericMatrix RecF;NumericMatrix RecG;NumericMatrix RecH;NumericMatrix RecI;NumericMatrix RecJ; NumericMatrix recA; NumericMatrix recB; NumericMatrix recC;NumericMatrix recD;NumericMatrix recE;NumericMatrix recF;NumericMatrix recG;NumericMatrix recH;NumericMatrix recI;NumericMatrix recJ; double rand_doub; int samp_mom; int samp_dad; int rand_allele; List tmpList_recsA(2); List tmpList_recsB(2); List tmpList_recsC(2);List tmpList_recsD(2);List tmpList_recsE(2);List tmpList_recsF(2);List tmpList_recsG(2);List tmpList_recsH(2);List tmpList_recsI(2);List tmpList_recsJ(2); //NumericVector FSH; //IntegerVector FSH_int; List tmpdoneA(2); List tmpdoneB(2); List tmpdoneC(2);List tmpdoneD(2);List tmpdoneE(2);List tmpdoneF(2);List tmpdoneG(2);List tmpdoneH(2);List tmpdoneI(2);List tmpdoneJ(2); NumericMatrix tmpageA_fem; NumericMatrix tmpageA_mal; NumericMatrix tmpageB_fem; NumericMatrix tmpageB_mal; NumericMatrix tmpageC_fem; NumericMatrix tmpageC_mal; NumericMatrix tmpageD_fem; NumericMatrix tmpageD_mal; NumericMatrix tmpageE_fem; NumericMatrix tmpageE_mal; NumericMatrix tmpageF_fem; NumericMatrix tmpageF_mal; NumericMatrix tmpageG_fem; NumericMatrix tmpageG_mal; NumericMatrix tmpageH_fem; NumericMatrix tmpageH_mal; NumericMatrix tmpageI_fem; NumericMatrix tmpageI_mal; NumericMatrix tmpageJ_fem; NumericMatrix tmpageJ_mal; NumericVector all_richA; NumericVector all_richB; NumericVector all_richC;NumericVector all_richD;NumericVector all_richE;NumericVector all_richF;NumericVector all_richG;NumericVector all_richH;NumericVector all_richI;NumericVector all_richJ; NumericVector arichA(Nyrs); NumericVector arichB(Nyrs); NumericVector arichC(Nyrs);NumericVector arichD(Nyrs);NumericVector arichE(Nyrs);NumericVector arichF(Nyrs);NumericVector arichG(Nyrs);NumericVector arichH(Nyrs);NumericVector arichI(Nyrs);NumericVector arichJ(Nyrs); NumericMatrix allspawnA(Nyrs,21); NumericMatrix allspawnB(Nyrs,21); NumericMatrix allspawnC(Nyrs,21);NumericMatrix allspawnD(Nyrs,21);NumericMatrix allspawnE(Nyrs,21);NumericMatrix allspawnF(Nyrs,21);NumericMatrix allspawnG(Nyrs,21);NumericMatrix allspawnH(Nyrs,21);NumericMatrix allspawnI(Nyrs,21);NumericMatrix allspawnJ(Nyrs,21); NumericVector N_aA(21); NumericVector N_aB(21); NumericVector N_aC(21);NumericVector N_aD(21);NumericVector N_aE(21);NumericVector N_aF(21);NumericVector N_aG(21);NumericVector N_aH(21);NumericVector N_aI(21);NumericVector N_aJ(21); List spawn_malA(Nyrs); List spawn_femA(Nyrs); List spawn_malB(Nyrs); List spawn_femB(Nyrs); List spawn_malC(Nyrs); List spawn_femC(Nyrs); List spawn_malD(Nyrs); List spawn_femD(Nyrs); List spawn_malE(Nyrs); List spawn_femE(Nyrs); List spawn_malF(Nyrs); List spawn_femF(Nyrs); List spawn_malG(Nyrs); List spawn_femG(Nyrs); List spawn_malH(Nyrs); List spawn_femH(Nyrs); List spawn_malI(Nyrs); List spawn_femI(Nyrs); List spawn_malJ(Nyrs); List spawn_femJ(Nyrs); NumericMatrix tmpMatrixA_1; NumericMatrix tmpMatrixA_2; NumericMatrix tmpMatrixB_1; NumericMatrix tmpMatrixB_2; NumericMatrix tmpMatrixC_1; NumericMatrix tmpMatrixC_2; NumericMatrix tmpMatrixD_1; NumericMatrix tmpMatrixD_2; NumericMatrix tmpMatrixE_1; NumericMatrix tmpMatrixE_2; NumericMatrix tmpMatrixF_1; NumericMatrix tmpMatrixF_2; NumericMatrix tmpMatrixG_1; NumericMatrix tmpMatrixG_2; NumericMatrix tmpMatrixH_1; NumericMatrix tmpMatrixH_2; NumericMatrix tmpMatrixI_1; NumericMatrix tmpMatrixI_2; NumericMatrix tmpMatrixJ_1; NumericMatrix tmpMatrixJ_2; IntegerVector NspawnersA(Nyrs); IntegerVector NspawnersB(Nyrs); IntegerVector NspawnersC(Nyrs); IntegerVector NspawnersD(Nyrs); IntegerVector NspawnersE(Nyrs); IntegerVector NspawnersF(Nyrs); IntegerVector NspawnersG(Nyrs); IntegerVector NspawnersH(Nyrs); IntegerVector NspawnersI(Nyrs); IntegerVector NspawnersJ(Nyrs); NumericMatrix spawn_malA_thisyr; //twice as many genotypes as spawners NumericMatrix spawn_femA_thisyr; NumericMatrix spawn_malB_thisyr; NumericMatrix spawn_femB_thisyr; NumericMatrix spawn_malC_thisyr; //twice as many genotypes as spawners NumericMatrix spawn_femC_thisyr; NumericMatrix spawn_malD_thisyr; NumericMatrix spawn_femD_thisyr; NumericMatrix spawn_malE_thisyr; //twice as many genotypes as spawners NumericMatrix spawn_femE_thisyr; NumericMatrix spawn_malF_thisyr; NumericMatrix spawn_femF_thisyr; NumericMatrix spawn_malG_thisyr; //twice as many genotypes as spawners NumericMatrix spawn_femG_thisyr; NumericMatrix spawn_malH_thisyr; NumericMatrix spawn_femH_thisyr; NumericMatrix spawn_malI_thisyr; //twice as many genotypes as spawners NumericMatrix spawn_femI_thisyr; NumericMatrix spawn_malJ_thisyr; NumericMatrix spawn_femJ_thisyr; NumericMatrix tmpMatrixA_1plus; NumericMatrix tmpMatrixA_2plus; NumericMatrix tmpMatrixB_1plus; NumericMatrix tmpMatrixB_2plus; NumericMatrix tmpMatrixC_1plus; NumericMatrix tmpMatrixC_2plus; NumericMatrix tmpMatrixD_1plus; NumericMatrix tmpMatrixD_2plus; NumericMatrix tmpMatrixE_1plus; NumericMatrix tmpMatrixE_2plus; NumericMatrix tmpMatrixF_1plus; NumericMatrix tmpMatrixF_2plus; NumericMatrix tmpMatrixG_1plus; NumericMatrix tmpMatrixG_2plus; NumericMatrix tmpMatrixH_1plus; NumericMatrix tmpMatrixH_2plus; NumericMatrix tmpMatrixI_1plus; NumericMatrix tmpMatrixI_2plus; NumericMatrix tmpMatrixJ_1plus; NumericMatrix tmpMatrixJ_2plus; List tmpList_plusgroupA(2); List tmpList_plusgroupB(2); List tmpList_plusgroupC(2); List tmpList_plusgroupD(2); List tmpList_plusgroupE(2); List tmpList_plusgroupF(2); List tmpList_plusgroupG(2); List tmpList_plusgroupH(2); List tmpList_plusgroupI(2); List tmpList_plusgroupJ(2); IntegerVector mal_spawnersA; //after random mating **also make sure you get the right row IntegerVector fem_spawnersA; //after random mating IntegerVector mal_spawnersB; //after random mating **also make sure you get the right row IntegerVector fem_spawnersB; //after random mating IntegerVector mal_spawnersC; //after random mating **also make sure you get the right row IntegerVector fem_spawnersC; //after random mating IntegerVector mal_spawnersD; //after random mating **also make sure you get the right row IntegerVector fem_spawnersD; //after random mating IntegerVector mal_spawnersE; //after random mating **also make sure you get the right row IntegerVector fem_spawnersE; //after random mating IntegerVector mal_spawnersF; //after random mating **also make sure you get the right row IntegerVector fem_spawnersF; //after random mating IntegerVector mal_spawnersG; //after random mating **also make sure you get the right row IntegerVector fem_spawnersG; //after random mating IntegerVector mal_spawnersH; //after random mating **also make sure you get the right row IntegerVector fem_spawnersH; //after random mating IntegerVector mal_spawnersI; //after random mating **also make sure you get the right row IntegerVector fem_spawnersI; //after random mating IntegerVector mal_spawnersJ; //after random mating **also make sure you get the right row IntegerVector fem_spawnersJ; //after random mating IntegerVector spawner_listM; //list of ones to choose from IntegerVector spawner_listF; NumericVector N_init_a1_star(21); NumericVector N_init_a2_star(21); NumericVector N_init_a3_star(21); NumericVector N_init_a4_star(21); NumericVector N_init_a5_star(21); NumericVector N_init_a6_star(21); NumericVector N_init_a7_star(21); NumericVector N_init_a8_star(21); NumericVector N_init_a9_star(21); NumericVector N_init_a10_star(21); NumericVector NtotalA(Nyrs+20); NumericVector NtotalB(Nyrs+20); NumericVector NtotalC(Nyrs+20); NumericVector NtotalD(Nyrs+20); NumericVector NtotalE(Nyrs+20); NumericVector NtotalF(Nyrs+20); NumericVector NtotalG(Nyrs+20); NumericVector NtotalH(Nyrs+20); NumericVector NtotalI(Nyrs+20); NumericVector NtotalJ(Nyrs+20); NumericVector SB40_vecA(Nyrs); NumericVector SB40_vecB(Nyrs); NumericVector SB40_vecC(Nyrs); NumericVector SB40_vecD(Nyrs); NumericVector SB40_vecE(Nyrs); NumericVector SB40_vecF(Nyrs); NumericVector SB40_vecG(Nyrs); NumericVector SB40_vecH(Nyrs); NumericVector SB40_vecI(Nyrs); NumericVector SB40_vecJ(Nyrs); int migs; int migs_leftovers; //END TRY FST CALCS //********************************************************************************************** //mgmt_counter //0=combined //1=separated (totally) //2=two mgmt units between 5|6 //3=split in two between 2|3 //4=tier 5 and combined //5=catch cascading //6=split at 8|9 //for (int mgmt_counter=0;mgmt_counter=1){ eta_jyA=as(rnorm(1,0,pow(sigma_ja2A,0.5))); eta_jyB=as(rnorm(1,0,pow(sigma_ja2B,0.5))); eta_jyC=as(rnorm(1,0,pow(sigma_ja2C,0.5))); eta_jyD=as(rnorm(1,0,pow(sigma_ja2D,0.5))); eta_jyE=as(rnorm(1,0,pow(sigma_ja2E,0.5))); eta_jyF=as(rnorm(1,0,pow(sigma_ja2F,0.5))); eta_jyG=as(rnorm(1,0,pow(sigma_ja2G,0.5))); eta_jyH=as(rnorm(1,0,pow(sigma_ja2H,0.5))); eta_jyI=as(rnorm(1,0,pow(sigma_ja2I,0.5))); eta_jyJ=as(rnorm(1,0,pow(sigma_ja2J,0.5))); eta_jy=as(rnorm(1,0,pow(sigma_ja2_tot,0.5))); epsilonA[k]=rho*epsilonA[k-1]+pow((1-pow(rho,2)),0.5)*eta_jyA; epsilonB[k]=rho*epsilonB[k-1]+pow((1-pow(rho,2)),0.5)*eta_jyB; epsilonC[k]=rho*epsilonC[k-1]+pow((1-pow(rho,2)),0.5)*eta_jyC; epsilonD[k]=rho*epsilonD[k-1]+pow((1-pow(rho,2)),0.5)*eta_jyD; epsilonE[k]=rho*epsilonE[k-1]+pow((1-pow(rho,2)),0.5)*eta_jyE; epsilonF[k]=rho*epsilonF[k-1]+pow((1-pow(rho,2)),0.5)*eta_jyF; epsilonG[k]=rho*epsilonG[k-1]+pow((1-pow(rho,2)),0.5)*eta_jyG; epsilonH[k]=rho*epsilonH[k-1]+pow((1-pow(rho,2)),0.5)*eta_jyH; epsilonI[k]=rho*epsilonI[k-1]+pow((1-pow(rho,2)),0.5)*eta_jyI; epsilonJ[k]=rho*epsilonJ[k-1]+pow((1-pow(rho,2)),0.5)*eta_jyJ; epsilon[k]=rho*epsilon[k-1]+pow((1-pow(rho,2)),0.5)*eta_jy; } //end of k>=1 loop psiA=as(rnorm(1,0,pow(sigma_js2A,0.5))); psiB=as(rnorm(1,0,pow(sigma_js2B,0.5))); psiC=as(rnorm(1,0,pow(sigma_js2C,0.5))); psiD=as(rnorm(1,0,pow(sigma_js2D,0.5))); psiE=as(rnorm(1,0,pow(sigma_js2E,0.5))); psiF=as(rnorm(1,0,pow(sigma_js2F,0.5))); psiG=as(rnorm(1,0,pow(sigma_js2G,0.5))); psiH=as(rnorm(1,0,pow(sigma_js2H,0.5))); psiI=as(rnorm(1,0,pow(sigma_js2I,0.5))); psiJ=as(rnorm(1,0,pow(sigma_js2J,0.5))); N_hat_aA = N_init_a1*exp(epsilonA[k]-sigma_ja2A/2); N_hat_aB = N_init_a2*exp(epsilonB[k]-sigma_ja2B/2); N_hat_aC = N_init_a3*exp(epsilonC[k]-sigma_ja2C/2); N_hat_aD = N_init_a4*exp(epsilonD[k]-sigma_ja2D/2); N_hat_aE = N_init_a5*exp(epsilonE[k]-sigma_ja2E/2); N_hat_aF = N_init_a6*exp(epsilonF[k]-sigma_ja2F/2); N_hat_aG = N_init_a7*exp(epsilonG[k]-sigma_ja2G/2); N_hat_aH = N_init_a8*exp(epsilonH[k]-sigma_ja2H/2); N_hat_aI = N_init_a9*exp(epsilonI[k]-sigma_ja2I/2); N_hat_aJ = N_init_a10*exp(epsilonJ[k]-sigma_ja2J/2); R_l0A_hat[k]=N_hat_aA[0]; //incorporate error into recruitment observation R_l0B_hat[k]=N_hat_aB[0];//if problem go back to old way this should work though R_l0C_hat[k]=N_hat_aC[0];//because above N_hat_aA = N_init_a1*exp(epsilonA[k]-sigma_ja2/2); R_l0D_hat[k]=N_hat_aD[0]; R_l0E_hat[k]=N_hat_aE[0]; R_l0F_hat[k]=N_hat_aF[0]; R_l0G_hat[k]=N_hat_aG[0]; R_l0H_hat[k]=N_hat_aH[0]; R_l0I_hat[k]=N_hat_aI[0]; R_l0J_hat[k]=N_hat_aJ[0]; N_hat_aA_S = N_init_a1*exp(psiA-(sigma_js2A/2)); N_hat_aB_S = N_init_a2*exp(psiB-(sigma_js2B/2)); N_hat_aC_S = N_init_a3*exp(psiC-(sigma_js2C/2)); N_hat_aD_S = N_init_a4*exp(psiD-(sigma_js2D/2)); N_hat_aE_S = N_init_a5*exp(psiE-(sigma_js2E/2)); N_hat_aF_S = N_init_a6*exp(psiF-(sigma_js2F/2)); N_hat_aG_S = N_init_a7*exp(psiG-(sigma_js2G/2)); N_hat_aH_S = N_init_a8*exp(psiH-(sigma_js2H/2)); N_hat_aI_S = N_init_a9*exp(psiI-(sigma_js2I/2)); N_hat_aJ_S = N_init_a10*exp(psiJ-(sigma_js2J/2)); N_init_a=N_init_a1+N_init_a2+N_init_a3+N_init_a4+N_init_a5+N_init_a6+N_init_a7+N_init_a8+N_init_a9+N_init_a10; N_hat_a=N_hat_aA+N_hat_aB+N_hat_aC+N_hat_aD+N_hat_aE+N_hat_aF+N_hat_aG+N_hat_aH+N_hat_aI+N_hat_aJ; psi_star=as(rnorm(1,0,pow(sigma_js2_tot,0.5))); //these are SSB with assessment error tempA=(N_hat_aA*Wt_fem*Q); SSB_hatA[k]=0.5*std::accumulate(tempA.begin(),tempA.end(),0.0); //spawning biomass of popA tempB=(N_hat_aB*Wt_fem*Q); SSB_hatB[k]=0.5*std::accumulate(tempB.begin(),tempB.end(),0.0); //spawning biomass of popB tempC=(N_hat_aC*Wt_fem*Q); SSB_hatC[k]=0.5*std::accumulate(tempC.begin(),tempC.end(),0.0); //spawning biomass of popC tempD=(N_hat_aD*Wt_fem*Q); SSB_hatD[k]=0.5*std::accumulate(tempD.begin(),tempD.end(),0.0); //spawning biomass of popD tempE=(N_hat_aE*Wt_fem*Q); SSB_hatE[k]=0.5*std::accumulate(tempE.begin(),tempE.end(),0.0); //spawning biomass of popE tempF=(N_hat_aF*Wt_fem*Q); SSB_hatF[k]=0.5*std::accumulate(tempF.begin(),tempF.end(),0.0); //spawning biomass of popF tempG=(N_hat_aG*Wt_fem*Q); SSB_hatG[k]=0.5*std::accumulate(tempG.begin(),tempG.end(),0.0); //spawning biomass of popG tempH=(N_hat_aH*Wt_fem*Q); SSB_hatH[k]=0.5*std::accumulate(tempH.begin(),tempH.end(),0.0); //spawning biomass of popH tempI=(N_hat_aI*Wt_fem*Q); SSB_hatI[k]=0.5*std::accumulate(tempI.begin(),tempI.end(),0.0); //spawning biomass of popI tempJ=(N_hat_aJ*Wt_fem*Q); SSB_hatJ[k]=0.5*std::accumulate(tempJ.begin(),tempJ.end(),0.0); //spawning biomass of popJ SSB_hat[k]= (SSB_hatA[k]+SSB_hatB[k]+SSB_hatC[k]+SSB_hatD[k]+SSB_hatE[k]+ SSB_hatF[k]+SSB_hatG[k]+SSB_hatH[k]+SSB_hatI[k]+SSB_hatJ[k]); //these are SSB with survey error tempA=(N_hat_aA_S*Wt_fem*Q); SSB_hatA_S[k]=0.5*std::accumulate(tempA.begin(),tempA.end(),0.0); //spawning biomass of popA tempB=(N_hat_aB_S*Wt_fem*Q); SSB_hatB_S[k]=0.5*std::accumulate(tempB.begin(),tempB.end(),0.0); //spawning biomass of popB tempC=(N_hat_aC_S*Wt_fem*Q); SSB_hatC_S[k]=0.5*std::accumulate(tempC.begin(),tempC.end(),0.0); //spawning biomass of popC tempD=(N_hat_aD_S*Wt_fem*Q); SSB_hatD_S[k]=0.5*std::accumulate(tempD.begin(),tempD.end(),0.0); //spawning biomass of popD tempE=(N_hat_aE_S*Wt_fem*Q); SSB_hatE_S[k]=0.5*std::accumulate(tempE.begin(),tempE.end(),0.0); //spawning biomass of popE tempF=(N_hat_aF_S*Wt_fem*Q); SSB_hatF_S[k]=0.5*std::accumulate(tempF.begin(),tempF.end(),0.0); //spawning biomass of popF tempG=(N_hat_aG_S*Wt_fem*Q); SSB_hatG_S[k]=0.5*std::accumulate(tempG.begin(),tempG.end(),0.0); //spawning biomass of popG tempH=(N_hat_aH_S*Wt_fem*Q); SSB_hatH_S[k]=0.5*std::accumulate(tempH.begin(),tempH.end(),0.0); //spawning biomass of popH tempI=(N_hat_aI_S*Wt_fem*Q); SSB_hatI_S[k]=0.5*std::accumulate(tempI.begin(),tempI.end(),0.0); //spawning biomass of popI tempJ=(N_hat_aJ_S*Wt_fem*Q); SSB_hatJ_S[k]=0.5*std::accumulate(tempJ.begin(),tempJ.end(),0.0); //spawning biomass of popJ SSB_hat_S[k]= (SSB_hatA_S[k]+SSB_hatB_S[k]+SSB_hatC_S[k]+SSB_hatD_S[k]+SSB_hatE_S[k]+ SSB_hatF_S[k]+SSB_hatG_S[k]+SSB_hatH_S[k]+SSB_hatI_S[k]+SSB_hatJ_S[k]); SB40_vecA[k]=SPR_init_F40*(std::accumulate(R_l0A_hat.begin(),R_l0A_hat.begin()+k,0.0)/k); SB40_vecB[k]=SPR_init_F40*(std::accumulate(R_l0B_hat.begin(),R_l0B_hat.begin()+k,0.0)/k); SB40_vecC[k]=SPR_init_F40*(std::accumulate(R_l0C_hat.begin(),R_l0C_hat.begin()+k,0.0)/k); SB40_vecD[k]=SPR_init_F40*(std::accumulate(R_l0D_hat.begin(),R_l0D_hat.begin()+k,0.0)/k); SB40_vecE[k]=SPR_init_F40*(std::accumulate(R_l0E_hat.begin(),R_l0E_hat.begin()+k,0.0)/k); SB40_vecF[k]=SPR_init_F40*(std::accumulate(R_l0F_hat.begin(),R_l0F_hat.begin()+k,0.0)/k); SB40_vecG[k]=SPR_init_F40*(std::accumulate(R_l0G_hat.begin(),R_l0G_hat.begin()+k,0.0)/k); SB40_vecH[k]=SPR_init_F40*(std::accumulate(R_l0H_hat.begin(),R_l0H_hat.begin()+k,0.0)/k); SB40_vecI[k]=SPR_init_F40*(std::accumulate(R_l0I_hat.begin(),R_l0I_hat.begin()+k,0.0)/k); SB40_vecJ[k]=SPR_init_F40*(std::accumulate(R_l0J_hat.begin(),R_l0J_hat.begin()+k,0.0)/k); if (k<10){ FishA=fishmort; FishB=fishmort; FishC=fishmort; FishD=fishmort; FishE=fishmort; FishF=fishmort; FishG=fishmort; FishH=fishmort; FishI=fishmort; FishJ=fishmort; TAC[k]=0; seprule_comb_vec[k]=0; } //end of if k<10 if (k>110){ FishA=0; FishB=0; FishC=0; FishD=0; FishE=0; FishF=0; FishG=0; FishH=0; FishI=0; FishJ=0; TAC[k]=0; seprule_comb_vec[k]=0; } //end of if k<110 if ((k>=10)&(k<=110)){ //cases where mgmt is "sep". Here I use 40:10 ctrl rule (correct version). if(mgmt_counter==1){ sepruleA = SSB_hatA[k]/(SPR_init_F40*(std::accumulate(R_l0A_hat.begin(),R_l0A_hat.begin()+k,0.0)/k)); sepruleB = SSB_hatB[k]/(SPR_init_F40*(std::accumulate(R_l0B_hat.begin(),R_l0B_hat.begin()+k,0.0)/k)); sepruleC = SSB_hatC[k]/(SPR_init_F40*(std::accumulate(R_l0C_hat.begin(),R_l0C_hat.begin()+k,0.0)/k)); sepruleD = SSB_hatD[k]/(SPR_init_F40*(std::accumulate(R_l0D_hat.begin(),R_l0D_hat.begin()+k,0.0)/k)); sepruleE = SSB_hatE[k]/(SPR_init_F40*(std::accumulate(R_l0E_hat.begin(),R_l0E_hat.begin()+k,0.0)/k)); sepruleF = SSB_hatF[k]/(SPR_init_F40*(std::accumulate(R_l0F_hat.begin(),R_l0F_hat.begin()+k,0.0)/k)); sepruleG = SSB_hatG[k]/(SPR_init_F40*(std::accumulate(R_l0G_hat.begin(),R_l0G_hat.begin()+k,0.0)/k)); sepruleH = SSB_hatH[k]/(SPR_init_F40*(std::accumulate(R_l0H_hat.begin(),R_l0H_hat.begin()+k,0.0)/k)); sepruleI = SSB_hatI[k]/(SPR_init_F40*(std::accumulate(R_l0I_hat.begin(),R_l0I_hat.begin()+k,0.0)/k)); sepruleJ = SSB_hatJ[k]/(SPR_init_F40*(std::accumulate(R_l0J_hat.begin(),R_l0J_hat.begin()+k,0.0)/k)); if (sepruleA>=1){FishA = fishmort;} if (sepruleA<=0.5){FishA = 0;} if ((sepruleA<1)&(sepruleA>0.5)) {FishA = fishmort*(sepruleA-0.05)*(1/0.95);}//CHANGED HERE if (sepruleB>=1){FishB = fishmort;} if (sepruleB<=0.5){FishB = 0;} if ((sepruleB<1)&(sepruleB>0.5)){FishB = fishmort*(sepruleB-0.05)*(1/0.95);}//CHANGED HERE if (sepruleC>=1){FishC = fishmort;} if (sepruleC<=0.5){FishC = 0;} if ((sepruleC<1)&(sepruleC>0.5)){FishC = fishmort*(sepruleC-0.05)*(1/0.95);}//CHANGED HERE if (sepruleD>=1){FishD = fishmort;} if (sepruleD<=0.5){FishD = 0;} if ((sepruleD<1)&(sepruleD>0.5)){FishD= fishmort*(sepruleD-0.05)*(1/0.95);}//CHANGED HERE if (sepruleE>=1){FishE = fishmort;} if (sepruleE<=0.5){FishE = 0;} if ((sepruleE<1)&(sepruleE>0.5)){FishE = fishmort*(sepruleE-0.05)*(1/0.95);}//CHANGED HERE if (sepruleF>=1){FishF = fishmort;} if (sepruleF<=0.5){FishF = 0;} if ((sepruleF<1)&(sepruleF>0.5)){FishF= fishmort*(sepruleF-0.05)*(1/0.95);}//CHANGED HERE if (sepruleG>=1){FishG = fishmort;} if (sepruleG<=0.5){FishG = 0;} if ((sepruleG<1)&(sepruleG>0.5)){FishG = fishmort*(sepruleG-0.05)*(1/0.95);}//CHANGED HERE if (sepruleH>=1){FishH = fishmort;} if (sepruleH<=0.5){FishH = 0;} if ((sepruleH<1)&(sepruleH>0.5)){FishH = fishmort*(sepruleH-0.05)*(1/0.95);}//CHANGED HERE if (sepruleI>=1){FishI = fishmort;} if (sepruleI<=0.5){FishI = 0;} if ((sepruleI<1)&(sepruleI>0.5)){FishI = fishmort*(sepruleI-0.05)*(1/0.95);}//CHANGED HERE if (sepruleJ>=1){FishJ = fishmort;} if (sepruleJ<=0.5){FishJ = 0;} if ((sepruleJ<1)&(sepruleJ>0.5)){FishJ = fishmort*(sepruleJ-0.05)*(1/0.95);}//CHANGED HERE //in Pcod only starting here }//end of if mgmt==sep separated by all 10 spatial areas //here is combined management if((mgmt_counter==0)|(mgmt_counter==4)){ seprule = (SSB_hat[k])/(SPR_init_F40*( (std::accumulate(R_l0A_hat.begin(),R_l0A_hat.begin()+k,0.0) +std::accumulate(R_l0B_hat.begin(),R_l0B_hat.begin()+k,0.0) +std::accumulate(R_l0C_hat.begin(),R_l0C_hat.begin()+k,0.0) +std::accumulate(R_l0D_hat.begin(),R_l0D_hat.begin()+k,0.0) +std::accumulate(R_l0E_hat.begin(),R_l0E_hat.begin()+k,0.0) +std::accumulate(R_l0F_hat.begin(),R_l0F_hat.begin()+k,0.0) +std::accumulate(R_l0G_hat.begin(),R_l0G_hat.begin()+k,0.0) +std::accumulate(R_l0H_hat.begin(),R_l0H_hat.begin()+k,0.0) +std::accumulate(R_l0I_hat.begin(),R_l0I_hat.begin()+k,0.0) +std::accumulate(R_l0J_hat.begin(),R_l0J_hat.begin()+k,0.0))/k)); if (seprule>=1){Fish = fishmort;} if (seprule<=0.5){Fish = 0;} if ((seprule<1)&(seprule>0.5)) {Fish = fishmort*(seprule-0.05)*(1/0.95);}//CHANGED HERE seprule_comb_vec[k]=seprule; temp14 = ((Wt_fem+Wt_mal)/2)*(N_hat_a)*((fishsel*Fish)/(fishsel*Fish+M))*(1-exp(-(fishsel*Fish+M)));//N_hat_a is all areas TAC[k]=std::accumulate(temp14.begin(),temp14.end(),0.0); //Total TAC just based on total number assessed if (mgmt_counter==4){TAC[k]=M*0.21*SSB_hat[k];} //assessment 2 function below Wt=(Wt_fem+Wt_mal)/2; //establish d which is the effort vector based on true number in each area //for (int i=0; i<10;i++){d[i]=(i+1)*100;} temp16a=fishsel*Wt*N_init_a1; Effort[0]=(1/pow(d[0],dpow))*std::accumulate(temp16a.begin(),temp16a.end(),0.0); temp16b=fishsel*Wt*N_init_a2; Effort[1]=(1/pow(d[1],dpow))*std::accumulate(temp16b.begin(),temp16b.end(),0.0); temp16c=fishsel*Wt*N_init_a3; Effort[2]=(1/pow(d[2],dpow))*std::accumulate(temp16c.begin(),temp16c.end(),0.0); temp16d=fishsel*Wt*N_init_a4; Effort[3]=(1/pow(d[3],dpow))*std::accumulate(temp16d.begin(),temp16d.end(),0.0); temp16e=fishsel*Wt*N_init_a5; Effort[4]=(1/pow(d[4],dpow))*std::accumulate(temp16e.begin(),temp16e.end(),0.0); temp16f=fishsel*Wt*N_init_a6; Effort[5]=(1/pow(d[5],dpow))*std::accumulate(temp16f.begin(),temp16f.end(),0.0); temp16g=fishsel*Wt*N_init_a7; Effort[6]=(1/pow(d[6],dpow))*std::accumulate(temp16g.begin(),temp16g.end(),0.0); temp16h=fishsel*Wt*N_init_a8; Effort[7]=(1/pow(d[7],dpow))*std::accumulate(temp16h.begin(),temp16h.end(),0.0); temp16i=fishsel*Wt*N_init_a9; Effort[8]=(1/pow(d[8],dpow))*std::accumulate(temp16i.begin(),temp16i.end(),0.0); temp16j=fishsel*Wt*N_init_a10; Effort[9]=(1/pow(d[9],dpow))*std::accumulate(temp16j.begin(),temp16j.end(),0.0); Effort_tot=std::accumulate(Effort.begin(),Effort.end(),0.0); if(Effort_tot>0){//account for the possibility that all fish will die for (int i=0;i<10;i++){ Effort_norm[i]=Effort[i]/Effort_tot;}} else{Effort_norm=rep(0,10);} //REMOVE //Effort_norm=rep(0.1,10); //end of only in Pcod for(int i=0;i<10;i++){ vec[i]=Effort_norm[i]*TAC[k];} Effort_comb(k,_)=Effort_norm; //vec tells you how much fish you are actually going to catch in each area //so find the TAC in each area that gets you that many fish. //TAC in each area is the lower of the TAC calculated by the Effort*totalTAC //OR the exploitable biomass in that area Wt*fishsel*N_init_a1 TACA[k]=vec[0]; TACB[k]=vec[1]; TACC[k]=vec[2]; TACD[k]=vec[3]; TACE[k]=vec[4]; TACF[k]=vec[5]; TACG[k]=vec[6]; TACH[k]=vec[7]; TACI[k]=vec[8]; TACJ[k]=vec[9]; //calculate FishA through FishJ this is fishing effort that results in desired total TAC //first create vector IntegerVector FSH_int = seq_len(600)-1; NumericVector FSH(600); for (int i=0;i<600;i++){ FSH[i]=FSH_int[i]*0.01;} int best_TAC_A; int best_TAC_B; int best_TAC_C; int best_TAC_D; int best_TAC_E; int best_TAC_F; int best_TAC_G; int best_TAC_H; int best_TAC_I; int best_TAC_J; int high_TAC_A; int high_TAC_B; int high_TAC_C; int high_TAC_D; int high_TAC_E; int high_TAC_F; int high_TAC_G; int high_TAC_H; int high_TAC_I; int high_TAC_J; for (int i=0;i<600;i++){ temp20A = Wt*(N_hat_aA)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestA[i]=std::accumulate(temp20A.begin(),temp20A.end(),0.0); } double highestA = *std::max_element(TACtestA.begin(),TACtestA.end()); for(int i=0; i<600;i++){ if(TACtestA[i]<=TACA[k]){best_TAC_A = i;}if(TACtestA[i]==highestA){high_TAC_A = i;} } for (int i=0;i<600;i++){ temp20B = Wt*(N_hat_aB)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestB[i]=std::accumulate(temp20B.begin(),temp20B.end(),0.0); } double highestB = *std::max_element(TACtestB.begin(),TACtestB.end()); for(int i=0; i<600;i++){ if(TACtestB[i]<=TACB[k]){best_TAC_B = i;}if(TACtestB[i]==highestB){high_TAC_B = i;} } for (int i=0;i<600;i++){ temp20C = Wt*(N_hat_aC)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestC[i]=std::accumulate(temp20C.begin(),temp20C.end(),0.0); } double highestC = *std::max_element(TACtestC.begin(),TACtestC.end()); for(int i=0; i<600;i++){ if(TACtestC[i]<=TACC[k]){best_TAC_C = i;}if(TACtestC[i]==highestC){high_TAC_C = i;} } for (int i=0;i<600;i++){ temp20D = Wt*(N_hat_aD)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestD[i]=std::accumulate(temp20D.begin(),temp20D.end(),0.0); } double highestD = *std::max_element(TACtestD.begin(),TACtestD.end()); for(int i=0; i<600;i++){ if(TACtestD[i]<=TACD[k]){best_TAC_D = i;}if(TACtestD[i]==highestD){high_TAC_D = i;} } for (int i=0;i<600;i++){ temp20E = Wt*(N_hat_aE)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestE[i]=std::accumulate(temp20E.begin(),temp20E.end(),0.0); } double highestE = *std::max_element(TACtestE.begin(),TACtestE.end()); for(int i=0; i<600;i++){ if(TACtestE[i]<=TACE[k]){best_TAC_E = i;}if(TACtestE[i]==highestE){high_TAC_E = i;} } for (int i=0;i<600;i++){ temp20F = Wt*(N_hat_aF)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestF[i]=std::accumulate(temp20F.begin(),temp20F.end(),0.0);//TACtestF is a range of catches under different fishing effort } double highestF = *std::max_element(TACtestF.begin(),TACtestF.end());//highestF is the highest TAC for(int i=0; i<600;i++){ if(TACtestF[i]<=TACF[k]){best_TAC_F = i;}if(TACtestF[i]==highestF){high_TAC_F = i;} } for (int i=0;i<600;i++){ temp20G = Wt*(N_hat_aG)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestG[i]=std::accumulate(temp20G.begin(),temp20G.end(),0.0); } double highestG = *std::max_element(TACtestG.begin(),TACtestG.end()); for(int i=0; i<600;i++){ if(TACtestG[i]<=TACG[k]){best_TAC_G = i;}if(TACtestG[i]==highestG){high_TAC_G = i;} } for (int i=0;i<600;i++){ temp20H = Wt*(N_hat_aH)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestH[i]=std::accumulate(temp20H.begin(),temp20H.end(),0.0); } double highestH = *std::max_element(TACtestH.begin(),TACtestH.end()); for(int i=0; i<600;i++){ if(TACtestH[i]<=TACH[k]){best_TAC_H = i;}if(TACtestH[i]==highestH){high_TAC_H = i;} } for (int i=0;i<600;i++){ temp20I = Wt*(N_hat_aI)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestI[i]=std::accumulate(temp20I.begin(),temp20I.end(),0.0); } double highestI = *std::max_element(TACtestI.begin(),TACtestI.end()); for(int i=0; i<600;i++){ if(TACtestI[i]<=TACI[k]){best_TAC_I = i;}if(TACtestI[i]==highestI){high_TAC_I = i;} } for (int i=0;i<600;i++){ temp20J = Wt*(N_hat_aJ)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestJ[i]=std::accumulate(temp20J.begin(),temp20J.end(),0.0); } double highestJ = *std::max_element(TACtestJ.begin(),TACtestJ.end()); for(int i=0; i<600;i++){ if(TACtestJ[i]<=TACJ[k]){best_TAC_J = i;}if(TACtestJ[i]==highestJ){high_TAC_J = i;} } if(TACA[k]==0){FishA=0;}else{ if(FSH[best_TAC_A]<=highestA){FishA = FSH[best_TAC_A];} else {FishA = FSH[high_TAC_A];}} if(TACB[k]==0){FishB=0;}else{ if(FSH[best_TAC_B]<=highestB){FishB = FSH[best_TAC_B];} else {FishB = FSH[high_TAC_B];}} if(TACC[k]==0){FishC=0;}else{ if(FSH[best_TAC_C]<=highestC){FishC = FSH[best_TAC_C];} else {FishC = FSH[high_TAC_C];}} if(TACD[k]==0){FishD=0;}else{ if(FSH[best_TAC_D]<=highestD){FishD = FSH[best_TAC_D];} else {FishD = FSH[high_TAC_D];}} if(TACE[k]==0){FishE=0;}else{ if(FSH[best_TAC_E]<=highestE){FishE = FSH[best_TAC_E];} else {FishE = FSH[high_TAC_E];}} if(TACF[k]==0){FishF=0;}else{ if(FSH[best_TAC_F]<=highestF){FishF = FSH[best_TAC_F];} else {FishF = FSH[high_TAC_F];}} if(TACG[k]==0){FishG=0;}else{ if(FSH[best_TAC_G]<=highestG){FishG = FSH[best_TAC_G];} else {FishG = FSH[high_TAC_G];}} if(TACH[k]==0){FishH=0;}else{ if(FSH[best_TAC_H]<=highestH){FishH = FSH[best_TAC_H];} else {FishH = FSH[high_TAC_H];}} if(TACI[k]==0){FishI=0;}else{ if(FSH[best_TAC_I]<=highestI){FishI = FSH[best_TAC_I];} else {FishI = FSH[high_TAC_I];}} if(TACJ[k]==0){FishJ=0;}else{ if(FSH[best_TAC_J]<=highestJ){FishJ = FSH[best_TAC_J];} else {FishJ = FSH[high_TAC_J];}} } // end of combined all 10 areas managed as one //here is separated into 2 areas split 5|6 if(mgmt_counter==2){ sepruleA = (SSB_hatA[k]+SSB_hatB[k]+SSB_hatC[k]+SSB_hatD[k]+SSB_hatE[k])/ (SPR_init_F40*((std::accumulate(R_l0A_hat.begin(),R_l0A_hat.begin()+k,0.0)+ std::accumulate(R_l0B_hat.begin(),R_l0B_hat.begin()+k,0.0) +std::accumulate(R_l0C_hat.begin(),R_l0C_hat.begin()+k,0.0) +std::accumulate(R_l0D_hat.begin(),R_l0D_hat.begin()+k,0.0) +std::accumulate(R_l0E_hat.begin(),R_l0E_hat.begin()+k,0.0))/k)); sepruleF = (SSB_hatF[k]+SSB_hatG[k]+SSB_hatH[k]+SSB_hatI[k]+SSB_hatJ[k])/ (SPR_init_F40*((std::accumulate(R_l0F_hat.begin(),R_l0F_hat.begin()+k,0.0)+ std::accumulate(R_l0G_hat.begin(),R_l0G_hat.begin()+k,0.0) +std::accumulate(R_l0H_hat.begin(),R_l0H_hat.begin()+k,0.0) +std::accumulate(R_l0I_hat.begin(),R_l0I_hat.begin()+k,0.0) +std::accumulate(R_l0J_hat.begin(),R_l0J_hat.begin()+k,0.0))/k)); if (sepruleA>=1){Fish1 = fishmort;} if (sepruleA<=0.5){Fish1 = 0;} if ((sepruleA<1)&(sepruleA>0.5)) {Fish1 = fishmort*(sepruleA-0.05)*(1/0.95);}//Fish 1 is fishing pressure in whole area 1 if (sepruleF>=1){Fish2 = fishmort;} if (sepruleF<=0.5){Fish2 = 0;} if ((sepruleF<1)&(sepruleF>0.5)) {Fish2 = fishmort*(sepruleF-0.05)*(1/0.95);}//Fish 2 is fishing pressure in whole area 2 temp14 = ((Wt_fem+Wt_mal)/2)*(N_hat_aA+N_hat_aB+N_hat_aC+N_hat_aD+N_hat_aE)*((fishsel*Fish1)/(fishsel*Fish1+M))*(1-exp(-(fishsel*Fish1+M)));//number in left half TAC1[k]=std::accumulate(temp14.begin(),temp14.end(),0.0); temp14 = ((Wt_fem+Wt_mal)/2)*(N_hat_aF+N_hat_aG+N_hat_aH+N_hat_aI+N_hat_aJ)*((fishsel*Fish2)/(fishsel*Fish2+M))*(1-exp(-(fishsel*Fish2+M)));//number in left half TAC2[k]=std::accumulate(temp14.begin(),temp14.end(),0.0); //assessment 2 function below Wt=(Wt_fem+Wt_mal)/2; //establish d which is the effort vector // for (int i=0; i<10;i++){d[i]=(i+1)*100;} temp16a=fishsel*Wt*N_init_a1; Effort[0]=(1/pow(d[0],dpow))*std::accumulate(temp16a.begin(),temp16a.end(),0.0); temp16b=fishsel*Wt*N_init_a2; Effort[1]=(1/pow(d[1],dpow))*std::accumulate(temp16b.begin(),temp16b.end(),0.0); temp16c=fishsel*Wt*N_init_a3; Effort[2]=(1/pow(d[2],dpow))*std::accumulate(temp16c.begin(),temp16c.end(),0.0); temp16d=fishsel*Wt*N_init_a4; Effort[3]=(1/pow(d[3],dpow))*std::accumulate(temp16d.begin(),temp16d.end(),0.0); temp16e=fishsel*Wt*N_init_a5; Effort[4]=(1/pow(d[4],dpow))*std::accumulate(temp16e.begin(),temp16e.end(),0.0); temp16f=fishsel*Wt*N_init_a6; Effort[5]=(1/pow(d[5],dpow))*std::accumulate(temp16f.begin(),temp16f.end(),0.0); temp16g=fishsel*Wt*N_init_a7; Effort[6]=(1/pow(d[6],dpow))*std::accumulate(temp16g.begin(),temp16g.end(),0.0); temp16h=fishsel*Wt*N_init_a8; Effort[7]=(1/pow(d[7],dpow))*std::accumulate(temp16h.begin(),temp16h.end(),0.0); temp16i=fishsel*Wt*N_init_a9; Effort[8]=(1/pow(d[8],dpow))*std::accumulate(temp16i.begin(),temp16i.end(),0.0); temp16j=fishsel*Wt*N_init_a10; Effort[9]=(1/pow(d[9],dpow))*std::accumulate(temp16j.begin(),temp16j.end(),0.0); double Effort_tot1=Effort[0]+Effort[1]+Effort[2]+Effort[3]+Effort[4]; double Effort_tot2=Effort[5]+Effort[6]+Effort[7]+Effort[8]+Effort[9]; if(Effort_tot1>0){//account for the possibility that all fish will die for (int i=0;i<5;i++){ Effort_norm[i]=Effort[i]/Effort_tot1;}} else{Effort_norm[0]=0; Effort_norm[1]=0; Effort_norm[2]=0; Effort_norm[3]=0; Effort_norm[4]=0;} if(Effort_tot2>0){ for (int i=5;i<10;i++){ Effort_norm[i]=Effort[i]/Effort_tot2;}} else{Effort_norm[5]=0; Effort_norm[6]=0; Effort_norm[7]=0; Effort_norm[8]=0; Effort_norm[9]=0;} //REMOVE //Effort_norm=rep(0.2,10); vec[0]=Effort_norm[0]*TAC1[k];vec[1]=Effort_norm[1]*TAC1[k];vec[2]=Effort_norm[2]*TAC1[k]; vec[3]=Effort_norm[3]*TAC1[k]; vec[4]=Effort_norm[4]*TAC1[k]; vec[5]=Effort_norm[5]*TAC2[k];vec[6]=Effort_norm[6]*TAC2[k];vec[7]=Effort_norm[7]*TAC2[k]; vec[8]=Effort_norm[8]*TAC2[k]; vec[9]=Effort_norm[9]*TAC2[k]; Effort_56(k,_)=Effort_norm; TACA[k]=vec[0]; TACB[k]=vec[1]; TACC[k]=vec[2]; TACD[k]=vec[3]; TACE[k]=vec[4]; TACF[k]=vec[5]; TACG[k]=vec[6]; TACH[k]=vec[7]; TACI[k]=vec[8]; TACJ[k]=vec[9]; //calculate FishA-J this is fishing effort that results in desired TACA and TACF //first create vector IntegerVector FSH_int = seq_len(600)-1; NumericVector FSH(600); for (int i=0;i<600;i++){ FSH[i]=FSH_int[i]*0.01;} int best_TAC_A; int best_TAC_B; int best_TAC_C; int best_TAC_D; int best_TAC_E; int best_TAC_F; int best_TAC_G; int best_TAC_H; int best_TAC_I; int best_TAC_J; int high_TAC_A; int high_TAC_B; int high_TAC_C; int high_TAC_D; int high_TAC_E; int high_TAC_F; int high_TAC_G; int high_TAC_H; int high_TAC_I; int high_TAC_J; for (int i=0;i<600;i++){ temp20A = Wt*(N_hat_aA)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestA[i]=std::accumulate(temp20A.begin(),temp20A.end(),0.0); } double highestA = *std::max_element(TACtestA.begin(),TACtestA.end()); for(int i=0; i<600;i++){ if(TACtestA[i]<=TACA[k]){best_TAC_A = i;}if(TACtestA[i]==highestA){high_TAC_A = i;} } for (int i=0;i<600;i++){ temp20B = Wt*(N_hat_aB)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestB[i]=std::accumulate(temp20B.begin(),temp20B.end(),0.0); } double highestB = *std::max_element(TACtestB.begin(),TACtestB.end()); for(int i=0; i<600;i++){ if(TACtestB[i]<=TACB[k]){best_TAC_B = i;}if(TACtestB[i]==highestB){high_TAC_B = i;} } for (int i=0;i<600;i++){ temp20C = Wt*(N_hat_aC)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestC[i]=std::accumulate(temp20C.begin(),temp20C.end(),0.0); } double highestC = *std::max_element(TACtestC.begin(),TACtestC.end()); for(int i=0; i<600;i++){ if(TACtestC[i]<=TACC[k]){best_TAC_C = i;}if(TACtestC[i]==highestC){high_TAC_C = i;} } for (int i=0;i<600;i++){ temp20D = Wt*(N_hat_aD)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestD[i]=std::accumulate(temp20D.begin(),temp20D.end(),0.0); } double highestD = *std::max_element(TACtestD.begin(),TACtestD.end()); for(int i=0; i<600;i++){ if(TACtestD[i]<=TACD[k]){best_TAC_D = i;}if(TACtestD[i]==highestD){high_TAC_D = i;} } for (int i=0;i<600;i++){ temp20E = Wt*(N_hat_aE)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestE[i]=std::accumulate(temp20E.begin(),temp20E.end(),0.0); } double highestE = *std::max_element(TACtestE.begin(),TACtestE.end()); for(int i=0; i<600;i++){ if(TACtestE[i]<=TACE[k]){best_TAC_E = i;}if(TACtestE[i]==highestE){high_TAC_E = i;} } for (int i=0;i<600;i++){ temp20F = Wt*(N_hat_aF)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestF[i]=std::accumulate(temp20F.begin(),temp20F.end(),0.0); } double highestF = *std::max_element(TACtestF.begin(),TACtestF.end()); for(int i=0; i<600;i++){ if(TACtestF[i]<=TACF[k]){best_TAC_F = i;}if(TACtestF[i]==highestF){high_TAC_F = i;} } for (int i=0;i<600;i++){ temp20G = Wt*(N_hat_aG)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestG[i]=std::accumulate(temp20G.begin(),temp20G.end(),0.0); } double highestG = *std::max_element(TACtestG.begin(),TACtestG.end()); for(int i=0; i<600;i++){ if(TACtestG[i]<=TACG[k]){best_TAC_G = i;}if(TACtestG[i]==highestG){high_TAC_G = i;} } for (int i=0;i<600;i++){ temp20H = Wt*(N_hat_aH)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestH[i]=std::accumulate(temp20H.begin(),temp20H.end(),0.0); } double highestH = *std::max_element(TACtestH.begin(),TACtestH.end()); for(int i=0; i<600;i++){ if(TACtestH[i]<=TACH[k]){best_TAC_H = i;}if(TACtestH[i]==highestH){high_TAC_H = i;} } for (int i=0;i<600;i++){ temp20I = Wt*(N_hat_aI)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestI[i]=std::accumulate(temp20I.begin(),temp20I.end(),0.0); } double highestI = *std::max_element(TACtestI.begin(),TACtestI.end()); for(int i=0; i<600;i++){ if(TACtestI[i]<=TACI[k]){best_TAC_I = i;}if(TACtestI[i]==highestI){high_TAC_I = i;} } for (int i=0;i<600;i++){ temp20J = Wt*(N_hat_aJ)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestJ[i]=std::accumulate(temp20J.begin(),temp20J.end(),0.0); } double highestJ = *std::max_element(TACtestJ.begin(),TACtestJ.end()); for(int i=0; i<600;i++){ if(TACtestJ[i]<=TACJ[k]){best_TAC_J = i;}if(TACtestJ[i]==highestJ){high_TAC_J = i;} } if(TACA[k]==0){FishA=0;}else{ if(FSH[best_TAC_A]<=highestA){FishA = FSH[best_TAC_A];} else {FishA = FSH[high_TAC_A];}} if(TACB[k]==0){FishB=0;}else{ if(FSH[best_TAC_B]<=highestB){FishB = FSH[best_TAC_B];} else {FishB = FSH[high_TAC_B];}} if(TACC[k]==0){FishC=0;}else{ if(FSH[best_TAC_C]<=highestC){FishC = FSH[best_TAC_C];} else {FishC = FSH[high_TAC_C];}} if(TACD[k]==0){FishD=0;}else{ if(FSH[best_TAC_D]<=highestD){FishD = FSH[best_TAC_D];} else {FishD = FSH[high_TAC_D];}} if(TACE[k]==0){FishE=0;}else{ if(FSH[best_TAC_E]<=highestE){FishE = FSH[best_TAC_E];} else {FishE = FSH[high_TAC_E];}} if(TACF[k]==0){FishF=0;}else{ if(FSH[best_TAC_F]<=highestF){FishF = FSH[best_TAC_F];} else {FishF = FSH[high_TAC_F];}} if(TACG[k]==0){FishG=0;}else{ if(FSH[best_TAC_G]<=highestG){FishG = FSH[best_TAC_G];} else {FishG = FSH[high_TAC_G];}} if(TACH[k]==0){FishH=0;}else{ if(FSH[best_TAC_H]<=highestH){FishH = FSH[best_TAC_H];} else {FishH = FSH[high_TAC_H];}} if(TACI[k]==0){FishI=0;}else{ if(FSH[best_TAC_I]<=highestI){FishI = FSH[best_TAC_I];} else {FishI = FSH[high_TAC_I];}} if(TACJ[k]==0){FishJ=0;}else{ if(FSH[best_TAC_J]<=highestJ){FishJ = FSH[best_TAC_J];} else {FishJ = FSH[high_TAC_J];}} } // end of two pops //here is separated into 2 areas by numbers not spatially if(mgmt_counter==3|mgmt_counter==6){ if (mgmt_counter==3){ //this is 2|3 sepruleA = (SSB_hatA[k]+SSB_hatB[k])/ (SPR_init_F40*((std::accumulate(R_l0A_hat.begin(),R_l0A_hat.begin()+k,0.0)+ std::accumulate(R_l0B_hat.begin(),R_l0B_hat.begin()+k,0.0))/k)); sepruleF = (SSB_hatC[k]+SSB_hatD[k]+SSB_hatE[k]+SSB_hatF[k]+SSB_hatG[k]+SSB_hatH[k]+SSB_hatI[k]+SSB_hatJ[k])/ (SPR_init_F40*((std::accumulate(R_l0C_hat.begin(),R_l0C_hat.begin()+k,0.0) +std::accumulate(R_l0D_hat.begin(),R_l0D_hat.begin()+k,0.0) +std::accumulate(R_l0E_hat.begin(),R_l0E_hat.begin()+k,0.0) +std::accumulate(R_l0F_hat.begin(),R_l0F_hat.begin()+k,0.0) +std::accumulate(R_l0G_hat.begin(),R_l0G_hat.begin()+k,0.0) +std::accumulate(R_l0H_hat.begin(),R_l0H_hat.begin()+k,0.0) +std::accumulate(R_l0I_hat.begin(),R_l0I_hat.begin()+k,0.0) +std::accumulate(R_l0J_hat.begin(),R_l0J_hat.begin()+k,0.0))/k)); if (sepruleA>=1){Fish1 = fishmort;} if (sepruleA<=0.5){Fish1 = 0;} if ((sepruleA<1)&(sepruleA>0.5)) {Fish1 = fishmort*(sepruleA-0.05)*(1/0.95);}//Fish 1 is fishing pressure in whole area 1 if (sepruleF>=1){Fish2 = fishmort;} if (sepruleF<=0.5){Fish2 = 0;} if ((sepruleF<1)&(sepruleF>0.5)) {Fish2 = fishmort*(sepruleF-0.05)*(1/0.95);}//Fish 2 is fishing pressure in whole area 2 temp14 = ((Wt_fem+Wt_mal)/2)*(N_hat_aA+N_hat_aB)*((fishsel*Fish1)/(fishsel*Fish1+M))*(1-exp(-(fishsel*Fish1+M)));//number in left half TAC1[k]=std::accumulate(temp14.begin(),temp14.end(),0.0); temp14 = ((Wt_fem+Wt_mal)/2)*(N_hat_aC+N_hat_aD+N_hat_aE+N_hat_aF+N_hat_aG+N_hat_aH+N_hat_aI+N_hat_aJ)*((fishsel*Fish2)/(fishsel*Fish2+M))*(1-exp(-(fishsel*Fish2+M)));//number in left half TAC2[k]=std::accumulate(temp14.begin(),temp14.end(),0.0); //assessment 2 function below Wt=(Wt_fem+Wt_mal)/2; //establish d which is the effort vector // for (int i=0; i<10;i++){d[i]=(i+1)*100;} temp16a=fishsel*Wt*N_init_a1; Effort[0]=(1/pow(d[0],dpow))*std::accumulate(temp16a.begin(),temp16a.end(),0.0); temp16b=fishsel*Wt*N_init_a2; Effort[1]=(1/pow(d[1],dpow))*std::accumulate(temp16b.begin(),temp16b.end(),0.0); temp16c=fishsel*Wt*N_init_a3; Effort[2]=(1/pow(d[2],dpow))*std::accumulate(temp16c.begin(),temp16c.end(),0.0); temp16d=fishsel*Wt*N_init_a4; Effort[3]=(1/pow(d[3],dpow))*std::accumulate(temp16d.begin(),temp16d.end(),0.0); temp16e=fishsel*Wt*N_init_a5; Effort[4]=(1/pow(d[4],dpow))*std::accumulate(temp16e.begin(),temp16e.end(),0.0); temp16f=fishsel*Wt*N_init_a6; Effort[5]=(1/pow(d[5],dpow))*std::accumulate(temp16f.begin(),temp16f.end(),0.0); temp16g=fishsel*Wt*N_init_a7; Effort[6]=(1/pow(d[6],dpow))*std::accumulate(temp16g.begin(),temp16g.end(),0.0); temp16h=fishsel*Wt*N_init_a8; Effort[7]=(1/pow(d[7],dpow))*std::accumulate(temp16h.begin(),temp16h.end(),0.0); temp16i=fishsel*Wt*N_init_a9; Effort[8]=(1/pow(d[8],dpow))*std::accumulate(temp16i.begin(),temp16i.end(),0.0); temp16j=fishsel*Wt*N_init_a10; Effort[9]=(1/pow(d[9],dpow))*std::accumulate(temp16j.begin(),temp16j.end(),0.0); double Effort_tot1=Effort[0]+Effort[1]; for (int i=0;i<2;i++){ if(Effort_tot1>0){//account for the possibility that all fish will die Effort_norm[i]=Effort[i]/Effort_tot1;} else{Effort_norm[i]=0;}} double Effort_tot2=Effort[2]+Effort[3]+Effort[4]+Effort[5]+Effort[6]+Effort[7]+Effort[8]+Effort[9]; for (int i=2;i<10;i++){ if(Effort_tot2>0){//account for the possibility that all fish will die Effort_norm[i]=Effort[i]/Effort_tot2;} else{Effort_norm[i]=0;}} //REMOVE //Effort_norm[0]=0.5; //Effort_norm[1]=0.5; //Effort_norm[2]=0.125; //Effort_norm[3]=0.125; //Effort_norm[4]=0.125; //Effort_norm[5]=0.125; //Effort_norm[6]=0.125; //Effort_norm[7]=0.125; //Effort_norm[8]=0.125; //Effort_norm[9]=0.125; vec[0]=Effort_norm[0]*TAC1[k];vec[1]=Effort_norm[1]*TAC1[k];vec[2]=Effort_norm[2]*TAC2[k]; vec[3]=Effort_norm[3]*TAC2[k]; vec[4]=Effort_norm[4]*TAC2[k]; vec[5]=Effort_norm[5]*TAC2[k];vec[6]=Effort_norm[6]*TAC2[k];vec[7]=Effort_norm[7]*TAC2[k]; vec[8]=Effort_norm[8]*TAC2[k]; vec[9]=Effort_norm[9]*TAC2[k]; Effort_23(k,_)=Effort_norm; }//mgmt_counter==3 if (mgmt_counter==6){ //this is 8|9 sepruleA = (SSB_hatA[k]+SSB_hatB[k]+SSB_hatC[k]+SSB_hatD[k]+SSB_hatE[k]+SSB_hatF[k]+SSB_hatG[k]+SSB_hatH[k])/ (SPR_init_F40*((std::accumulate(R_l0A_hat.begin(),R_l0A_hat.begin()+k,0.0)+ std::accumulate(R_l0B_hat.begin(),R_l0B_hat.begin()+k,0.0) +std::accumulate(R_l0C_hat.begin(),R_l0C_hat.begin()+k,0.0) +std::accumulate(R_l0D_hat.begin(),R_l0D_hat.begin()+k,0.0) +std::accumulate(R_l0E_hat.begin(),R_l0E_hat.begin()+k,0.0) +std::accumulate(R_l0F_hat.begin(),R_l0F_hat.begin()+k,0.0) +std::accumulate(R_l0G_hat.begin(),R_l0G_hat.begin()+k,0.0) +std::accumulate(R_l0H_hat.begin(),R_l0H_hat.begin()+k,0.0))/k)); sepruleF = (SSB_hatI[k]+SSB_hatJ[k])/ (SPR_init_F40*((std::accumulate(R_l0I_hat.begin(),R_l0I_hat.begin()+k,0.0) +std::accumulate(R_l0J_hat.begin(),R_l0J_hat.begin()+k,0.0))/k)); if (sepruleA>=1){Fish1 = fishmort;} if (sepruleA<=0.5){Fish1 = 0;} if ((sepruleA<1)&(sepruleA>0.5)) {Fish1 = fishmort*(sepruleA-0.05)*(1/0.95);}//Fish 1 is fishing pressure in whole area 1 if (sepruleF>=1){Fish2 = fishmort;} if (sepruleF<=0.5){Fish2 = 0;} if ((sepruleF<1)&(sepruleF>0.5)) {Fish2 = fishmort*(sepruleF-0.05)*(1/0.95);}//Fish 2 is fishing pressure in whole area 2 temp14 = ((Wt_fem+Wt_mal)/2)*(N_hat_aA+N_hat_aB+N_hat_aC+N_hat_aD+N_hat_aE+N_hat_aF+N_hat_aG+N_hat_aH)*((fishsel*Fish1)/(fishsel*Fish1+M))*(1-exp(-(fishsel*Fish1+M)));//number in left half TAC1[k]=std::accumulate(temp14.begin(),temp14.end(),0.0); temp14 = ((Wt_fem+Wt_mal)/2)*(N_hat_aI+N_hat_aJ)*((fishsel*Fish2)/(fishsel*Fish2+M))*(1-exp(-(fishsel*Fish2+M)));//number in left half TAC2[k]=std::accumulate(temp14.begin(),temp14.end(),0.0); //assessment 2 function below Wt=(Wt_fem+Wt_mal)/2; //establish d which is the effort vector // for (int i=0; i<10;i++){d[i]=(i+1)*100;} temp16a=fishsel*Wt*N_init_a1; Effort[0]=(1/pow(d[0],dpow))*std::accumulate(temp16a.begin(),temp16a.end(),0.0); temp16b=fishsel*Wt*N_init_a2; Effort[1]=(1/pow(d[1],dpow))*std::accumulate(temp16b.begin(),temp16b.end(),0.0); temp16c=fishsel*Wt*N_init_a3; Effort[2]=(1/pow(d[2],dpow))*std::accumulate(temp16c.begin(),temp16c.end(),0.0); temp16d=fishsel*Wt*N_init_a4; Effort[3]=(1/pow(d[3],dpow))*std::accumulate(temp16d.begin(),temp16d.end(),0.0); temp16e=fishsel*Wt*N_init_a5; Effort[4]=(1/pow(d[4],dpow))*std::accumulate(temp16e.begin(),temp16e.end(),0.0); temp16f=fishsel*Wt*N_init_a6; Effort[5]=(1/pow(d[5],dpow))*std::accumulate(temp16f.begin(),temp16f.end(),0.0); temp16g=fishsel*Wt*N_init_a7; Effort[6]=(1/pow(d[6],dpow))*std::accumulate(temp16g.begin(),temp16g.end(),0.0); temp16h=fishsel*Wt*N_init_a8; Effort[7]=(1/pow(d[7],dpow))*std::accumulate(temp16h.begin(),temp16h.end(),0.0); temp16i=fishsel*Wt*N_init_a9; Effort[8]=(1/pow(d[8],dpow))*std::accumulate(temp16i.begin(),temp16i.end(),0.0); temp16j=fishsel*Wt*N_init_a10; Effort[9]=(1/pow(d[9],dpow))*std::accumulate(temp16j.begin(),temp16j.end(),0.0); double Effort_tot1=Effort[0]+Effort[1]+Effort[2]+Effort[3]+Effort[4]+Effort[5] +Effort[6]+Effort[7]; for (int i=0;i<8;i++){ if(Effort_tot1>0){//account for the possibility that all fish will die Effort_norm[i]=Effort[i]/Effort_tot1;} else{Effort_norm[i]=0;}} double Effort_tot2=Effort[8]+Effort[9]; for (int i=8;i<10;i++){ if(Effort_tot2>0){//account for the possibility that all fish will die Effort_norm[i]=Effort[i]/Effort_tot2;} else{Effort_norm[i]=0;}} //REMOVE //Effort_norm[0]=0.125; //Effort_norm[1]=0.125; //Effort_norm[2]=0.125; //Effort_norm[3]=0.125; //Effort_norm[4]=0.125; //Effort_norm[5]=0.125; //Effort_norm[6]=0.125; //Effort_norm[7]=0.125; //Effort_norm[8]=0.5; //Effort_norm[9]=0.5; vec[0]=Effort_norm[0]*TAC1[k];vec[1]=Effort_norm[1]*TAC1[k];vec[2]=Effort_norm[2]*TAC1[k]; vec[3]=Effort_norm[3]*TAC1[k]; vec[4]=Effort_norm[4]*TAC1[k]; vec[5]=Effort_norm[5]*TAC1[k];vec[6]=Effort_norm[6]*TAC1[k];vec[7]=Effort_norm[7]*TAC1[k]; vec[8]=Effort_norm[8]*TAC2[k]; vec[9]=Effort_norm[9]*TAC2[k]; Effort_89(k,_)=Effort_norm; }//end if mgmt_counter==6 TACA[k]=vec[0]; TACB[k]=vec[1]; TACC[k]=vec[2]; TACD[k]=vec[3]; TACE[k]=vec[4]; TACF[k]=vec[5]; TACG[k]=vec[6]; TACH[k]=vec[7]; TACI[k]=vec[8]; TACJ[k]=vec[9]; //calculate FishA-J this is fishing effort that results in desired TACA and TACF //first create vector IntegerVector FSH_int = seq_len(600)-1; NumericVector FSH(600); for (int i=0;i<600;i++){ FSH[i]=FSH_int[i]*0.01;} int best_TAC_A; int best_TAC_B; int best_TAC_C; int best_TAC_D; int best_TAC_E; int best_TAC_F; int best_TAC_G; int best_TAC_H; int best_TAC_I; int best_TAC_J; int high_TAC_A; int high_TAC_B; int high_TAC_C; int high_TAC_D; int high_TAC_E; int high_TAC_F; int high_TAC_G; int high_TAC_H; int high_TAC_I; int high_TAC_J; for (int i=0;i<600;i++){ temp20A = Wt*(N_hat_aA)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestA[i]=std::accumulate(temp20A.begin(),temp20A.end(),0.0); } double highestA = *std::max_element(TACtestA.begin(),TACtestA.end()); for(int i=0; i<600;i++){ if(TACtestA[i]<=TACA[k]){best_TAC_A = i;}if(TACtestA[i]==highestA){high_TAC_A = i;} } for (int i=0;i<600;i++){ temp20B = Wt*(N_hat_aB)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestB[i]=std::accumulate(temp20B.begin(),temp20B.end(),0.0); } double highestB = *std::max_element(TACtestB.begin(),TACtestB.end()); for(int i=0; i<600;i++){ if(TACtestB[i]<=TACB[k]){best_TAC_B = i;}if(TACtestB[i]==highestB){high_TAC_B = i;} } for (int i=0;i<600;i++){ temp20C = Wt*(N_hat_aC)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestC[i]=std::accumulate(temp20C.begin(),temp20C.end(),0.0); } double highestC = *std::max_element(TACtestC.begin(),TACtestC.end()); for(int i=0; i<600;i++){ if(TACtestC[i]<=TACC[k]){best_TAC_C = i;}if(TACtestC[i]==highestC){high_TAC_C = i;} } for (int i=0;i<600;i++){ temp20D = Wt*(N_hat_aD)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestD[i]=std::accumulate(temp20D.begin(),temp20D.end(),0.0); } double highestD = *std::max_element(TACtestD.begin(),TACtestD.end()); for(int i=0; i<600;i++){ if(TACtestD[i]<=TACD[k]){best_TAC_D = i;}if(TACtestD[i]==highestD){high_TAC_D = i;} } for (int i=0;i<600;i++){ temp20E = Wt*(N_hat_aE)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestE[i]=std::accumulate(temp20E.begin(),temp20E.end(),0.0); } double highestE = *std::max_element(TACtestE.begin(),TACtestE.end()); for(int i=0; i<600;i++){ if(TACtestE[i]<=TACE[k]){best_TAC_E = i;}if(TACtestE[i]==highestE){high_TAC_E = i;} } for (int i=0;i<600;i++){ temp20F = Wt*(N_hat_aF)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestF[i]=std::accumulate(temp20F.begin(),temp20F.end(),0.0); } double highestF = *std::max_element(TACtestF.begin(),TACtestF.end()); for(int i=0; i<600;i++){ if(TACtestF[i]<=TACF[k]){best_TAC_F = i;}if(TACtestF[i]==highestF){high_TAC_F = i;} } for (int i=0;i<600;i++){ temp20G = Wt*(N_hat_aG)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestG[i]=std::accumulate(temp20G.begin(),temp20G.end(),0.0); } double highestG = *std::max_element(TACtestG.begin(),TACtestG.end()); for(int i=0; i<600;i++){ if(TACtestG[i]<=TACG[k]){best_TAC_G = i;}if(TACtestG[i]==highestG){high_TAC_G = i;} } for (int i=0;i<600;i++){ temp20H = Wt*(N_hat_aH)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestH[i]=std::accumulate(temp20H.begin(),temp20H.end(),0.0); } double highestH = *std::max_element(TACtestH.begin(),TACtestH.end()); for(int i=0; i<600;i++){ if(TACtestH[i]<=TACH[k]){best_TAC_H = i;}if(TACtestH[i]==highestH){high_TAC_H = i;} } for (int i=0;i<600;i++){ temp20I = Wt*(N_hat_aI)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestI[i]=std::accumulate(temp20I.begin(),temp20I.end(),0.0); } double highestI = *std::max_element(TACtestI.begin(),TACtestI.end()); for(int i=0; i<600;i++){ if(TACtestI[i]<=TACI[k]){best_TAC_I = i;}if(TACtestI[i]==highestI){high_TAC_I = i;} } for (int i=0;i<600;i++){ temp20J = Wt*(N_hat_aJ)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestJ[i]=std::accumulate(temp20J.begin(),temp20J.end(),0.0); } double highestJ = *std::max_element(TACtestJ.begin(),TACtestJ.end()); for(int i=0; i<600;i++){ if(TACtestJ[i]<=TACJ[k]){best_TAC_J = i;}if(TACtestJ[i]==highestJ){high_TAC_J = i;} } if(TACA[k]==0){FishA=0;}else{ if(FSH[best_TAC_A]<=highestA){FishA = FSH[best_TAC_A];} else {FishA = FSH[high_TAC_A];}} if(TACB[k]==0){FishB=0;}else{ if(FSH[best_TAC_B]<=highestB){FishB = FSH[best_TAC_B];} else {FishB = FSH[high_TAC_B];}} if(TACC[k]==0){FishC=0;}else{ if(FSH[best_TAC_C]<=highestC){FishC = FSH[best_TAC_C];} else {FishC = FSH[high_TAC_C];}} if(TACD[k]==0){FishD=0;}else{ if(FSH[best_TAC_D]<=highestD){FishD = FSH[best_TAC_D];} else {FishD = FSH[high_TAC_D];}} if(TACE[k]==0){FishE=0;}else{ if(FSH[best_TAC_E]<=highestE){FishE = FSH[best_TAC_E];} else {FishE = FSH[high_TAC_E];}} if(TACF[k]==0){FishF=0;}else{ if(FSH[best_TAC_F]<=highestF){FishF = FSH[best_TAC_F];} else {FishF = FSH[high_TAC_F];}} if(TACG[k]==0){FishG=0;}else{ if(FSH[best_TAC_G]<=highestG){FishG = FSH[best_TAC_G];} else {FishG = FSH[high_TAC_G];}} if(TACH[k]==0){FishH=0;}else{ if(FSH[best_TAC_H]<=highestH){FishH = FSH[best_TAC_H];} else {FishH = FSH[high_TAC_H];}} if(TACI[k]==0){FishI=0;}else{ if(FSH[best_TAC_I]<=highestI){FishI = FSH[best_TAC_I];} else {FishI = FSH[high_TAC_I];}} if(TACJ[k]==0){FishJ=0;}else{ if(FSH[best_TAC_J]<=highestJ){FishJ = FSH[best_TAC_J];} else {FishJ = FSH[high_TAC_J];}} } // end of 2|3 split //here is catch cascading if(mgmt_counter==5){ seprule = (SSB_hat[k])/(SPR_init_F40*( (std::accumulate(R_l0A_hat.begin(),R_l0A_hat.begin()+k,0.0) +std::accumulate(R_l0B_hat.begin(),R_l0B_hat.begin()+k,0.0) +std::accumulate(R_l0C_hat.begin(),R_l0C_hat.begin()+k,0.0) +std::accumulate(R_l0D_hat.begin(),R_l0D_hat.begin()+k,0.0) +std::accumulate(R_l0E_hat.begin(),R_l0E_hat.begin()+k,0.0) +std::accumulate(R_l0F_hat.begin(),R_l0F_hat.begin()+k,0.0) +std::accumulate(R_l0G_hat.begin(),R_l0G_hat.begin()+k,0.0) +std::accumulate(R_l0H_hat.begin(),R_l0H_hat.begin()+k,0.0) +std::accumulate(R_l0I_hat.begin(),R_l0I_hat.begin()+k,0.0) +std::accumulate(R_l0J_hat.begin(),R_l0J_hat.begin()+k,0.0))/k)); if (seprule>=1){Fish = fishmort;} if (seprule<=0.5){Fish = 0;} if ((seprule<1)&(seprule>0.5)) {Fish = fishmort*(seprule-0.05)*(1/0.95);}//CHANGED HERE seprule_comb_vec[k]=seprule; temp14 = ((Wt_fem+Wt_mal)/2)*(N_init_a)*((fishsel*Fish)/(fishsel*Fish+M))*(1-exp(-(fishsel*Fish+M)));//N_init_a is all areas TAC[k]=std::accumulate(temp14.begin(),temp14.end(),0.0); //Total TAC just based on total number assessed //assessment 2 function below Wt=(Wt_fem+Wt_mal)/2; vec_CC[0]=SSB_hatA_S[k]/SSB_hat_S[k]; vec_CC[1]=SSB_hatB_S[k]/SSB_hat_S[k]; vec_CC[2]=SSB_hatC_S[k]/SSB_hat_S[k]; vec_CC[3]=SSB_hatD_S[k]/SSB_hat_S[k]; vec_CC[4]=SSB_hatE_S[k]/SSB_hat_S[k]; vec_CC[5]=SSB_hatF_S[k]/SSB_hat_S[k]; vec_CC[6]=SSB_hatG_S[k]/SSB_hat_S[k]; vec_CC[7]=SSB_hatH_S[k]/SSB_hat_S[k]; vec_CC[8]=SSB_hatI_S[k]/SSB_hat_S[k]; vec_CC[9]=SSB_hatJ_S[k]/SSB_hat_S[k]; //vec tells you how much fish you are actually going to catch in each area //so find the TAC in each area that gets you that many fish. //TAC in each area is the lower of the TAC calculated by the Effort*totalTAC //OR the exploitable biomass in that area Wt*fishsel*N_init_a1 TACA[k]=vec_CC[0]*TAC[k]; //TAC here is weighted by the estimated biomass but not the distance from fishing port TACB[k]=vec_CC[1]*TAC[k]; TACC[k]=vec_CC[2]*TAC[k]; TACD[k]=vec_CC[3]*TAC[k]; TACE[k]=vec_CC[4]*TAC[k]; TACF[k]=vec_CC[5]*TAC[k]; TACG[k]=vec_CC[6]*TAC[k]; TACH[k]=vec_CC[7]*TAC[k]; TACI[k]=vec_CC[8]*TAC[k]; TACJ[k]=vec_CC[9]*TAC[k]; //calculate FishA through FishJ this is fishing effort that results in desired total TAC //first create vector IntegerVector FSH_int = seq_len(600)-1; NumericVector FSH(600); for (int i=0;i<600;i++){ FSH[i]=FSH_int[i]*0.01;} int best_TAC_A; int best_TAC_B; int best_TAC_C; int best_TAC_D; int best_TAC_E; int best_TAC_F; int best_TAC_G; int best_TAC_H; int best_TAC_I; int best_TAC_J; int high_TAC_A; int high_TAC_B; int high_TAC_C; int high_TAC_D; int high_TAC_E; int high_TAC_F; int high_TAC_G; int high_TAC_H; int high_TAC_I; int high_TAC_J; for (int i=0;i<600;i++){ temp20A = Wt*(N_hat_aA_S)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestA[i]=std::accumulate(temp20A.begin(),temp20A.end(),0.0); } double highestA = *std::max_element(TACtestA.begin(),TACtestA.end()); for(int i=0; i<600;i++){ if(TACtestA[i]<=TACA[k]){best_TAC_A = i;}if(TACtestA[i]==highestA){high_TAC_A = i;} } for (int i=0;i<600;i++){ temp20B = Wt*(N_hat_aB_S)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestB[i]=std::accumulate(temp20B.begin(),temp20B.end(),0.0); } double highestB = *std::max_element(TACtestB.begin(),TACtestB.end()); for(int i=0; i<600;i++){ if(TACtestB[i]<=TACB[k]){best_TAC_B = i;}if(TACtestB[i]==highestB){high_TAC_B = i;} } for (int i=0;i<600;i++){ temp20C = Wt*(N_hat_aC_S)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestC[i]=std::accumulate(temp20C.begin(),temp20C.end(),0.0); } double highestC = *std::max_element(TACtestC.begin(),TACtestC.end()); for(int i=0; i<600;i++){ if(TACtestC[i]<=TACC[k]){best_TAC_C = i;}if(TACtestC[i]==highestC){high_TAC_C = i;} } for (int i=0;i<600;i++){ temp20D = Wt*(N_hat_aD_S)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestD[i]=std::accumulate(temp20D.begin(),temp20D.end(),0.0); } double highestD = *std::max_element(TACtestD.begin(),TACtestD.end()); for(int i=0; i<600;i++){ if(TACtestD[i]<=TACD[k]){best_TAC_D = i;}if(TACtestD[i]==highestD){high_TAC_D = i;} } for (int i=0;i<600;i++){ temp20E = Wt*(N_hat_aE_S)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestE[i]=std::accumulate(temp20E.begin(),temp20E.end(),0.0); } double highestE = *std::max_element(TACtestE.begin(),TACtestE.end()); for(int i=0; i<600;i++){ if(TACtestE[i]<=TACE[k]){best_TAC_E = i;}if(TACtestE[i]==highestE){high_TAC_E = i;} } for (int i=0;i<600;i++){ temp20F = Wt*(N_hat_aF_S)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestF[i]=std::accumulate(temp20F.begin(),temp20F.end(),0.0);//TACtestF is a range of catches under different fishing effort } double highestF = *std::max_element(TACtestF.begin(),TACtestF.end());//highestF is the highest TAC for(int i=0; i<600;i++){ if(TACtestF[i]<=TACF[k]){best_TAC_F = i;}if(TACtestF[i]==highestF){high_TAC_F = i;} } for (int i=0;i<600;i++){ temp20G = Wt*(N_hat_aG_S)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestG[i]=std::accumulate(temp20G.begin(),temp20G.end(),0.0); } double highestG = *std::max_element(TACtestG.begin(),TACtestG.end()); for(int i=0; i<600;i++){ if(TACtestG[i]<=TACG[k]){best_TAC_G = i;}if(TACtestG[i]==highestG){high_TAC_G = i;} } for (int i=0;i<600;i++){ temp20H = Wt*(N_hat_aH_S)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestH[i]=std::accumulate(temp20H.begin(),temp20H.end(),0.0); } double highestH = *std::max_element(TACtestH.begin(),TACtestH.end()); for(int i=0; i<600;i++){ if(TACtestH[i]<=TACH[k]){best_TAC_H = i;}if(TACtestH[i]==highestH){high_TAC_H = i;} } for (int i=0;i<600;i++){ temp20I = Wt*(N_hat_aI_S)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestI[i]=std::accumulate(temp20I.begin(),temp20I.end(),0.0); } double highestI = *std::max_element(TACtestI.begin(),TACtestI.end()); for(int i=0; i<600;i++){ if(TACtestI[i]<=TACI[k]){best_TAC_I = i;}if(TACtestI[i]==highestI){high_TAC_I = i;} } for (int i=0;i<600;i++){ temp20J = Wt*(N_hat_aJ_S)*((fishsel*FSH[i])/(fishsel*FSH[i]+M))*(1-exp(-(fishsel*FSH[i]+M))); TACtestJ[i]=std::accumulate(temp20J.begin(),temp20J.end(),0.0); } double highestJ = *std::max_element(TACtestJ.begin(),TACtestJ.end()); for(int i=0; i<600;i++){ if(TACtestJ[i]<=TACJ[k]){best_TAC_J = i;}if(TACtestJ[i]==highestJ){high_TAC_J = i;} } if(TACA[k]==0){FishA=0;}else{ if(FSH[best_TAC_A]<=highestA){FishA = FSH[best_TAC_A];} else {FishA = FSH[high_TAC_A];}} if(TACB[k]==0){FishB=0;}else{ if(FSH[best_TAC_B]<=highestB){FishB = FSH[best_TAC_B];} else {FishB = FSH[high_TAC_B];}} if(TACC[k]==0){FishC=0;}else{ if(FSH[best_TAC_C]<=highestC){FishC = FSH[best_TAC_C];} else {FishC = FSH[high_TAC_C];}} if(TACD[k]==0){FishD=0;}else{ if(FSH[best_TAC_D]<=highestD){FishD = FSH[best_TAC_D];} else {FishD = FSH[high_TAC_D];}} if(TACE[k]==0){FishE=0;}else{ if(FSH[best_TAC_E]<=highestE){FishE = FSH[best_TAC_E];} else {FishE = FSH[high_TAC_E];}} if(TACF[k]==0){FishF=0;}else{ if(FSH[best_TAC_F]<=highestF){FishF = FSH[best_TAC_F];} else {FishF = FSH[high_TAC_F];}} if(TACG[k]==0){FishG=0;}else{ if(FSH[best_TAC_G]<=highestG){FishG = FSH[best_TAC_G];} else {FishG = FSH[high_TAC_G];}} if(TACH[k]==0){FishH=0;}else{ if(FSH[best_TAC_H]<=highestH){FishH = FSH[best_TAC_H];} else {FishH = FSH[high_TAC_H];}} if(TACI[k]==0){FishI=0;}else{ if(FSH[best_TAC_I]<=highestI){FishI = FSH[best_TAC_I];} else {FishI = FSH[high_TAC_I];}} if(TACJ[k]==0){FishJ=0;}else{ if(FSH[best_TAC_J]<=highestJ){FishJ = FSH[best_TAC_J];} else {FishJ = FSH[high_TAC_J];}} } // end of catch cascading } //end of if k>=10 FishAvec[k]=FishA; FishBvec[k]=FishB; FishCvec[k]=FishC; FishDvec[k]=FishD; FishEvec[k]=FishE; FishFvec[k]=FishF; FishGvec[k]=FishG; FishHvec[k]=FishH; FishIvec[k]=FishI; FishJvec[k]=FishJ; //Fishvec[k]=Fish; might not need this //take fishing pressure in each area and calculate TAC for each area temp15=((Wt_fem+Wt_mal)/2)*N_hat_aA*((fishsel*FishA)/(fishsel*FishA+M))*(1-exp(-(fishsel*FishA+M))); TACA[k]=std::accumulate(temp15.begin(),temp15.end(),0.0); temp16=((Wt_fem+Wt_mal)/2)*N_hat_aB*((fishsel*FishB)/(fishsel*FishB+M))*(1-exp(-(fishsel*FishB+M))); TACB[k]=std::accumulate(temp16.begin(),temp16.end(),0.0); temp16=((Wt_fem+Wt_mal)/2)*N_hat_aC*((fishsel*FishC)/(fishsel*FishC+M))*(1-exp(-(fishsel*FishC+M))); TACC[k]=std::accumulate(temp16.begin(),temp16.end(),0.0); temp16=((Wt_fem+Wt_mal)/2)*N_hat_aD*((fishsel*FishD)/(fishsel*FishD+M))*(1-exp(-(fishsel*FishD+M))); TACD[k]=std::accumulate(temp16.begin(),temp16.end(),0.0); temp16=((Wt_fem+Wt_mal)/2)*N_hat_aE*((fishsel*FishE)/(fishsel*FishE+M))*(1-exp(-(fishsel*FishE+M))); TACE[k]=std::accumulate(temp16.begin(),temp16.end(),0.0); temp16=((Wt_fem+Wt_mal)/2)*N_hat_aF*((fishsel*FishF)/(fishsel*FishF+M))*(1-exp(-(fishsel*FishF+M))); TACF[k]=std::accumulate(temp16.begin(),temp16.end(),0.0); temp16=((Wt_fem+Wt_mal)/2)*N_hat_aG*((fishsel*FishG)/(fishsel*FishG+M))*(1-exp(-(fishsel*FishG+M))); TACG[k]=std::accumulate(temp16.begin(),temp16.end(),0.0); temp16=((Wt_fem+Wt_mal)/2)*N_hat_aH*((fishsel*FishH)/(fishsel*FishH+M))*(1-exp(-(fishsel*FishH+M))); TACH[k]=std::accumulate(temp16.begin(),temp16.end(),0.0); temp16=((Wt_fem+Wt_mal)/2)*N_hat_aI*((fishsel*FishI)/(fishsel*FishI+M))*(1-exp(-(fishsel*FishI+M))); TACI[k]=std::accumulate(temp16.begin(),temp16.end(),0.0); temp16=((Wt_fem+Wt_mal)/2)*N_hat_aJ*((fishsel*FishJ)/(fishsel*FishJ+M))*(1-exp(-(fishsel*FishJ+M))); TACJ[k]=std::accumulate(temp16.begin(),temp16.end(),0.0); TAC[k]=TACA[k]+TACB[k]+TACC[k]+TACD[k]+TACE[k]+TACF[k]+TACG[k]+TACH[k]+TACI[k]+TACJ[k]; //get spawning fish numbers spawnersA = Q*N_init_a1*0.5; spawnersB = Q*N_init_a2*0.5; spawnersC = Q*N_init_a3*0.5; spawnersD = Q*N_init_a4*0.5; spawnersE = Q*N_init_a5*0.5; spawnersF = Q*N_init_a6*0.5; spawnersG = Q*N_init_a7*0.5; spawnersH = Q*N_init_a8*0.5; spawnersI = Q*N_init_a9*0.5; spawnersJ = Q*N_init_a10*0.5; for (int i=0; i<21; i++){ temp=int(spawnersA[i]); if((spawnersA[i]+0.5)>=(int(spawnersA[i])+1)){ temp=int(spawnersA[i])+1; } spawnersA_rounded[i]=temp; temp=int(spawnersB[i]); if((spawnersB[i]+0.5)>=(int(spawnersB[i])+1)){ temp=int(spawnersB[i])+1; } spawnersB_rounded[i]=temp; temp=int(spawnersC[i]); if((spawnersC[i]+0.5)>=(int(spawnersC[i])+1)){ temp=int(spawnersC[i])+1; } spawnersC_rounded[i]=temp; temp=int(spawnersD[i]); if((spawnersD[i]+0.5)>=(int(spawnersD[i])+1)){ temp=int(spawnersD[i])+1; } spawnersD_rounded[i]=temp; temp=int(spawnersE[i]); if((spawnersE[i]+0.5)>=(int(spawnersE[i])+1)){ temp=int(spawnersE[i])+1; } spawnersE_rounded[i]=temp; temp=int(spawnersF[i]); if((spawnersF[i]+0.5)>=(int(spawnersF[i])+1)){ temp=int(spawnersF[i])+1; } spawnersF_rounded[i]=temp; temp=int(spawnersG[i]); if((spawnersG[i]+0.5)>=(int(spawnersG[i])+1)){ temp=int(spawnersG[i])+1; } spawnersG_rounded[i]=temp; temp=int(spawnersH[i]); if((spawnersH[i]+0.5)>=(int(spawnersH[i])+1)){ temp=int(spawnersH[i])+1; } spawnersH_rounded[i]=temp; temp=int(spawnersI[i]); if((spawnersI[i]+0.5)>=(int(spawnersI[i])+1)){ temp=int(spawnersI[i])+1; } spawnersI_rounded[i]=temp; temp=int(spawnersJ[i]); if((spawnersJ[i]+0.5)>=(int(spawnersJ[i])+1)){ temp=int(spawnersJ[i])+1; } spawnersJ_rounded[i]=temp; } //make sure there are not more spawners than actual fish. for(int i=0;i<21;i++){ if((N_init_a1[i]/2)(rnorm(1,0,pow(sigma_jr2,0.5))); //Beverton Holt //N_aA and N_aB are numbers next year*****could double check this part variance may be too high N_aA[0]=((4*h*R_l0A[0]*S_lyA)/(S_l0*(1-h)+S_lyA*(5*h-1))*exp(sigma_rec-(sigma_jr2/2))); //S_l0 is the number of spawners in all pops - all same N_aB[0]=((4*h*R_l0B[0]*S_lyB)/(S_l0*(1-h)+S_lyB*(5*h-1))*exp(sigma_rec-(sigma_jr2/2))); //BH spawner recruit relationship N_aC[0]=((4*h*R_l0C[0]*S_lyC)/(S_l0*(1-h)+S_lyC*(5*h-1))*exp(sigma_rec-(sigma_jr2/2))); //do not add sigma_rec until age 3 N_aD[0]=((4*h*R_l0D[0]*S_lyD)/(S_l0*(1-h)+S_lyD*(5*h-1))*exp(sigma_rec-(sigma_jr2/2))); N_aE[0]=((4*h*R_l0E[0]*S_lyE)/(S_l0*(1-h)+S_lyE*(5*h-1))*exp(sigma_rec-(sigma_jr2/2))); N_aF[0]=((4*h*R_l0F[0]*S_lyF)/(S_l0*(1-h)+S_lyF*(5*h-1))*exp(sigma_rec-(sigma_jr2/2))); N_aG[0]=((4*h*R_l0G[0]*S_lyG)/(S_l0*(1-h)+S_lyG*(5*h-1))*exp(sigma_rec-(sigma_jr2/2))); N_aH[0]=((4*h*R_l0H[0]*S_lyH)/(S_l0*(1-h)+S_lyH*(5*h-1))*exp(sigma_rec-(sigma_jr2/2))); N_aI[0]=((4*h*R_l0I[0]*S_lyI)/(S_l0*(1-h)+S_lyI*(5*h-1))*exp(sigma_rec-(sigma_jr2/2))); N_aJ[0]=((4*h*R_l0J[0]*S_lyJ)/(S_l0*(1-h)+S_lyJ*(5*h-1))*exp(sigma_rec-(sigma_jr2/2))); if(NspawnersA[k]==0){N_aA[0]=0;} if(NspawnersB[k]==0){N_aB[0]=0;} if(NspawnersC[k]==0){N_aC[0]=0;} if(NspawnersD[k]==0){N_aD[0]=0;} if(NspawnersE[k]==0){N_aE[0]=0;} if(NspawnersF[k]==0){N_aF[0]=0;} if(NspawnersG[k]==0){N_aG[0]=0;} if(NspawnersH[k]==0){N_aH[0]=0;} if(NspawnersI[k]==0){N_aI[0]=0;} if(NspawnersJ[k]==0){N_aJ[0]=0;} for(int i=1; i<20; i++){ N_aA[i]=N_init_a1[i-1]*exp(-(fishsel[i-1]*FishA+.34)); //correct because you want the selectivity of the next youngest age N_aB[i]=N_init_a2[i-1]*exp(-(fishsel[i-1]*FishB+.34)); //and fishsel now goes from 21 to 1 N_aC[i]=N_init_a3[i-1]*exp(-(fishsel[i-1]*FishC+.34)); //correct because you want the selectivity of the next youngest age N_aD[i]=N_init_a4[i-1]*exp(-(fishsel[i-1]*FishD+.34)); //and fishsel now goes from 21 to 1 N_aE[i]=N_init_a5[i-1]*exp(-(fishsel[i-1]*FishE+.34)); //correct because you want the selectivity of the next youngest age N_aF[i]=N_init_a6[i-1]*exp(-(fishsel[i-1]*FishF+.34)); //and fishsel now goes from 21 to 1 N_aG[i]=N_init_a7[i-1]*exp(-(fishsel[i-1]*FishG+.34)); //correct because you want the selectivity of the next youngest age N_aH[i]=N_init_a8[i-1]*exp(-(fishsel[i-1]*FishH+.34)); //and fishsel now goes from 21 to 1 N_aI[i]=N_init_a9[i-1]*exp(-(fishsel[i-1]*FishI+.34)); //correct because you want the selectivity of the next youngest age N_aJ[i]=N_init_a10[i-1]*exp(-(fishsel[i-1]*FishJ+.34)); //and fishsel now goes from 21 to 1 } N_aA[20]=N_init_a1[19]*exp(-(fishsel[19]*FishA+M))+N_init_a1[20]*exp(-(fishsel[20]*FishA+M)); //dont forget about plus group N_aB[20]=N_init_a2[19]*exp(-(fishsel[19]*FishB+M))+N_init_a2[20]*exp(-(fishsel[20]*FishB+M)); N_aC[20]=N_init_a3[19]*exp(-(fishsel[19]*FishC+M))+N_init_a3[20]*exp(-(fishsel[20]*FishC+M)); N_aD[20]=N_init_a4[19]*exp(-(fishsel[19]*FishD+M))+N_init_a4[20]*exp(-(fishsel[20]*FishD+M)); N_aE[20]=N_init_a5[19]*exp(-(fishsel[19]*FishE+M))+N_init_a5[20]*exp(-(fishsel[20]*FishE+M)); N_aF[20]=N_init_a6[19]*exp(-(fishsel[19]*FishF+M))+N_init_a6[20]*exp(-(fishsel[20]*FishF+M)); N_aG[20]=N_init_a7[19]*exp(-(fishsel[19]*FishG+M))+N_init_a7[20]*exp(-(fishsel[20]*FishG+M)); N_aH[20]=N_init_a8[19]*exp(-(fishsel[19]*FishH+M))+N_init_a8[20]*exp(-(fishsel[20]*FishH+M)); N_aI[20]=N_init_a9[19]*exp(-(fishsel[19]*FishI+M))+N_init_a9[20]*exp(-(fishsel[20]*FishI+M)); N_aJ[20]=N_init_a10[19]*exp(-(fishsel[19]*FishJ+M))+N_init_a10[20]*exp(-(fishsel[20]*FishJ+M)); //now round N_aA and N_aB to integers for (int i=0; i<21; i++){ temp18=int(N_aA[i]); if((N_aA[i]+0.5)>=(int(N_aA[i])+1)){ temp18=int(N_aA[i])+1; } N_aA[i]=temp18; if (temp18%2!=0){ temp2=as(runif(1,-1,1)); if(temp2>=0){N_aA[i]=N_aA[i]+1;} else{N_aA[i]=N_aA[i]-1; } } } for (int i=0; i<21; i++){ temp19=int(N_aB[i]); if((N_aB[i]+0.5)>=(int(N_aB[i])+1)){ temp19=int(N_aB[i])+1;} N_aB[i]=temp19; if (temp19%2!=0){ temp2=as(runif(1,-1,1)); if(temp2>=0){N_aB[i]=N_aB[i]+1;} else{N_aB[i]=N_aB[i]-1; } } } for (int i=0; i<21; i++){ temp19=int(N_aC[i]); if((N_aC[i]+0.5)>=(int(N_aC[i])+1)){ temp19=int(N_aC[i])+1;} N_aC[i]=temp19; if (temp19%2!=0){ temp2=as(runif(1,-1,1)); if(temp2>=0){N_aC[i]=N_aC[i]+1;} else{N_aC[i]=N_aC[i]-1; }}} for (int i=0; i<21; i++){ temp19=int(N_aD[i]); if((N_aD[i]+0.5)>=(int(N_aD[i])+1)){ temp19=int(N_aD[i])+1;} N_aD[i]=temp19; if (temp19%2!=0){ temp2=as(runif(1,-1,1)); if(temp2>=0){N_aD[i]=N_aD[i]+1;} else{N_aD[i]=N_aD[i]-1; }}} for (int i=0; i<21; i++){ temp19=int(N_aE[i]); if((N_aE[i]+0.5)>=(int(N_aE[i])+1)){ temp19=int(N_aE[i])+1;} N_aE[i]=temp19; if (temp19%2!=0){ temp2=as(runif(1,-1,1)); if(temp2>=0){N_aE[i]=N_aE[i]+1;} else{N_aE[i]=N_aE[i]-1; }}} for (int i=0; i<21; i++){ temp19=int(N_aF[i]); if((N_aF[i]+0.5)>=(int(N_aF[i])+1)){ temp19=int(N_aF[i])+1;} N_aF[i]=temp19; if (temp19%2!=0){ temp2=as(runif(1,-1,1)); if(temp2>=0){N_aF[i]=N_aF[i]+1;} else{N_aF[i]=N_aF[i]-1; }}} for (int i=0; i<21; i++){ temp19=int(N_aG[i]); if((N_aG[i]+0.5)>=(int(N_aG[i])+1)){ temp19=int(N_aG[i])+1;} N_aG[i]=temp19; if (temp19%2!=0){ temp2=as(runif(1,-1,1)); if(temp2>=0){N_aG[i]=N_aG[i]+1;} else{N_aG[i]=N_aG[i]-1; }}} for (int i=0; i<21; i++){ temp19=int(N_aH[i]); if((N_aH[i]+0.5)>=(int(N_aH[i])+1)){ temp19=int(N_aH[i])+1;} N_aH[i]=temp19; if (temp19%2!=0){ temp2=as(runif(1,-1,1)); if(temp2>=0){N_aH[i]=N_aH[i]+1;} else{N_aH[i]=N_aH[i]-1; }}} for (int i=0; i<21; i++){ temp19=int(N_aI[i]); if((N_aI[i]+0.5)>=(int(N_aI[i])+1)){ temp19=int(N_aI[i])+1;} N_aI[i]=temp19; if (temp19%2!=0){ temp2=as(runif(1,-1,1)); if(temp2>=0){N_aI[i]=N_aI[i]+1;} else{N_aI[i]=N_aI[i]-1; }}} for (int i=0; i<21; i++){ temp19=int(N_aJ[i]); if((N_aJ[i]+0.5)>=(int(N_aJ[i])+1)){ temp19=int(N_aJ[i])+1;} N_aJ[i]=temp19; if (temp19%2!=0){ temp2=as(runif(1,-1,1)); if(temp2>=0){N_aJ[i]=N_aJ[i]+1;} else{N_aJ[i]=N_aJ[i]-1; }}} List Nvec(10); //Nvec[0]=N_init_a1;Nvec[1]=N_init_a2;Nvec[2]=N_init_a3;Nvec[3]=N_init_a4;Nvec[4]=N_init_a5; //Nvec[5]=N_init_a6;Nvec[6]=N_init_a7;Nvec[7]=N_init_a8;Nvec[8]=N_init_a9;Nvec[9]=N_init_a10; Nvec[0]=N_aA;Nvec[1]=N_aB;Nvec[2]=N_aC;Nvec[3]=N_aD;Nvec[4]=N_aE; Nvec[5]=N_aF;Nvec[6]=N_aG;Nvec[7]=N_aH;Nvec[8]=N_aI;Nvec[9]=N_aJ; //for each year so do not use this: //if (k%8==zcounter%8&k>1){ for (int i=0;i<10;i++){ //here calculate the number of migrants in a 10x10 matrix only off diagonals are greater than zero for (int j=0;j<10;j++){ if((i+1==j)|(i-1==j)){ if(i+1==j){ Nvec_temp=Nvec[i];} //note all N_init_a1_FIRSTs are the same so it does not matter which goes here if(i-1==j){ Nvec_temp=Nvec[i];} double NMIG=nmig*(1/(std::accumulate(N_init_a1_FIRST.begin(),N_init_a1_FIRST.end(),0.0)))*std::accumulate(Nvec_temp.begin(),Nvec_temp.end(),0.0); temp=int(NMIG); if((temp+0.5)>=(int(NMIG)+1)){ temp=int(NMIG)+1; } int nmig2=temp; if (nmig2%2!=0){ double temp2=as(runif(1,-1,1)); if(temp2>=0){nmig2=nmig2+1;} else{nmig2=nmig2-1;} } migrants(i,j)=nmig2; } }//iloop }//jloop miglist[k]=clone(migrants); N_a_vec[0]=N_aA;N_a_vec[1]=N_aB;N_a_vec[2]=N_aC;N_a_vec[3]=N_aD;N_a_vec[4]=N_aE; N_a_vec[5]=N_aF;N_a_vec[6]=N_aG;N_a_vec[7]=N_aH;N_a_vec[8]=N_aI;N_a_vec[9]=N_aJ; for(int i=0;i<10;i++){ NumericVector N_a=N_a_vec[i]; int migdiff; if((i>0)&(i<9)){ migdiff=10*((migrants(i+1,i)+migrants(i-1,i))-(migrants(i,i-1)+migrants(i,i+1)));} //into minus out of if(i==0){migdiff=10*(migrants(1,0)-migrants(0,1));} //into minus out of if(i==9){migdiff=10*(migrants(8,9)-migrants(9,8));} //into minus out of if (mgmt_counter==0){ migdifflist_comb(k,i)=migdiff;} if (mgmt_counter==1){ migdifflist_sep(k,i)=migdiff;} if (mgmt_counter==2){ migdifflist_56(k,i)=migdiff;} if (mgmt_counter==3){ migdifflist_23(k,i)=migdiff;} if (mgmt_counter==5){ migdifflist_CC(k,i)=migdiff;} if (mgmt_counter==6){ migdifflist_89(k,i)=migdiff;} if (migdiff>0){ int mig_which; for (int i=0;i(runif(1,0,12818)); //works, I checked it int pick3=int(prob); if((prob+0.5)>=(int(prob)+1)){ pick3=int(prob)+1; } mig_which=PROBS[pick3]; N_a[mig_which]=N_a[mig_which]+2; } }//end of if migdiff>0 so add individuals if (migdiff<0){ int mig_which; for (int i=0;i<(-1*migdiff/2);i++){ //figure out which age group to add the new migrants to, make migdiff positive double prob=as(runif(1,0,12818)); //works, I checked it int pick3=int(prob); if((prob+0.5)>=(int(prob)+1)){ pick3=int(prob)+1; } mig_which=PROBS[pick3]; if (N_a[mig_which]>=2){ N_a[mig_which]=N_a[mig_which]-2;}} } N_a_vec[i]=N_a;//put N_a back into the vector to save it }//end of i loop //restore N_aA through N_aJ N_aA=as(N_a_vec[0]);N_aB=as(N_a_vec[1]);N_aC=as(N_a_vec[2]);N_aD=as(N_a_vec[3]);N_aE=as(N_a_vec[4]); N_aF=as(N_a_vec[5]);N_aG=as(N_a_vec[6]);N_aH=as(N_a_vec[7]);N_aI=as(N_a_vec[8]);N_aJ=as(N_a_vec[9]); N_init_a1=N_aA; N_init_a2=N_aB; N_init_a3=N_aC; N_init_a4=N_aD; N_init_a5=N_aE; N_init_a6=N_aF; N_init_a7=N_aG; N_init_a8=N_aH; N_init_a9=N_aI; N_init_a10=N_aJ; //how many individuals are there? R_l0A[k]=N_init_a1[0]; //R_l0A is the number of recruits R_l0B[k]=N_init_a2[0]; R_l0C[k]=N_init_a3[0]; R_l0D[k]=N_init_a4[0]; R_l0E[k]=N_init_a5[0]; R_l0F[k]=N_init_a6[0]; R_l0G[k]=N_init_a7[0]; R_l0H[k]=N_init_a8[0]; R_l0I[k]=N_init_a9[0]; R_l0J[k]=N_init_a10[0]; NtotalA[k] = std::accumulate(N_init_a1.begin(),N_init_a1.end(),0.0); NtotalB[k] = std::accumulate(N_init_a2.begin(),N_init_a2.end(),0.0); NtotalC[k] = std::accumulate(N_init_a3.begin(),N_init_a3.end(),0.0); NtotalD[k] = std::accumulate(N_init_a4.begin(),N_init_a4.end(),0.0); NtotalE[k] = std::accumulate(N_init_a5.begin(),N_init_a5.end(),0.0); NtotalF[k] = std::accumulate(N_init_a6.begin(),N_init_a6.end(),0.0); NtotalG[k] = std::accumulate(N_init_a7.begin(),N_init_a7.end(),0.0); NtotalH[k] = std::accumulate(N_init_a8.begin(),N_init_a8.end(),0.0); NtotalI[k] = std::accumulate(N_init_a9.begin(),N_init_a9.end(),0.0); NtotalJ[k] = std::accumulate(N_init_a10.begin(),N_init_a10.end(),0.0); allA.row(k)=rev(N_init_a1); //each row is the number of each age during that year - reversed so it is 21-1 for next section allB.row(k)=rev(N_init_a2); allC.row(k)=rev(N_init_a3); allD.row(k)=rev(N_init_a4); allE.row(k)=rev(N_init_a5); allF.row(k)=rev(N_init_a6); allG.row(k)=rev(N_init_a7); allH.row(k)=rev(N_init_a8); allI.row(k)=rev(N_init_a9); allJ.row(k)=rev(N_init_a10); allspawnA.row(k)=rev(spawnersA_rounded); //each row is the number of each age during that year - reversed so it is 21-1 for next section allspawnB.row(k)=rev(spawnersB_rounded); allspawnC.row(k)=rev(spawnersC_rounded); allspawnD.row(k)=rev(spawnersD_rounded); allspawnE.row(k)=rev(spawnersE_rounded); allspawnF.row(k)=rev(spawnersF_rounded); allspawnG.row(k)=rev(spawnersG_rounded); allspawnH.row(k)=rev(spawnersH_rounded); allspawnI.row(k)=rev(spawnersI_rounded); allspawnJ.row(k)=rev(spawnersJ_rounded); //create separate matrices for each of the mgmt cases. Can return for separated case (left out here). if(mgmt_counter==0){ if (S_lyA_mat(zcounter,0)==0){SSBrel_init_combA(zcounter,k)=0;}else{SSBrel_init_combA(zcounter,k)=S_lyA/S_lyA_mat(zcounter,0);} //true SSB relative to initial SSB popA if (S_lyB_mat(zcounter,0)==0){SSBrel_init_combB(zcounter,k)=0;}else{SSBrel_init_combB(zcounter,k)=S_lyB/S_lyB_mat(zcounter,0);} if (S_lyC_mat(zcounter,0)==0){SSBrel_init_combC(zcounter,k)=0;}else{SSBrel_init_combC(zcounter,k)=S_lyC/S_lyC_mat(zcounter,0);} if (S_lyD_mat(zcounter,0)==0){SSBrel_init_combD(zcounter,k)=0;}else{SSBrel_init_combD(zcounter,k)=S_lyD/S_lyD_mat(zcounter,0);} if (S_lyE_mat(zcounter,0)==0){SSBrel_init_combE(zcounter,k)=0;}else{SSBrel_init_combE(zcounter,k)=S_lyE/S_lyE_mat(zcounter,0);} if (S_lyF_mat(zcounter,0)==0){SSBrel_init_combF(zcounter,k)=0;}else{SSBrel_init_combF(zcounter,k)=S_lyF/S_lyF_mat(zcounter,0);} if (S_lyG_mat(zcounter,0)==0){SSBrel_init_combG(zcounter,k)=0;}else{SSBrel_init_combG(zcounter,k)=S_lyG/S_lyG_mat(zcounter,0);} if (S_lyH_mat(zcounter,0)==0){SSBrel_init_combH(zcounter,k)=0;}else{SSBrel_init_combH(zcounter,k)=S_lyH/S_lyH_mat(zcounter,0);} if (S_lyI_mat(zcounter,0)==0){SSBrel_init_combI(zcounter,k)=0;}else{SSBrel_init_combI(zcounter,k)=S_lyI/S_lyI_mat(zcounter,0);} if (S_lyJ_mat(zcounter,0)==0){SSBrel_init_combJ(zcounter,k)=0;}else{SSBrel_init_combJ(zcounter,k)=S_lyJ/S_lyJ_mat(zcounter,0);} SSB_hatA_comb(zcounter,k)=SSB_hatA[k]; SSB_hatB_comb(zcounter,k)=SSB_hatB[k]; SSB_hatC_comb(zcounter,k)=SSB_hatC[k]; SSB_hatD_comb(zcounter,k)=SSB_hatD[k]; SSB_hatE_comb(zcounter,k)=SSB_hatE[k]; SSB_hatF_comb(zcounter,k)=SSB_hatF[k]; SSB_hatG_comb(zcounter,k)=SSB_hatG[k]; SSB_hatH_comb(zcounter,k)=SSB_hatH[k]; SSB_hatI_comb(zcounter,k)=SSB_hatI[k]; SSB_hatJ_comb(zcounter,k)=SSB_hatJ[k]; } if(mgmt_counter==1){ //this is totally separated each area different if (S_lyA_mat(zcounter,0)==0){SSBrel_init_sepA(zcounter,k)=0;}else{SSBrel_init_sepA(zcounter,k)=S_lyA/S_lyA_mat(zcounter,0);} //true SSB relative to initial SSB popA if (S_lyB_mat(zcounter,0)==0){SSBrel_init_sepB(zcounter,k)=0;}else{SSBrel_init_sepB(zcounter,k)=S_lyB/S_lyB_mat(zcounter,0);} if (S_lyC_mat(zcounter,0)==0){SSBrel_init_sepC(zcounter,k)=0;}else{SSBrel_init_sepC(zcounter,k)=S_lyC/S_lyC_mat(zcounter,0);} if (S_lyD_mat(zcounter,0)==0){SSBrel_init_sepD(zcounter,k)=0;}else{SSBrel_init_sepD(zcounter,k)=S_lyD/S_lyD_mat(zcounter,0);} if (S_lyE_mat(zcounter,0)==0){SSBrel_init_sepE(zcounter,k)=0;}else{SSBrel_init_sepE(zcounter,k)=S_lyE/S_lyE_mat(zcounter,0);} if (S_lyF_mat(zcounter,0)==0){SSBrel_init_sepF(zcounter,k)=0;}else{SSBrel_init_sepF(zcounter,k)=S_lyF/S_lyF_mat(zcounter,0);} if (S_lyG_mat(zcounter,0)==0){SSBrel_init_sepG(zcounter,k)=0;}else{SSBrel_init_sepG(zcounter,k)=S_lyG/S_lyG_mat(zcounter,0);} if (S_lyH_mat(zcounter,0)==0){SSBrel_init_sepH(zcounter,k)=0;}else{SSBrel_init_sepH(zcounter,k)=S_lyH/S_lyH_mat(zcounter,0);} if (S_lyI_mat(zcounter,0)==0){SSBrel_init_sepI(zcounter,k)=0;}else{SSBrel_init_sepI(zcounter,k)=S_lyI/S_lyI_mat(zcounter,0);} if (S_lyJ_mat(zcounter,0)==0){SSBrel_init_sepJ(zcounter,k)=0;}else{SSBrel_init_sepJ(zcounter,k)=S_lyJ/S_lyJ_mat(zcounter,0);} SSB_hatA_sep(zcounter,k)=SSB_hatA[k]; SSB_hatB_sep(zcounter,k)=SSB_hatB[k]; SSB_hatC_sep(zcounter,k)=SSB_hatC[k]; SSB_hatD_sep(zcounter,k)=SSB_hatD[k]; SSB_hatE_sep(zcounter,k)=SSB_hatE[k]; SSB_hatF_sep(zcounter,k)=SSB_hatF[k]; SSB_hatG_sep(zcounter,k)=SSB_hatG[k]; SSB_hatH_sep(zcounter,k)=SSB_hatH[k]; SSB_hatI_sep(zcounter,k)=SSB_hatI[k]; SSB_hatJ_sep(zcounter,k)=SSB_hatJ[k]; } if(mgmt_counter==2){ if (S_lyA_mat(zcounter,0)==0){SSBrel_init_2unitsA(zcounter,k)=0;}else{SSBrel_init_2unitsA(zcounter,k)=S_lyA/S_lyA_mat(zcounter,0);} //true SSB relative to initial SSB popA if (S_lyB_mat(zcounter,0)==0){SSBrel_init_2unitsB(zcounter,k)=0;}else{SSBrel_init_2unitsB(zcounter,k)=S_lyB/S_lyB_mat(zcounter,0);} if (S_lyC_mat(zcounter,0)==0){SSBrel_init_2unitsC(zcounter,k)=0;}else{SSBrel_init_2unitsC(zcounter,k)=S_lyC/S_lyC_mat(zcounter,0);} if (S_lyD_mat(zcounter,0)==0){SSBrel_init_2unitsD(zcounter,k)=0;}else{SSBrel_init_2unitsD(zcounter,k)=S_lyD/S_lyD_mat(zcounter,0);} if (S_lyE_mat(zcounter,0)==0){SSBrel_init_2unitsE(zcounter,k)=0;}else{SSBrel_init_2unitsE(zcounter,k)=S_lyE/S_lyE_mat(zcounter,0);} if (S_lyF_mat(zcounter,0)==0){SSBrel_init_2unitsF(zcounter,k)=0;}else{SSBrel_init_2unitsF(zcounter,k)=S_lyF/S_lyF_mat(zcounter,0);} if (S_lyG_mat(zcounter,0)==0){SSBrel_init_2unitsG(zcounter,k)=0;}else{SSBrel_init_2unitsG(zcounter,k)=S_lyG/S_lyG_mat(zcounter,0);} if (S_lyH_mat(zcounter,0)==0){SSBrel_init_2unitsH(zcounter,k)=0;}else{SSBrel_init_2unitsH(zcounter,k)=S_lyH/S_lyH_mat(zcounter,0);} if (S_lyI_mat(zcounter,0)==0){SSBrel_init_2unitsI(zcounter,k)=0;}else{SSBrel_init_2unitsI(zcounter,k)=S_lyI/S_lyI_mat(zcounter,0);} if (S_lyJ_mat(zcounter,0)==0){SSBrel_init_2unitsJ(zcounter,k)=0;}else{SSBrel_init_2unitsJ(zcounter,k)=S_lyJ/S_lyJ_mat(zcounter,0);} SSB_hatA_56(zcounter,k)=SSB_hatA[k]; SSB_hatB_56(zcounter,k)=SSB_hatB[k]; SSB_hatC_56(zcounter,k)=SSB_hatC[k]; SSB_hatD_56(zcounter,k)=SSB_hatD[k]; SSB_hatE_56(zcounter,k)=SSB_hatE[k]; SSB_hatF_56(zcounter,k)=SSB_hatF[k]; SSB_hatG_56(zcounter,k)=SSB_hatG[k]; SSB_hatH_56(zcounter,k)=SSB_hatH[k]; SSB_hatI_56(zcounter,k)=SSB_hatI[k]; SSB_hatJ_56(zcounter,k)=SSB_hatJ[k]; } if(mgmt_counter==3){//split between 2 and 3 if (S_lyA_mat(zcounter,0)==0){SSBrel_init_23_A(zcounter,k)=0;}else{SSBrel_init_23_A(zcounter,k)=S_lyA/S_lyA_mat(zcounter,0);} //true SSB relative to initial SSB popA if (S_lyB_mat(zcounter,0)==0){SSBrel_init_23_B(zcounter,k)=0;}else{SSBrel_init_23_B(zcounter,k)=S_lyB/S_lyB_mat(zcounter,0);} if (S_lyC_mat(zcounter,0)==0){SSBrel_init_23_C(zcounter,k)=0;}else{SSBrel_init_23_C(zcounter,k)=S_lyC/S_lyC_mat(zcounter,0);} if (S_lyD_mat(zcounter,0)==0){SSBrel_init_23_D(zcounter,k)=0;}else{SSBrel_init_23_D(zcounter,k)=S_lyD/S_lyD_mat(zcounter,0);} if (S_lyE_mat(zcounter,0)==0){SSBrel_init_23_E(zcounter,k)=0;}else{SSBrel_init_23_E(zcounter,k)=S_lyE/S_lyE_mat(zcounter,0);} if (S_lyF_mat(zcounter,0)==0){SSBrel_init_23_F(zcounter,k)=0;}else{SSBrel_init_23_F(zcounter,k)=S_lyF/S_lyF_mat(zcounter,0);} if (S_lyG_mat(zcounter,0)==0){SSBrel_init_23_G(zcounter,k)=0;}else{SSBrel_init_23_G(zcounter,k)=S_lyG/S_lyG_mat(zcounter,0);} if (S_lyH_mat(zcounter,0)==0){SSBrel_init_23_H(zcounter,k)=0;}else{SSBrel_init_23_H(zcounter,k)=S_lyH/S_lyH_mat(zcounter,0);} if (S_lyI_mat(zcounter,0)==0){SSBrel_init_23_I(zcounter,k)=0;}else{SSBrel_init_23_I(zcounter,k)=S_lyI/S_lyI_mat(zcounter,0);} if (S_lyJ_mat(zcounter,0)==0){SSBrel_init_23_J(zcounter,k)=0;}else{SSBrel_init_23_J(zcounter,k)=S_lyJ/S_lyJ_mat(zcounter,0);} SSB_hatA_23(zcounter,k)=SSB_hatA[k]; SSB_hatB_23(zcounter,k)=SSB_hatB[k]; SSB_hatC_23(zcounter,k)=SSB_hatC[k]; SSB_hatD_23(zcounter,k)=SSB_hatD[k]; SSB_hatE_23(zcounter,k)=SSB_hatE[k]; SSB_hatF_23(zcounter,k)=SSB_hatF[k]; SSB_hatG_23(zcounter,k)=SSB_hatG[k]; SSB_hatH_23(zcounter,k)=SSB_hatH[k]; SSB_hatI_23(zcounter,k)=SSB_hatI[k]; SSB_hatJ_23(zcounter,k)=SSB_hatJ[k]; } if(mgmt_counter==5){//catch cascading if (S_lyA_mat(zcounter,0)==0){SSBrel_init_tierCCA(zcounter,k)=0;}else{SSBrel_init_tierCCA(zcounter,k)=S_lyA/S_lyA_mat(zcounter,0);} //true SSB relative to initial SSB popA if (S_lyB_mat(zcounter,0)==0){SSBrel_init_tierCCB(zcounter,k)=0;}else{SSBrel_init_tierCCB(zcounter,k)=S_lyB/S_lyB_mat(zcounter,0);} if (S_lyC_mat(zcounter,0)==0){SSBrel_init_tierCCC(zcounter,k)=0;}else{SSBrel_init_tierCCC(zcounter,k)=S_lyC/S_lyC_mat(zcounter,0);} if (S_lyD_mat(zcounter,0)==0){SSBrel_init_tierCCD(zcounter,k)=0;}else{SSBrel_init_tierCCD(zcounter,k)=S_lyD/S_lyD_mat(zcounter,0);} if (S_lyE_mat(zcounter,0)==0){SSBrel_init_tierCCE(zcounter,k)=0;}else{SSBrel_init_tierCCE(zcounter,k)=S_lyE/S_lyE_mat(zcounter,0);} if (S_lyF_mat(zcounter,0)==0){SSBrel_init_tierCCF(zcounter,k)=0;}else{SSBrel_init_tierCCF(zcounter,k)=S_lyF/S_lyF_mat(zcounter,0);} if (S_lyG_mat(zcounter,0)==0){SSBrel_init_tierCCG(zcounter,k)=0;}else{SSBrel_init_tierCCG(zcounter,k)=S_lyG/S_lyG_mat(zcounter,0);} if (S_lyH_mat(zcounter,0)==0){SSBrel_init_tierCCH(zcounter,k)=0;}else{SSBrel_init_tierCCH(zcounter,k)=S_lyH/S_lyH_mat(zcounter,0);} if (S_lyI_mat(zcounter,0)==0){SSBrel_init_tierCCI(zcounter,k)=0;}else{SSBrel_init_tierCCI(zcounter,k)=S_lyI/S_lyI_mat(zcounter,0);} if (S_lyJ_mat(zcounter,0)==0){SSBrel_init_tierCCJ(zcounter,k)=0;}else{SSBrel_init_tierCCJ(zcounter,k)=S_lyJ/S_lyJ_mat(zcounter,0);} SSB_hatA_CC(zcounter,k)=SSB_hatA[k]; SSB_hatB_CC(zcounter,k)=SSB_hatB[k]; SSB_hatC_CC(zcounter,k)=SSB_hatC[k]; SSB_hatD_CC(zcounter,k)=SSB_hatD[k]; SSB_hatE_CC(zcounter,k)=SSB_hatE[k]; SSB_hatF_CC(zcounter,k)=SSB_hatF[k]; SSB_hatG_CC(zcounter,k)=SSB_hatG[k]; SSB_hatH_CC(zcounter,k)=SSB_hatH[k]; SSB_hatI_CC(zcounter,k)=SSB_hatI[k]; SSB_hatJ_CC(zcounter,k)=SSB_hatJ[k]; } if(mgmt_counter==6){//split between 8 and 9 if (S_lyA_mat(zcounter,0)==0){SSBrel_init_89A(zcounter,k)=0;}else{SSBrel_init_89A(zcounter,k)=S_lyA/S_lyA_mat(zcounter,0);} //true SSB relative to initial SSB popA if (S_lyB_mat(zcounter,0)==0){SSBrel_init_89B(zcounter,k)=0;}else{SSBrel_init_89B(zcounter,k)=S_lyB/S_lyB_mat(zcounter,0);} if (S_lyC_mat(zcounter,0)==0){SSBrel_init_89C(zcounter,k)=0;}else{SSBrel_init_89C(zcounter,k)=S_lyC/S_lyC_mat(zcounter,0);} if (S_lyD_mat(zcounter,0)==0){SSBrel_init_89D(zcounter,k)=0;}else{SSBrel_init_89D(zcounter,k)=S_lyD/S_lyD_mat(zcounter,0);} if (S_lyE_mat(zcounter,0)==0){SSBrel_init_89E(zcounter,k)=0;}else{SSBrel_init_89E(zcounter,k)=S_lyE/S_lyE_mat(zcounter,0);} if (S_lyF_mat(zcounter,0)==0){SSBrel_init_89F(zcounter,k)=0;}else{SSBrel_init_89F(zcounter,k)=S_lyF/S_lyF_mat(zcounter,0);} if (S_lyG_mat(zcounter,0)==0){SSBrel_init_89G(zcounter,k)=0;}else{SSBrel_init_89G(zcounter,k)=S_lyG/S_lyG_mat(zcounter,0);} if (S_lyH_mat(zcounter,0)==0){SSBrel_init_89H(zcounter,k)=0;}else{SSBrel_init_89H(zcounter,k)=S_lyH/S_lyH_mat(zcounter,0);} if (S_lyI_mat(zcounter,0)==0){SSBrel_init_89I(zcounter,k)=0;}else{SSBrel_init_89I(zcounter,k)=S_lyI/S_lyI_mat(zcounter,0);} if (S_lyJ_mat(zcounter,0)==0){SSBrel_init_89J(zcounter,k)=0;}else{SSBrel_init_89J(zcounter,k)=S_lyJ/S_lyJ_mat(zcounter,0);} SSB_hatA_89(zcounter,k)=SSB_hatA[k]; SSB_hatB_89(zcounter,k)=SSB_hatB[k]; SSB_hatC_89(zcounter,k)=SSB_hatC[k]; SSB_hatD_89(zcounter,k)=SSB_hatD[k]; SSB_hatE_89(zcounter,k)=SSB_hatE[k]; SSB_hatF_89(zcounter,k)=SSB_hatF[k]; SSB_hatG_89(zcounter,k)=SSB_hatG[k]; SSB_hatH_89(zcounter,k)=SSB_hatH[k]; SSB_hatI_89(zcounter,k)=SSB_hatI[k]; SSB_hatJ_89(zcounter,k)=SSB_hatJ[k]; }//89 split } //end of k0) {bxA[i] = spawnersA_rounded_star[i]*N_init_a1_star[0]/(2*NspawnersA[k]);}//div by2 because want half the number of N1s here else{bxA[i]=0;} if (NspawnersB[k]*newnumsB[i]>0) {bxB[i] = spawnersB_rounded_star[i]*N_init_a2_star[0]/(2*NspawnersB[k]);} else{bxB[i]=0;} if (NspawnersC[k]*newnumsC[i]>0) {bxC[i] = spawnersC_rounded_star[i]*N_init_a3_star[0]/(2*NspawnersC[k]);} else{bxC[i]=0;} if (NspawnersD[k]*newnumsD[i]>0) {bxD[i] = spawnersD_rounded_star[i]*N_init_a4_star[0]/(2*NspawnersD[k]);} else{bxD[i]=0;} if (NspawnersE[k]*newnumsE[i]>0) {bxE[i] = spawnersE_rounded_star[i]*N_init_a5_star[0]/(2*NspawnersE[k]);} else{bxE[i]=0;} if (NspawnersF[k]*newnumsF[i]>0) {bxF[i] = spawnersF_rounded_star[i]*N_init_a6_star[0]/(2*NspawnersF[k]);} else{bxF[i]=0;} if (NspawnersG[k]*newnumsG[i]>0) {bxG[i] = spawnersG_rounded_star[i]*N_init_a7_star[0]/(2*NspawnersG[k]);} else{bxG[i]=0;} if (NspawnersH[k]*newnumsH[i]>0) {bxH[i] = spawnersH_rounded_star[i]*N_init_a8_star[0]/(2*NspawnersH[k]);} else{bxH[i]=0;} if (NspawnersI[k]*newnumsI[i]>0) {bxI[i] = spawnersI_rounded_star[i]*N_init_a9_star[0]/(2*NspawnersI[k]);} else{bxI[i]=0;} if (NspawnersJ[k]*newnumsJ[i]>0) {bxJ[i] = spawnersJ_rounded_star[i]*N_init_a10_star[0]/(2*NspawnersJ[k]);} else{bxJ[i]=0;} } //Rcout<<"bxA"<0){ b_primexA[i] = 2*bxA[i]/std::accumulate(bxlxA.begin(),bxlxA.end(),0.0);} else{b_primexA[i] =0;} if(std::accumulate(bxlxB.begin(),bxlxB.end()-1,0.0)>0){ b_primexB[i] = 2*bxB[i]/std::accumulate(bxlxB.begin(),bxlxB.end(),0.0);} else{b_primexB[i]=0;} if(std::accumulate(bxlxC.begin(),bxlxC.end()-1,0.0)>0){ b_primexC[i] = 2*bxC[i]/std::accumulate(bxlxC.begin(),bxlxC.end(),0.0);} else{b_primexC[i]=0;} if(std::accumulate(bxlxD.begin(),bxlxD.end()-1,0.0)>0){ b_primexD[i] = 2*bxD[i]/std::accumulate(bxlxD.begin(),bxlxD.end(),0.0);} else{b_primexD[i]=0;} if(std::accumulate(bxlxE.begin(),bxlxE.end()-1,0.0)>0){ b_primexE[i] = 2*bxE[i]/std::accumulate(bxlxE.begin(),bxlxE.end(),0.0);} else{b_primexE[i]=0;} if(std::accumulate(bxlxF.begin(),bxlxF.end()-1,0.0)>0){ b_primexF[i] = 2*bxF[i]/std::accumulate(bxlxF.begin(),bxlxF.end(),0.0);} else{b_primexF[i]=0;} if(std::accumulate(bxlxG.begin(),bxlxG.end()-1,0.0)>0){ b_primexG[i] = 2*bxG[i]/std::accumulate(bxlxG.begin(),bxlxG.end(),0.0);} else{b_primexG[i]=0;} if(std::accumulate(bxlxH.begin(),bxlxH.end()-1,0.0)>0){ b_primexH[i] = 2*bxH[i]/std::accumulate(bxlxH.begin(),bxlxH.end(),0.0);} else{b_primexH[i]=0;} if(std::accumulate(bxlxI.begin(),bxlxI.end()-1,0.0)>0){ b_primexI[i] = 2*bxI[i]/std::accumulate(bxlxI.begin(),bxlxI.end(),0.0);} else{b_primexI[i]=0;} if(std::accumulate(bxlxJ.begin(),bxlxJ.end()-1,0.0)>0){ b_primexJ[i] = 2*bxJ[i]/std::accumulate(bxlxJ.begin(),bxlxJ.end(),0.0);} else{b_primexJ[i]=0;} } BxA=b_primexA*newnumsA; BxB=b_primexB*newnumsB; BxC=b_primexC*newnumsC; BxD=b_primexD*newnumsD; BxE=b_primexE*newnumsE; BxF=b_primexF*newnumsF; BxG=b_primexG*newnumsG; BxH=b_primexH*newnumsH; BxI=b_primexI*newnumsI; BxJ=b_primexJ*newnumsJ; L_popAvec = x_age*BxA; L_popBvec = x_age*BxB; L_popCvec = x_age*BxC; L_popDvec = x_age*BxD; L_popEvec = x_age*BxE; L_popFvec = x_age*BxF; L_popGvec = x_age*BxG; L_popHvec = x_age*BxH; L_popIvec = x_age*BxI; L_popJvec = x_age*BxJ; //Rcout<<"L_popAvec"<0){ L_popA = std::accumulate(L_popAvec.begin(),L_popAvec.end()-1,0.0)/N_init_a1_star[0];} else {L_popA=0;} if (N_init_a2_star[0]>0){ L_popB = std::accumulate(L_popBvec.begin(),L_popBvec.end()-1,0.0)/N_init_a2_star[0];} else{L_popB=0;} if (N_init_a3_star[0]>0){ L_popC = std::accumulate(L_popCvec.begin(),L_popCvec.end()-1,0.0)/N_init_a3_star[0];} else{L_popC=0;} if (N_init_a4_star[0]>0){ L_popD = std::accumulate(L_popDvec.begin(),L_popDvec.end()-1,0.0)/N_init_a4_star[0];} else{L_popD=0;} if (N_init_a5_star[0]>0){ L_popE = std::accumulate(L_popEvec.begin(),L_popEvec.end()-1,0.0)/N_init_a5_star[0];} else{L_popE=0;} if (N_init_a6_star[0]>0){ L_popF = std::accumulate(L_popFvec.begin(),L_popFvec.end()-1,0.0)/N_init_a6_star[0];} else{L_popF=0;} if (N_init_a7_star[0]>0){ L_popG = std::accumulate(L_popGvec.begin(),L_popGvec.end()-1,0.0)/N_init_a7_star[0];} else{L_popG=0;} if (N_init_a8_star[0]>0){ L_popH = std::accumulate(L_popHvec.begin(),L_popHvec.end()-1,0.0)/N_init_a8_star[0];} else{L_popH=0;} if (N_init_a9_star[0]>0){ L_popI = std::accumulate(L_popIvec.begin(),L_popIvec.end()-1,0.0)/N_init_a9_star[0];} else{L_popI=0;} if (N_init_a10_star[0]>0){ L_popJ = std::accumulate(L_popJvec.begin(),L_popJvec.end()-1,0.0)/N_init_a10_star[0];} else{L_popJ=0;} k_barA[0]=0;k_barB[0]=0;k_barC[0]=0;k_barD[0]=0;k_barE[0]=0; k_barF[0]=0;k_barG[0]=0;k_barH[0]=0;k_barI[0]=0;k_barJ[0]=0; for (int i=1;i<21;i++){ k_barA[i] = b_primexA[i]+k_barA[i-1]; k_barB[i] = b_primexB[i]+k_barB[i-1]; k_barC[i] = b_primexC[i]+k_barC[i-1]; k_barD[i] = b_primexD[i]+k_barD[i-1]; k_barE[i] = b_primexE[i]+k_barE[i-1]; k_barF[i] = b_primexF[i]+k_barF[i-1]; k_barG[i] = b_primexG[i]+k_barG[i-1]; k_barH[i] = b_primexH[i]+k_barH[i-1]; k_barI[i] = b_primexI[i]+k_barI[i-1]; k_barJ[i] = b_primexJ[i]+k_barJ[i-1]; } //Rcout<<"k_barA"<=0){DxA[i]=newnumsA[i]-newnumsA[i+1];} // else{DxA[i]=0;} // if (newnumsB[i]-newnumsB[i+1]>=0){DxB[i]=newnumsB[i]-newnumsB[i+1];} // else{DxB[i]=0;} // if (newnumsC[i]-newnumsC[i+1]>=0){DxC[i]=newnumsC[i]-newnumsC[i+1];} // else{DxC[i]=0;} // if (newnumsD[i]-newnumsD[i+1]>=0){DxD[i]=newnumsD[i]-newnumsD[i+1];} // else{DxD[i]=0;} // if (newnumsE[i]-newnumsE[i+1]>=0){DxE[i]=newnumsE[i]-newnumsE[i+1];} // else{DxE[i]=0;} // if (newnumsF[i]-newnumsF[i+1]>=0){DxF[i]=newnumsF[i]-newnumsF[i+1];} // else{DxF[i]=0;} // if (newnumsG[i]-newnumsG[i+1]>=0){DxG[i]=newnumsG[i]-newnumsG[i+1];} // else{DxG[i]=0;} // if (newnumsH[i]-newnumsH[i+1]>=0){DxH[i]=newnumsH[i]-newnumsH[i+1];} // else{DxH[i]=0;} // if (newnumsI[i]-newnumsI[i+1]>=0){DxI[i]=newnumsI[i]-newnumsI[i+1];} // else{DxI[i]=0;} // if (newnumsJ[i]-newnumsJ[i+1]>=0){DxJ[i]=newnumsJ[i]-newnumsJ[i+1];} // else{DxJ[i]=0;} // } k_barADxA=k_barA*DxA; k_barBDxB=k_barB*DxB; k_barCDxC=k_barC*DxC; k_barDDxD=k_barD*DxD; k_barEDxE=k_barE*DxE; k_barFDxF=k_barF*DxF; k_barGDxG=k_barG*DxG; k_barHDxH=k_barH*DxH; k_barIDxI=k_barI*DxI; k_barJDxJ=k_barJ*DxJ; //Rcout<<"k_barADxA"<0){ k_bardotA = 2*std::accumulate(k_barADxA.begin(),k_barADxA.end(),0.0)/N_init_a1_star[0];} else{k_bardotA=0;}//divide by number of newborns if (N_init_a2_star[0]>0){ k_bardotB = 2*std::accumulate(k_barBDxB.begin(),k_barBDxB.end(),0.0)/N_init_a2_star[0];} else{k_bardotB=0;} if (N_init_a3_star[0]>0){ k_bardotC = 2*std::accumulate(k_barCDxC.begin(),k_barCDxC.end(),0.0)/N_init_a3_star[0];} else{k_bardotC=0;} if (N_init_a4_star[0]>0){ k_bardotD = 2*std::accumulate(k_barDDxD.begin(),k_barDDxD.end(),0.0)/N_init_a4_star[0];} else{k_bardotD=0;} if (N_init_a5_star[0]>0){ k_bardotE = 2*std::accumulate(k_barEDxE.begin(),k_barEDxE.end(),0.0)/N_init_a5_star[0];} else{k_bardotE=0;} if (N_init_a6_star[0]>0){ k_bardotF = 2*std::accumulate(k_barFDxF.begin(),k_barFDxF.end(),0.0)/N_init_a6_star[0];} else{k_bardotF=0;} if (N_init_a7_star[0]>0){ k_bardotG = 2*std::accumulate(k_barGDxG.begin(),k_barGDxG.end(),0.0)/N_init_a7_star[0];} else{k_bardotG=0;} if (N_init_a8_star[0]>0){ k_bardotH = 2*std::accumulate(k_barHDxH.begin(),k_barHDxH.end(),0.0)/N_init_a8_star[0];} else{k_bardotH=0;} if (N_init_a9_star[0]>0){ k_bardotI = 2*std::accumulate(k_barIDxI.begin(),k_barIDxI.end(),0.0)/N_init_a9_star[0];} else{k_bardotI=0;} if (N_init_a10_star[0]>0){ k_bardotJ = 2*std::accumulate(k_barJDxJ.begin(),k_barJDxJ.end(),0.0)/N_init_a10_star[0];} else{k_bardotJ=0;} delxA=k_barA-k_bardotA; delxB=k_barB-k_bardotB; delxC=k_barC-k_bardotC; delxD=k_barD-k_bardotD; delxE=k_barE-k_bardotE; delxF=k_barF-k_bardotF; delxG=k_barG-k_bardotG; delxH=k_barH-k_bardotH; delxI=k_barI-k_bardotI; delxJ=k_barJ-k_bardotJ; //Rcout<<"k_bardotA"<0){VkA = SSDtA/N_init_a1_star[0];} //just for females or males but the values are the same else{VkA=0;} if (N_init_a2_star[0]>0){VkB = SSDtB/N_init_a2_star[0];} else{VkB=0;} if (N_init_a3_star[0]>0){VkC = SSDtC/N_init_a3_star[0];} else{VkC=0;} if (N_init_a4_star[0]>0){VkD = SSDtD/N_init_a4_star[0];} else{VkD=0;} if (N_init_a5_star[0]>0){VkE = SSDtE/N_init_a5_star[0];} else{VkE=0;} if (N_init_a6_star[0]>0){VkF = SSDtF/N_init_a6_star[0];} else{VkF=0;} if (N_init_a7_star[0]>0){VkG = SSDtG/N_init_a7_star[0];} else{VkG=0;} if (N_init_a8_star[0]>0){VkH = SSDtH/N_init_a8_star[0];} else{VkH=0;} if (N_init_a9_star[0]>0){VkI = SSDtI/N_init_a9_star[0];} else{VkI=0;} if (N_init_a10_star[0]>0){VkJ = SSDtJ/N_init_a10_star[0];} else{VkJ=0;} //VkA and VkB are total variance //Rcout<<"L_popA"<