//Chesapeake Bay sex-specific catch multiple survye model
//for blue crabs in Chesapeake Bay
//By: M. Wilberg
//Date:10/21/2010
//Last updated: 2/28/2011
//Winter dredge survey and VIMs trawl surveys are end of year surveys
TOP_OF_MAIN_SECTION
  //increase number of estimated parameters
  gradient_structure::set_NUM_DEPENDENT_VARIABLES(1000);
  gradient_structure::set_GRADSTACK_BUFFER_SIZE(200040);
  gradient_structure::set_CMPDIF_BUFFER_SIZE(1000000);
  arrmblsize = 10000000;
  
  //gradient_structure::set_ARRAY_MEMBLOCK_SIZE(150352);
GLOBALS_SECTION

DATA_SECTION  //Read in data
  init_int fyear   //first year of the model
  init_int lyear   //last year of the model

  //Catch data
  init_int ftcyear  //first year of total catch
  init_int ltcyear  //last year of total catch
  init_vector com_TC_obs(ftcyear,ltcyear)  //total catch
  init_int fscyear  //first year sex-specific catch
  init_int lscyear  //last year sex-specific catch
  init_matrix com_C_obs(1,2,fscyear,lscyear)  //Sex-specific catch
  init_number C_sds

  //Survey data
  init_int nsurveys  //number of surveys
  init_ivector fsyear(1,nsurveys)  //vector of first year of surveys
  init_ivector lsyear(1,nsurveys)  //vector of last survey years
  init_3darray ad_survey_obs(1,nsurveys,1,2,fsyear,lsyear)  //Adult survey CPUE
  init_matrix re_survey_obs(1,nsurveys,fsyear,lsyear)  //Recruit survey CPUE
  init_vector sa_be(1,nsurveys)    //indicate beginning or end of year for surveys
  init_vector sa_time(1,nsurveys)  //Timing of surveys
  init_vector sr_time(1,nsurveys)  //timing of surveys
  init_vector ad_survey_sds(1,nsurveys)  //Survey SDs for adults
  init_vector re_survey_sds(1,nsurveys)  //survey SDs for recruits
  init_ivector survey_dist(1,nsurveys)   //Distribution for survey likelihoods
  init_vector in_M(1,2)     //starting value for M
  init_int M_phase     //estimation phase for male M
  init_number in_sex_r    //sex ratio of recruits
  init_int sex_r_phase //phase for sex ratio estimation
  init_number salpha   //starting value for SR alpha
  init_number sbeta    //starting value for SR beta
  init_number in_rec_sd   //standard deviation for recruitment
  init_int rsd_phase   //phase for estimation of recruitment sd
  init_number sp_time  //proportion of the year before spawning occurs
  init_number p_F_sp   //proportion of female recruits that spawn at age-1

  init_number in_rf  //partial recruitment of recruits to the fishery
  init_int rf_phase  //phase for estimation of partial recruitment
  init_number p_rec  //Proportion of recreational harvest

  //number for reference point explorations
  init_number Frat_init  //lowest value for ratio of Male to female F
  init_number Frat_max   //highest value for ratio of male to female F
  init_number Frat_inc   //Increment between values
  
  init_number FSPR_init  //lowest value of female F used in SPR calcs
  init_number FSPR_max   //highest value of female F used in SPR calcs
  init_number FSPR_inc   //increment for F

  int Frat_num
  int FSPR_num

  !!Frat_num=(Frat_max-Frat_init)/Frat_inc+1;
  !!FSPR_num=(FSPR_max-FSPR_init)/FSPR_inc+1;

  init_int nspr          //number of values FSPR will be calculated at
  init_vector SPR_targ(1,nspr)  //Values of SPR for FSPR reference point calculations
  init_int niter         //number of iterations for Newton method of calculating FSPR
  
  init_int test

  vector nyrs(1,nsurveys)  //vector for number of years in each survey
  ivector fy(1,nsurveys)   //vector of alternative first years for survey

  //Total harvest including recreational
  init_vector TC_obs(ftcyear,ltcyear)  //total catch
  init_matrix C_obs(1,2,fscyear,lscyear)  //Sex-specific catch

  //Variances for data sets
  number C_var  //variance of catch
  vector ad_survey_var(1,nsurveys)  //variances for adult surveys
  vector re_survey_var(1,nsurveys)  //variances for recruitment surveys

  //Define index varaibles
  int x  //index variable for sex (1 female, 2 male)
  int y  //index variable for year
  int s  //index variable for survey
  int i  //index variable for ratio of male to female F
  int j  //index variable for F in SPR and MSY calculations

  int ispr
  int iter

 LOCAL_CALCS
  if(test!=12345)  //check to make sure end of file number is correct
  {
    //if not correct, output the data and exit.
    cout << "Data not reading properly" << endl;
    cout << "fyear, lyear: " << fyear << ", " << lyear << endl;
    cout << "ftcyear, ltcyear: " << ftcyear << ", " << ltcyear << endl;
    cout << "Total catch" << endl << TC_obs << endl;
    cout << "fscyear, lscyear: " << ftcyear << ", " << ltcyear << endl;
    cout << "Sex-specific catch" << endl << C_obs << endl;
    cout << "SD of catch: " << C_sds << endl;
    cout << "number of surveys: " << nsurveys << endl;
    cout << "fsyear, lsyear: " << fsyear << ", " << lsyear << endl;
    cout << "Adult survey indices" << endl << ad_survey_obs << endl;
    cout << "Recruit survey indices" << endl << re_survey_obs << endl;
    cout << "M: " << in_M << endl;
    cout << "sex ratio (proportion female): " << in_sex_r << endl;
    cout << "starting values for alpha, beta: " << salpha << ", " << sbeta << endl;
    cout << "EOF test: " << test << endl;
    exit(1);
  }

  //Calculate SDs from variances
  C_var=square(C_sds);  //variance of catch
  ad_survey_var=square(ad_survey_sds);  //variances for adult surveys
  re_survey_var=square(re_survey_sds);  //variances for recruitment surveys

  TC_obs=com_TC_obs*(1.+p_rec);  //total catch
  C_obs=com_C_obs*(1.+p_rec);    //Sex-specific catch

 END_CALCS

