//ADModel Builder template for Goose Island Gully Pacific ocean perch
// Laura Richards, July 1996 (richardsl@pbs.dfo.ca)


// This template was used to produce results in
//    Richards, Schnute and Olsen, 1997, CJFAS 54: July issue
// based on the model described by
//    Schnute and Richards, 1995, CJFAS 52: 2963-2077

// Data notes: Age data are missing for 1986 and 1988; Survey 1
//        (research vessel) ended in 1984; Survey 2 (commercial vessel)
//        covers 1984-95, with two surveys in 1995


DATA_SECTION
  init_int NYear                       // number of years of data
  init_int NAge                        // number of age classes
  init_int NAge1                       // max. age for surface reading
  init_number PMin                     // min p_at for grouping classes
  init_matrix data_agewgt(1,759,1,4)   // yr, age, prop, wgt
  init_matrix data_cpue(1,NYear,1,5)   // yr, catch, suv1, surv2, cat_frac

  int NITerms                          // counter re grouping classes
  int NPTerms                          // counter re grouping classes

  vector cf(1,NYear)                   // proportion catch caught at survey
  vector D_t(1,NYear)                  // observed biomass in year t

  matrix w_at(1,NAge,1,NYear)          // observed weight at age in year t

PARAMETER_SECTION
  init_bounded_number F(0.0,0.3,-1)    // historic fishing mortality
  init_bounded_number M(0.02,0.08)     // natural mortality
  init_bounded_number R(0.1,10.0)      // historic mean recruitment
  init_bounded_number gamma(-1.0,1.0)  // recruitment correlation
  init_bounded_number alpha(0.02,50.0) // selectivity exponent
  init_bounded_number beta(0.0,1.0)    // selectivity at age class 1
  init_bounded_number q(0.1,10.0)      // catchability for research
  init_bounded_number q2(0.1,10.0)     // catchability for commerc.
  init_number rho(-1)                  // ratio of sig1^2 and kap^2
// Recruitments at time t: estimate in phase 2,
  init_bounded_vector Rt(2-NAge1,NYear,1.0E-03,50.0,2)

  number kap_2                         // rec+index var. Eq. L.11
  number sig2_2                        // process error (survival, lognormal)
  number tau2_2                        // measurement error (logit-normal)

  vector B_t(1,NYear)                  // Biomass at time t
  vector beta_a(1,NAge)                // Selectivity at age a
  vector C_t(1,NYear)                  // Catch in numbers at time t
  vector P_t(1,NYear)                  // Expl. pop. at time t
  vector Pred_It(1,NYear)              // Abund. index at time t
  vector Pred_I2(1,NYear)              // Abund. index 2 at time t

  matrix Nat(1,NAge,1,NYear)           // Numbers at age a, time t
  matrix u_at(1,NAge,1,NYear)          // Expl. proportion at age a, time t

// State moments and dynamics
  number S                             // exp(-M), survival
  vector betaNat(1,NAge)               // elem_prod of beta_a*N_at

// Residuals
  vector xi_t(2-NAge1,NYear)           // xi_t in Eq. L.4
  number xiSS                          // sum of squares of xi_t in Eq. L.4
  vector eta_0t(1,NYear)               // eta_0t in Eq. L.5
  vector eta_t(1,NAge)                 // eta_at, for age a
  number eta0tSS                       // sum of squares of eta_0t
  number sumu                          // Collects terms in Eq. L.6
  number sumeta                        // Collects terms in Eq. L.6

// Likelihood terms

  vector like_terms(1,5)               // Eq L.18 terms
  number penalty                       // a penalty value
  objective_function_value lf          // the objective function value

// Standard Deviation Report

  sdreport_vector B_Final(1,NYear)

PRELIMINARY_CALCS_SECTION
  // extract proportion of catch taken prior to index.
  cf=column(data_cpue,5);
  D_t=column(data_cpue,2);             // Catch biomass at time t
  int cnt=1;
  for (int t=1;t<=NYear;t++)
    for (int a=1;a<=NAge;a++,cnt++)
      w_at[a][t]=data_agewgt[cnt][4];  // Observed weight at age a, time t

PROCEDURE_SECTION
  Calc_Selectivity();
  Calc_Nat_and_Moments();
  Calc_Pred_It();
  Calc_tau2_2();
  Calc_kap_2();
  Calc_Objective_Function();


FUNCTION Calc_Selectivity
  for (int a=1;a<NAge;a++)
    beta_a[a]=1.0-(1.0-beta)*pow(double(NAge-a)/double(NAge-1.0),alpha);
  beta_a[NAge]=1.0;


FUNCTION Calc_Nat_and_Moments
  int mage;
  S = exp( -M );

  for (int t=1;t<=NYear;t++)
  {
    if ( t==1 )
    {
      Nat[1][1] = Rt[1];
      for (int a=2;a<=NAge1;a++)
      {
        Nat[a][1] = Rt[2-a] * exp( -(F+M) * (a-1.0) );
      }
      Nat[NAge1][1] = Nat[NAge1][1] / (1.0 - exp( -(F+M) ));
    }
    else
    {
      mage=NAge1+t-1;
      if ( mage > NAge )
        mage = NAge;  // add age class each year to NAge
      Nat[1][t]=Rt[t];
      for (int a=2;a<=mage;a++)
        Nat[a][t]=S*( Nat[a-1][t-1]-u_at[a-1][t-1]*C_t[t-1] );
      Nat[mage][t]=S*( Nat[mage-1][t-1]+Nat[mage][t-1] -
                   (u_at[mage-1][t-1]+u_at[mage][t-1])*C_t[t-1]);
    }

    // Calculate the state moments.
    betaNat = elem_prod( beta_a, extract_column( Nat,t ) );
    P_t[t] = sum( betaNat );
    u_at.colfill(t, (betaNat / P_t[t]) );
    B_t[t] = sum( elem_prod( betaNat,extract_column(w_at,t) ) );
    C_t[t] = D_t[t] / sum( elem_prod(extract_column(u_at,t),
                                     extract_column(w_at,t)) );

  }     // for t=1,NYear
  if ( sd_phase() )
  {
    B_Final = log(B_t);
  }