PARAMETER_SECTION
  //initial R and N
  init_bounded_number log_init_N(0.,20.,1)
  init_bounded_number log_init_R(0.,20.,1)

  //Fishing mortality for each year
  init_bounded_vector log_mean_F(1,2,-10.,10.,1)
  init_bounded_dev_vector log_F_dev_1(fyear,lyear-1,-10.,10.,2)
  init_bounded_dev_vector log_F_dev_2(fyear,lyear-1,-10.,10.,2)

  //Stock-recruitment parameters
  init_bounded_number log_alpha(-10.,10.,1)
  init_bounded_number log_beta(-10.,10.,4)
  init_bounded_dev_vector log_rec_devs(fyear,lyear-1,-10,10,1)
  init_bounded_number log_rec_sd(-5,3,rsd_phase)
  
  //Natural mortality rate
  init_bounded_number log_M(-10.,2.,M_phase)

  //Partial recruitment (i.e., selectivity) of recruits to the fishery
  init_bounded_number log_rf(-10.,0.,rf_phase)

  init_bounded_number log_sex_r(-10,0,sex_r_phase)    //sex ratio of recruits

  //Calculated values
  matrix N(1,2,fyear,lyear)   //sex-specific abundance
  matrix R(1,2,fyear,lyear)   //sex-specific recruitment
  matrix SP(1,2,fyear,lyear)      //number of spawners
  matrix C(1,2,fyear,lyear-1)   //sex-specific catch
  vector TC(fyear,lyear-1)      //total catch
  matrix F(1,2,fyear,lyear-1)   //sex-specific fishing mortality rate
  3darray ad_survey_est(1,nsurveys,1,2,fyear,lyear)  //estimated adult survey indices
  matrix re_survey_est(1,nsurveys,fyear,lyear)  //estimated recruitment survey indices
  sdreport_matrix qa(1,nsurveys,1,2)   //catchability for adult surveys (sex-specific)
  sdreport_vector qr(1,nsurveys)   //catchability for recruitment surveys
  vector M(1,2)                    //Instantaneous natural mortality rate
  number rf                //selectivity (partial recruitment) of recruits to the fishery
    
  number alpha                //Alpha of the S-R relationship
  number beta                 //Beta of the S-R relationship
  number rec_sd               //standard deviation for recruitment
  number rec_var                    //variance for recruitment deviations
  number sex_r    //sex ratio of recruits

  //Variables for reference point calculations 
  number FSPR  //F for SPR calculations
  number Frat  //Ratio of male to female F

  matrix SPRf(1,Frat_num,1,FSPR_num)
  matrix SPRm(1,Frat_num,1,FSPR_num)

  matrix YPRf(1,Frat_num,1,FSPR_num)
  matrix YPRm(1,Frat_num,1,FSPR_num)

  matrix N_eq(1,Frat_num,1,FSPR_num)
  matrix Nf_eq(1,Frat_num,1,FSPR_num)
  matrix Nm_eq(1,Frat_num,1,FSPR_num)
  matrix Nfs_eq(1,Frat_num,1,FSPR_num)
  matrix Nms_eq(1,Frat_num,1,FSPR_num)

  matrix R_eq(1,Frat_num,1,FSPR_num)
  matrix Rf_eq(1,Frat_num,1,FSPR_num)
  matrix Rm_eq(1,Frat_num,1,FSPR_num)

  matrix C_eq(1,Frat_num,1,FSPR_num)
  matrix Cf_eq(1,Frat_num,1,FSPR_num)
  matrix Cm_eq(1,Frat_num,1,FSPR_num)

  matrix uf_eq(1,Frat_num,1,FSPR_num)
  matrix uf0_eq(1,Frat_num,1,FSPR_num)
  matrix um_eq(1,Frat_num,1,FSPR_num)
  matrix urat_eq(1,Frat_num,1,FSPR_num)

  vector FSPR_ref(1,nspr)
  number SPR_temp
  number dSPR


  //variables for likelihood function
  matrix Lsa(1,nsurveys,1,2)  //likelihood components for adult surveys
  vector Lsr(1,nsurveys)  //likelihood components for recruit surveys
  vector Lc(1,3)         //likelihood components for catch time series (1=female, 2-male, 3-total)
  number Lrdev           //likelihood for recruitment deviations
  number F_pen           //penalty for deviations between male and female F

  number F_pen_mu        //mean for F_pen
  number F_pen_var       //var for F_pen

  //Negative log likelihood
  objective_function_value negLL

 LOCAL_CALCS
  //specify starting values for parameters
  log_alpha=log(salpha);
  log_beta=log(sbeta);
  log_mean_F=log(0.6);
  log_init_N=log(73.7);
  log_init_R=log(150.);
  log_M=log(in_M(2));
  log_rf=log(in_rf);
  log_rec_sd=log(in_rec_sd);
  log_sex_r=log(in_sex_r);

  //calculate number of years in each survey (needs to be here because 
  //functions are not visible in the DATA_SECTION LOCAL_CALCS)
  for(s=1;s<=nsurveys;s++)
  {
    if(fsyear(s)>fyear)
    {
      fy(s)=fsyear(s);
    }
    else
    {
      fy(s)=fyear;
    }
    nyrs(s)=double(lyear-fy(s)+1);
  }

 END_CALCS

PROCEDURE_SECTION
  //Put in values for the first year and convert parameters to the arithmetic scale
  set_initial_conditions();  
  //calculate abundance for each year
  calculate_abundance_and_catch();
  calculate_predicted_indices();
  calculate_objective_function();

FUNCTION set_initial_conditions
  //convert parameters from the log scale
  alpha=exp(log_alpha);
  beta=exp(log_beta);
  rec_sd=exp(log_rec_sd);
  rec_var=square(rec_sd);
  sex_r=exp(log_sex_r);
  M(1)=in_M(1);
  M(2)=exp(log_M);
  rf=exp(log_rf);

  //Calculate F for each year
  for(y=fyear;y<lyear;y++)
  {
    F(1,y)=exp(log_mean_F(1)+log_F_dev_1(y));
    F(2,y)=exp(log_mean_F(2)+log_F_dev_2(y));
  }

  //assign starting abundance and recruitment
  for(x=1;x<=2;x++)
  {
    N(x,fyear)=exp(log_init_N);
    R(x,fyear)=exp(log_init_R);
  }

FUNCTION calculate_abundance_and_catch  
  for(y=fyear;y<lyear;y++)
  {
    for(x=1;x<=2;x++)  //loop over sexes
    {
      //abundance for the next year
      N(x,y+1)=R(x,y)*exp(-(M(x)+rf*F(x,y)))+N(x,y)*exp(-(M(x)+F(x,y)));
      //spawners also include some animals that were recruits in the beginning of the year
      SP(x,y)=nspawn(F(x,y),M(x),N(x,y),R(x,y));

      //Baranov catch equation
      C(x,y)=baranov(F(x,y),M(x),N(x,y))+baranov(rf*F(x,y),M(x),R(x,y));
    }

    //Calculate next year's recruitment
    //Ricker S-R function uses only females for egg production but total abundance for density dependence
    R(1,y+1)=sex_r*alpha*SP(1,y)*exp(-beta*(SP(1,y)+SP(2,y)))*exp(log_rec_devs(y)); 
    R(2,y+1)=R(1,y+1)*(1.-sex_r)/sex_r;
  }

  //Total catch
  for(y=fyear;y<lyear;y++)
  {
    TC(y)=C(1,y)+C(2,y);
  }