FUNCTION Calc_Pred_It
  for (int t=1;t<=NYear;t++)
  {
    Pred_It[t] = q  * ( B_t[t] - cf[t] * D_t[t] );
    Pred_I2[t] = q2 * ( B_t[t] - cf[t] * D_t[t] );
  }


FUNCTION Calc_tau2_2
  int a;
  NPTerms = 0;
  tau2_2 = 0.0;
  for (int t=1;t<=NYear;t++)
  {
    if ( (t!=24) && (t!=26) )           // Missing data 1986 and 1988
    {
       int nages=1;
       double sumprop = 0.0, remprop=1.0;
       sumu = 0.0; sumeta = 0.0;

       for (a=1;a<NAge;a++)
       {
         remprop -= data_agewgt[ (t-1)*NAge+a ][3];
         sumprop += data_agewgt[ (t-1)*NAge+a ][3];
         sumu += u_at[a][t];
         if ((sumprop>=PMin) && (remprop>=PMin) && ((t>14) || (a<NAge1)) )
         {
           eta_t[a] = log( sumprop+1.0E-10 ) - log( sumu+1.0E-10 );
           sumeta += eta_t[a];
           sumu = 0.0;
           sumprop = 0.0;
           nages++;
           NPTerms++;
         }
         else
           eta_t[a] = 0.0;
       }    // for a=1,NAge

       // Ensure last age class.
       sumprop += data_agewgt[ (t-1)*NAge+NAge ][3];
       sumu += u_at[NAge][t];

       eta_t[NAge] = log( sumprop ) - log( sumu );
       sumeta += eta_t[NAge];

       for (a=1;a<=NAge;a++)
       {
         if ( value( eta_t[a] ) != 0.0 )
           eta_t[a] = eta_t[a] - ( sumeta / (double)nages );
           tau2_2 += eta_t[a] * eta_t[a];
       }     // for a=1,NAge

    }     // if t<>24 and t<>26
  }     // for t=1,NYear
  tau2_2 =  tau2_2 / (double)NPTerms;


FUNCTION Calc_kap_2
  int t;
  xiSS = 0.0; eta0tSS = 0.0;
  kap_2 = 0.0;

  xi_t[2-NAge1] = log( Rt[2-NAge1] ) - log( R );
  for (t=3-NAge1;t<=NYear;t++)
  {
    xi_t[t] = log(Rt[t]) - ((1.0-gamma)*log(R)) - (gamma*log(Rt[t-1]));
    xiSS += xi_t[t] * xi_t[t];
  }
  NITerms = 0;
  for (t=1;t<=NYear;t++)
  {
  // if the survey biomass estimate positive...
    if ( data_cpue[t][3] > 0.0 )
    {
      // fix for extra survey in 1995, data stuffed into survey 1 column.
      if ( t==NYear )
        eta_0t[t] = log( data_cpue[t][3] ) - log( Pred_I2[t] );
      else
        eta_0t[t] = log( data_cpue[t][3] ) - log( Pred_It[t] );
      eta0tSS += eta_0t[t] * eta_0t[t];
      NITerms++;
    }

    if ( data_cpue[t][4] > 0.0 )
    {
      eta_0t[t] = log( data_cpue[t][4] ) - log( Pred_I2[t] );
      eta0tSS += eta_0t[t] * eta_0t[t];
      NITerms++;
    }
  }     // for t=1,NYear

  kap_2 = (((1.0-gamma*gamma)*square(xi_t[2-NAge1]) + xiSS)/rho) +
          (eta0tSS / (1.0-rho) );
  kap_2 = kap_2 / (double)(NAge1+NYear-2+NITerms);


FUNCTION Calc_Objective_Function
  // NOTE: twice negative log-likelihood: adm requires neg log-lk
  //                so mulitply by 0.5
  like_terms[1] = (double)(NAge1+NYear-2) * log( rho );
  like_terms[2] = (double)(NITerms) * log(1.0-rho);
  like_terms[3] = -log(1.0-gamma*gamma);
  like_terms[4] = (double)(NAge1+NYear-2+NITerms) * log(kap_2);
  like_terms[5] = (double)(NPTerms) * log(tau2_2);
  lf = 0.5*sum( like_terms ) + penalty;


RUNTIME_SECTION
  convergence_criteria 1.0E-06


REPORT_SECTION
  report << "NYear   : " << NYear      << endl;
  report << "NAge1   : " << NAge1      << endl;
  report << "NAge    : " << NAge       << endl;
  report << "PMin    : " << PMin       << endl;
  report << endl;
  report << "ObjFun  : " << lf         << endl;
  report << endl;
  report << "M       : " << M          << endl;
  report << "R       : " << R          << endl;
  report << "gamma   : " << gamma      << endl;
  report << "alpha   : " << alpha      << endl;
  report << "beta    : " << beta       << endl;
  report << "q       : " << q          << endl;
  report << "q2      : " << q2         << endl;
  report << "rho     : " << rho        << endl;
  report << "kap_2   : " << kap_2      << endl;
  report << "tau2_2  : " << tau2_2     << endl;
  report << "B_t[T]  : " << B_t[NYear] << endl;