FUNCTION calculate_predicted_indices
  for(s=1;s<=nsurveys;s++)
  {
    re_survey_est(s)=0.0;
    qr(s)=0.0;
    double missy=0.0;
    for(y=fy(s);y<=lyear;y++)
    {
      if(re_survey_obs(s,y)!=-99.)  //check to make sure year is not missing
      {
        if(!last_phase())
        {
         //small constant added to recruitment in earlier stages to 
         //increase numerical stability
         qr(s)+=log(re_survey_obs(s,y))-log((R(1,y)+.01)/sex_r);
        }
        else
        {
          //small constant not included in last estimation stage
          qr(s)+=log(re_survey_obs(s,y))-log(R(1,y)/sex_r);
        } 
     }
      else
      {
        //survey index missed a year
        missy+=1.;
      }
    }
    //calculate geometric mean
    qr(s)=exp(qr(s)/(nyrs(s)-missy));

    for(x=1;x<=2;x++)
    {
      //calculate catchability for each sex-index combination
      double missy=0.0;
      qa(s,x)=0.;
      for(y=fy(s);y<=lyear;y++)
      {
      if(s<nsurveys)  
      {
        if(ad_survey_obs(s,x,y)!=-99.)  //check to make sure year is not missing
        {
          if(sa_be(s)==0) //beginning of year survey
          {
            qa(s,x)+=log(ad_survey_obs(s,x,y))-log(N(x,y));
          }
          else
          {
            qa(s,x)+=log(ad_survey_obs(s,x,y))-log(N(x,y)*exp(-sa_time(s)*(M(x)+F(x,y)))+R(x,y)*exp(-sa_time(s)*(M(x)+rf*F(x,y))));
           }
        }
        else
        {
          //survey index missed a year
          missy+=1.;
        }
      }
      else  //last adult survey is the winter dredge survey
      {
        qa(s,x)=1.0;  //Set adult q for winter dredge survey to 1.
      }
      }
      //calculate geometric mean
      if(s<nsurveys) qa(s,x)=exp(qa(s,x)/(nyrs(s)-missy));      

      //Calculate each predicted index of abundance
        if(sa_be(s)==0) //beginning of year survey
        {
          ad_survey_est(s,x)=qa(s,x)*N(x);
        }
        else
        {
          ad_survey_est(s,x)(fy(s),lyear-1)=qa(s,x)*(elem_prod(N(x)(fy(s),lyear-1),exp(-sa_time(s)*(M(x)+F(x)(fy(s),lyear-1))))+
                              elem_prod(R(x)(fy(s),lyear-1),exp(-sa_time(s)*(M(x)+rf*F(x)(fy(s),lyear-1)))));
        }
    }
      re_survey_est(s)=qr(s)*R(1)/sex_r;
  }

FUNCTION calculate_objective_function
  //Calculate survey likelihood components
  for(s=1;s<=nsurveys;s++)
  {
    for(x=1;x<=2;x++)
    {
      //Adult survey
      if(survey_dist(s)==1)
      {
        Lsa(s,x)=lognorm(ad_survey_obs(s,x)(fy(s),lyear),ad_survey_est(s,x)(fy(s),lyear),ad_survey_var(s));
      }
      if(survey_dist(s)==0)
      {
        Lsa(s,x)=normal(ad_survey_obs(s,x)(fy(s),lyear),ad_survey_est(s,x)(fy(s),lyear),ad_survey_sds(s));
      }
    }
    //Recruit survey
    if(survey_dist(s)==1)
    {
      Lsr(s)=lognorm(re_survey_obs(s)(fy(s),lyear),re_survey_est(s)(fy(s),lyear),re_survey_var(s));
    }
    if(survey_dist(s)==0)
    {    
      Lsr(s)=normal(re_survey_obs(s)(fy(s),lyear),re_survey_est(s)(fy(s),lyear),re_survey_sds(s));
    }
  }
  
  //calculate sex-specific catch likelihood components
  for(x=1;x<=2;x++)
  {
    Lc(x)=lognorm(C_obs(x)(fscyear,lyear-1),C(x)(fscyear,lyear-1),C_var);
  }

  //calculate total catch likelihood components
  if(fscyear>fyear)  //not all catch is sex-specific
  {
    Lc(3)=lognorm(TC_obs(fyear,fscyear-1),TC(fyear,fscyear-1),C_var);
  }
  else
  {
    Lc(3)=0.0;
  }

  //calculate likelihood component for recruitment deviations
  Lrdev=0.5*size(log_rec_devs)*log(rec_var)+0.5*norm2(log_rec_devs)/rec_var;

  //Calculate penalty for differences between male and female F for years 
  //before sex-specific catch
  //Calculate mean F ratio
  F_pen_mu=sum(elem_div(F(2)(fscyear,2006),F(1)(fscyear,2006)))/double(2006-fscyear+1);
  //calculate variance of F_ratio
  F_pen_var=norm2(elem_div(F(2)(fscyear,2006),F(1)(fscyear,2006))-F_pen_mu)/double(2006-fscyear-1);
  //calculate penalty based on normal distribution of F ratio
  F_pen=0.5*norm2(elem_div(F(2)(fyear,fscyear-1),F(1)(fyear,fscyear-1))-F_pen_mu)/(F_pen_var+.0001);  //small constant added so quotient is defined
  
  negLL=sum(Lsa)+sum(Lsr)+sum(Lc)+Lrdev+F_pen;

FUNCTION calculate_reference_points 
  //Calculate equilibrium abundance and catch under each of the Fs
  //Modified for a two-sex model following Shepherd (1982)

  Frat=Frat_init;
  i=0;
  j=0;

  for(i=1;i<=Frat_num;i++)
  {
    //i++;  //counter to keep track of results
    FSPR=FSPR_init;
    for(j=1;j<=FSPR_num;j++)
    {
      //calculate SPR
      SPRf(i,j)=calc_SPR(FSPR,M(1),sex_r);
      SPRm(i,j)=calc_SPR(FSPR*Frat,M(2),1.-sex_r);

      //calculate YPR
      YPRf(i,j)=calc_YPR(FSPR,M(1),sex_r);   //female
      YPRm(i,j)=calc_YPR(FSPR*Frat,M(2),1.-sex_r);  //male

      //calculate equilibrium total N
      N_eq(i,j)=(log(SPRf(i,j))+log_alpha+0.5*rec_var)/beta;

      //calculate equilibrium female N
      Nf_eq(i,j)=N_eq(i,j)*SPRf(i,j)/(SPRf(i,j)+SPRm(i,j));
      Nfs_eq(i,j)=Nf_eq(i,j)*exp(sp_time*(M(1)+FSPR));   

      //calculate equilibrium male N
      Nm_eq(i,j)=N_eq(i,j)-Nf_eq(i,j); 
      Nms_eq(i,j)=Nm_eq(i,j)*exp(sp_time*(M(2)+FSPR*Frat));    
 
      //calculate equilibrium recruitment
      R_eq(i,j)=Nf_eq(i,j)/SPRf(i,j);

      //calculate female equilibrium recruitment
      Rf_eq(i,j)=sex_r*R_eq(i,j);

      //calculate male equilibrium recruitment
      Rm_eq(i,j)=R_eq(i,j)-Rf_eq(i,j);

      //calculate equilibrium catch in numbers
      C_eq(i,j)=R_eq(i,j)*(YPRf(i,j)+YPRm(i,j));
      
      //calculate female equilibrium catch
      Cf_eq(i,j)=R_eq(i,j)*YPRf(i,j);

      //calculate male equilibrium catch
      Cm_eq(i,j)=R_eq(i,j)*YPRm(i,j);

      //calculate female exploitation rate
      uf_eq(i,j)=FSPR/(FSPR+M(1))*(1.-exp(-(FSPR+M(1))));
      uf0_eq(i,j)=Cf_eq(i,j)/(Nfs_eq(i,j)+Rf_eq(i,j));

      //calculate male exploitation rate
      um_eq(i,j)=FSPR*Frat/(FSPR*Frat+M(2))*(1.-exp(-(FSPR*Frat+M(2))));

      //calculate ratio of male to female exploitation rate
      urat_eq(i,j)=Frat;

      //increment the female F for the SPR
      FSPR+=FSPR_inc;
    }

    //increment for ratio of female to male F
    Frat+=Frat_inc;
  }

    //Find F that produces specific %SPR
    for(ispr=1;ispr<=nspr;ispr++)
    {
      FSPR_ref(ispr)=0.5;
      for(iter=1;iter<=niter;iter++)
      {
        SPR_temp=(exp(M(1))-1.)/(exp(M(1)+FSPR_ref(ispr))-1.)-SPR_targ(ispr);
        dSPR=-(exp(2.*M(1)+FSPR_ref(ispr))-exp(M(1)+FSPR_ref(ispr)))/(exp(2.*(M(1)+FSPR_ref(ispr)))-2.*exp(M(1)+FSPR_ref(ispr))+1.);
        FSPR_ref(ispr)-=SPR_temp/dSPR;
      }
    }


FUNCTION dvariable lognorm(dvector obs, dvar_vector est, double var)
  //Function to calculate a lognormal log likelihood
  dvariable L;
  //need to make sure not to include missing data (coded as -99)
  L=0.0;
  double missy=0.0;
  for(y=obs.indexmin();y<=obs.indexmax();y++)
  {
    if(obs(y)!=-99)
    {
      //calculate sum of squares
      L+=square(log(obs(y))-log(est(y)));
    }
    else
    {
      //ass one to the number of missing data points
      missy+=1.0;
    }
  }
  L=0.5*(size(obs)-missy)*log(var)+0.5*L/var;
  return(L);

FUNCTION dvariable normal(dvector obs, dvar_vector est, double cv)
  //Function to calculate a lognormal log likelihood
  dvariable L;
  //need to make sure not to include missing data (coded as -99)
  L=0.0;
  double missy=0.0;
  for(y=obs.indexmin();y<=obs.indexmax();y++)
  {
    if(obs(y)!=-99)
    {
      //calculate sum of squares
      L+=log(cv*est(y))+0.5*square((obs(y)-est(y))/(cv*est(y)));
    }
  }
  return(L);

FUNCTION dvariable calc_SD(dvector obs, dvar_vector est)
  //Function to calculate a standard deviation
  dvariable V;
  //need to make sure not to include missing data (coded as -99)
  V=0.0;
  double missy=0.0;
  for(y=obs.indexmin();y<=obs.indexmax();y++)
  {
    if(obs(y)!=-99)
    {
      //calculate sum of squares
      V+=square(log(obs(y))-log(est(y)));
    }
    else
    {
      //ass one to the number of missing data points
      missy+=1.0;
    }
  }
  V/=(size(obs)-missy);
  return(sqrt(V));

//Overload Baranov catch equation function
FUNCTION dvariable baranov(dvariable F, dvariable M, dvariable N)
  //Baranov catch equation
  return(F/(F+M)*(1.-exp(-(F+M)))*N);

FUNCTION dvariable baranov(prevariable F, dvariable M, dvariable N)
  //Baranov catch equation
  return(F/(F+M)*(1.-exp(-(F+M)))*N);

FUNCTION dvariable baranov(prevariable F, dvariable M, data_number N)
  //Baranov catch equation
  return(F/(F+M)*(1.-exp(-(F+M)))*N);

//Overload SPR and YPR functions so they can take a variable from the DATA_SECTION
FUNCTION dvariable calc_SPR(dvariable F, dvariable M, dvariable srat)
  //SPR based on infinite series solution
  return(srat*(p_F_sp*exp(-sp_time*(M+rf*F))+exp(-((1.+sp_time)*M+rf*F+sp_time*F))/(1.-exp(-(M+F)))));

FUNCTION dvariable calc_SPR(dvariable F, dvariable M, data_number srat)
  //SPR based on infinite series solution
  return(srat*(p_F_sp*exp(-sp_time*(M+rf*F))+exp(-((1.+sp_time)*M+rf*F+sp_time*F))/(1.-exp(-(M+F)))));

FUNCTION dvariable calc_YPR(dvariable F, dvariable M, dvariable srat)
  //YPR based on infinite series solution
  dvariable NPR=srat*exp(-(M+rf*F))/(1.-exp(-(M+F)));
  return(baranov(rf*F,M,srat)+baranov(F,M,NPR));

FUNCTION dvariable calc_YPR(dvariable F, dvariable M, data_number srat)
  //YPR based on infinite series solution
  dvariable NPR=srat*exp(-(M+rf*F))/(1.-exp(-(M+F)));
  return(baranov(rf*F,M,srat)+baranov(F,M,NPR));

FUNCTION dvariable nspawn(dvariable F, dvariable M, dvariable N, dvariable R)
  return (N*exp(-sp_time*(M+F))+p_F_sp*R*exp(-sp_time*(M+rf*F)));

FUNCTION double size(dvector obs)
  return(double(obs.indexmax()-obs.indexmin()+1));

FUNCTION double size(param_init_bounded_dev_vector obs)  //overload size function
  return(double(obs.indexmax()-obs.indexmin()+1));

//******************************************************************************
//Functions to write report files
//******************************************************************************
  
FUNCTION general_report
  ofstream ofs_gen("gen_results.dat");
  {
    ofs_gen << "Name Value" << endl;
    ofs_gen << "negLL " << negLL <<endl;
    ofs_gen << "Lsa_1f " << Lsa(1,1) <<endl;
    ofs_gen << "Lsa_1m " << Lsa(1,2) <<endl;
    ofs_gen << "Lsa_2f " << Lsa(2,1) <<endl;
    ofs_gen << "Lsa_2m " << Lsa(2,2) <<endl;
    ofs_gen << "Lsa_3f " << Lsa(3,1) <<endl;
    ofs_gen << "Lsa_3m " << Lsa(3,2) <<endl;
    ofs_gen << "Lsr_1 " << Lsr(1) << endl;
    ofs_gen << "Lsr_2 " << Lsr(2) << endl;
    ofs_gen << "Lsr_3 " << Lsr(3) << endl;
    ofs_gen << "Lc_f " <<  Lc(1) << endl;
    ofs_gen << "Lc_m " <<  Lc(2) << endl;
    ofs_gen << "Lc_t " <<  Lc(3) << endl;
    ofs_gen << "Lrdev " << Lrdev << endl;
    ofs_gen << "F_pen " << F_pen << endl;
    ofs_gen << "sa_be_1 " << sa_be(1) << endl;
    ofs_gen << "sa_be_2 " << sa_be(2) << endl;
    ofs_gen << "sa_be_3 " << sa_be(3) << endl;
    ofs_gen << "sa_time_1 " << sa_time(1) << endl;
    ofs_gen << "sa_time_2 " << sa_time(2) << endl;
    ofs_gen << "sa_time_3 " << sa_time(3) << endl;
    ofs_gen << "ads_sds_1 " << ad_survey_sds(1) << endl;
    ofs_gen << "ads_sds_2 " << ad_survey_sds(2) << endl;
    ofs_gen << "ads_sds_3 " << ad_survey_sds(3) << endl;
    s=1;
    x=1;
    ofs_gen << "adsf_sds_1 " << calc_SD(ad_survey_obs(s,x)(fy(s),lyear),ad_survey_est(s,x)(fy(s),lyear)) << endl;
    x=2;
    ofs_gen << "adsm_sds_1 " << calc_SD(ad_survey_obs(s,x)(fy(s),lyear),ad_survey_est(s,x)(fy(s),lyear)) << endl;
    s=2;
    x=1;
    ofs_gen << "adsf_sds_2 " << calc_SD(ad_survey_obs(s,x)(fy(s),lyear),ad_survey_est(s,x)(fy(s),lyear)) << endl;
    x=2;
    ofs_gen << "adsm_sds_2 " << calc_SD(ad_survey_obs(s,x)(fy(s),lyear),ad_survey_est(s,x)(fy(s),lyear)) << endl;
    s=3;
    x=1;
    ofs_gen << "adsf_sds_3 " << calc_SD(ad_survey_obs(s,x)(fy(s),lyear),ad_survey_est(s,x)(fy(s),lyear)) << endl;
    x=2;
    ofs_gen << "adsm_sds_3 " << calc_SD(ad_survey_obs(s,x)(fy(s),lyear),ad_survey_est(s,x)(fy(s),lyear)) << endl;
    ofs_gen << "res_sds_1 " << re_survey_sds(1) << endl;
    ofs_gen << "res_sds_2 " << re_survey_sds(2) << endl;
    ofs_gen << "res_sds_3 " << re_survey_sds(3) << endl;
    s=1;
    ofs_gen << "adsf_sds_1 " << calc_SD(re_survey_obs(s)(fy(s),lyear),re_survey_est(s)(fy(s),lyear)) << endl;
    s=2;
    ofs_gen << "adsf_sds_2 " << calc_SD(re_survey_obs(s)(fy(s),lyear),re_survey_est(s)(fy(s),lyear)) << endl;
    s=3;
    ofs_gen << "adsf_sds_3 " << calc_SD(re_survey_obs(s)(fy(s),lyear),re_survey_est(s)(fy(s),lyear)) << endl;
    ofs_gen << "p_rec " << p_rec << endl;
    ofs_gen << "Mf " << M(1) << endl;
    ofs_gen << "Mm " << M(2) << endl;
    ofs_gen << "rec_sd " << rec_sd << endl;
    ofs_gen << "sp_time " << sp_time <<endl;
    ofs_gen << "rf " << rf <<endl;
    ofs_gen << "lyear " << lyear <<endl;
    //catchability
    ofs_gen << "qa_1f " << qa(1,1) << endl;
    ofs_gen << "qa_1m " << qa(1,2) << endl;
    ofs_gen << "qa_2f " << qa(2,1) << endl;
    ofs_gen << "qa_2m " << qa(2,2) << endl;
    ofs_gen << "qa_3f " << qa(3,1) << endl;
    ofs_gen << "qa_3m " << qa(3,2) << endl;
    ofs_gen << "qr_1 " << qr(1) << endl;
    ofs_gen << "qr_2 " << qr(2) << endl;
    ofs_gen << "qr_3 " << qr(3) << endl;
    //SR parameters
    ofs_gen << "alpha " << alpha << endl;
    ofs_gen << "beta " << beta << endl;
    //sex ratio
    ofs_gen << "sex_r " << sex_r << endl;
  }

FUNCTION obs_pred
  ofstream ofs_op("obs_pred_results.dat");
  {
    ofs_op << "year sex a_r s_c snum obs pred" << endl;
    for(y=fyear;y<=lyear-1;y++)
    {
      //total observed and predicted catch
      ofs_op << y << " t a c 0 " << TC_obs(y) << " " << TC(y) <<  endl;
    }

    //sex-specific catch
    for(y=fscyear;y<=lyear-1;y++)
    {
      for(x=1;x<=2;x++)
      {
        ofs_op << y << " " << x << " a c 0 " << C_obs(x,y) << " " << C(x,y) << endl;
      }
    }
    
    //sex-specific survey CPUE
    for(s=1;s<=nsurveys;s++)
    {
      for(y=fy(s);y<=lyear;y++)
      {
        for(x=1;x<=2;x++)
        {
          //adult surveys
          ofs_op << y << " " << x << " a s " << s << " " << ad_survey_obs(s,x,y) << " " << ad_survey_est(s,x,y) << endl;
        }
        //recruit surveys
          ofs_op << y << " 0 r s " << s << " " << re_survey_obs(s,y) << " " << re_survey_est(s,y) << endl;
      }
    }
    //recruitment deviations
    for(y=fyear;y<=lyear-1;y++)
    {
          ofs_op << y+1 << " r r r 0 " << alpha*SP(1,y)*exp(-beta*(SP(1,y)+SP(2,y)))  << " " << R(1,y+1)/sex_r  <<  endl;
    }
  }

FUNCTION HPD_estimates
  ofstream ofs_hpd("HPD_results.dat");
  {
    ofs_hpd << "Year Adult_F Adult_M Rec u_F u_M Frat N_SP_F N_SP_T u_Ft" << endl;
    for(y=fyear;y<=lyear;y++)
    {
      if(y<lyear)
      {
      ofs_hpd << y << " " << N(1,y) << " " << N(2,y) << " " << R(1,y)/sex_r << " " << F(1,y)/(F(1,y)+M(1))*(1.-exp(-(F(1,y)+M(1)))) << " " << F(2,y)/(F(2,y)+M(2))*(1.-exp(-(F(2,y)+M(2)))) << " " << F(2,y)/F(1,y) << " " << SP(1,y)<< " " <<SP(1,y)+SP(2,y) << " " << C(1,y)/(R(1,y)+N(1,y)) << endl;
      }
      else
      {
        ofs_hpd << y << " " << N(1,y) << " " << N(2,y) << " " << R(1,y)/sex_r << " " << "NA" << " " << "NA" << " " << "NA" << " " << "NA"<< " " << "NA" << " NA" << endl;
      }
    }
  }

FUNCTION MSY_estimates
  ofstream ofs_msy("MSY_results.dat");
  {
    //Column headings
    ofs_msy << "n u_fem F_rat SY SYf SYm Nf Nm R SPRf u_fem0" << endl;

    Frat=Frat_init;
    for(i=1;i<=Frat_num;i++)
    {
      FSPR=FSPR_init;
      for(j=1;j<=FSPR_num;j++)
      {
        ofs_msy << j << " " << uf_eq(i,j) << " " << urat_eq(i,j) << " " << C_eq(i,j) << " " << Cf_eq(i,j) << " " << Cm_eq(i,j) << " " << Nfs_eq(i,j) << " " << Nms_eq(i,j) << " " << R_eq(i,j) << " " << SPRf(i,j)/SPRf(1,1) << " " << uf0_eq(i,j) << endl;
      }

      //increment for ratio of female to male F
      Frat+=Frat_inc;
    }
  }

REPORT_SECTION
  //Call reporting functions
  //cout << "Starting report section" << endl;
  general_report();
  obs_pred();
  //cout << "wrote model fit" << endl;
  //write out results to file
  calculate_reference_points();
  //cout << "reference point calculations" << endl;
  MSY_estimates();
  //cout << "output reference points" << endl;
  HPD_estimates();
  //cout << "output HPD estimates" << endl;
  
  report << "Beginning of report section" << endl;
  report << "Likelihood components" << endl;
  report << "Lsa" << endl;
  report << Lsa << endl;
  report << "Lsr" << endl << Lsr <<endl;
  report << "Lc" << endl << Lc << endl;
  report << "Lrdev" << endl << Lrdev << endl;
  report << "F_pen" << endl << F_pen << endl;
  report << endl;
  report << "Variables that define scenarios" << endl;
  report << "VIMS_s_time MDT_s_time WDS_time" << endl;
  report <<  sa_time << endl;
  report << "CVs for indices of abundance" << endl;
  report << "VIMS_a MDT_a WDS_a VIMS_r MDT_r WDS_r" << endl;
  report << ad_survey_sds << " " << re_survey_sds << endl;
  report << "Mf Mm R_sd SP_m prf lyear" << endl;
  report << M << " " << rec_sd << " " << sp_time << " " << rf << " " << lyear  << endl;
  report << "Adult catchability" << endl;
  report << qa << endl;
  report << "recruit catchability" << endl;
  report << qr << endl;

RUNTIME_SECTION
  maximum_function_evaluations 5000, 25000, 20000, 20000, 20000, 20000

//Leave space below this line

