DATA_SECTION
  init_number lambda
  init_int styr  
  init_int endyr  
  init_int nages
  init_vector obs_rpn(styr,endyr)
  init_vector obs_totcat(styr,endyr)
  init_int nyrs_age
  init_ivector yrs_age(1,nyrs_age)
  init_int obs_nages
  init_matrix obs_p_age(1,nyrs_age,1,obs_nages)
  init_int nlenbins_m
  init_int nlenbins_f
  init_matrix trans_m(1,obs_nages,1,nlenbins_m)
  init_matrix trans_f(1,obs_nages,1,nlenbins_f)
  init_matrix obs_p_len_m(styr,endyr,1,nlenbins_m)
  init_matrix obs_p_len_f(styr,endyr,1,nlenbins_f)
  init_vector wt_age(1,nages)
  init_vector maturity(1,nages)

INITIALIZATION_SECTION
  log_gamma -20
  log_beta 0.
  log_a50 1.5
  m .1
  log_avg_rec 8
  log_avg_init_age 4
  log_q 1.8

PARAMETER_SECTION
  init_number log_q(1)
  init_number log_gamma(-3)
  init_number log_a50(3)
  init_number log_beta(3)
  init_number m(-1)

  init_number log_avg_rec(1)
  init_bounded_dev_vector log_rec_dev(styr,endyr,-25,25,2)
  init_number log_avg_init_age(1)
  init_bounded_dev_vector init_age_dev(2,nages,-15,15,2)
  init_vector catch_dev(1986,1990,4)
  vector new_totcat(styr,endyr)
  matrix pred_p_len_m(styr,endyr,1,nlenbins_m)
  matrix pred_p_len_f(styr,endyr,1,nlenbins_f)
  matrix N(styr,endyr,1,nages)
  matrix pred_N(styr,endyr,1,obs_nages)
  vector u(styr,endyr)
  vector pred_rpn(styr,endyr)
  matrix pred_p_age(styr,endyr,1,obs_nages)
  vector sel(1,nages)
  vector N_f(styr,endyr)
  number avg_hr
  number surv
  number a50
  number gamma
  number beta
  number like_rec
  number like_init
  number like_age
  number like_len
  number like_hr
  number like_rpn
  sdreport_number biom96

  objective_function_value f

PRELIMINARY_CALCS_SECTION
  // this is just to ``invent'' some relative average
  // weight at age numbers
  new_totcat=obs_totcat;
PROCEDURE_SECTION
  // example of using FUNCTION to structure the procedure section
  if (last_phase())
    {
     for (int i=1986;i<=1990;i++)
       {
       new_totcat(i)=obs_totcat(i)*exp(catch_dev(i));
       }
   }
  get_mortality_and_survivial_rates();
  get_numbers_at_age();
  get_catch_at_age();
  evaluate_the_objective_function();

FUNCTION get_mortality_and_survivial_rates

  surv=exp(-m);

  // calculate the selectivity from the Thompson
  beta=exp(log_beta);
  a50=exp(log_a50);
  gamma=exp(log_gamma);

  for (int j=1;j<=nages;j++)
  {
    sel(j)=(1/(1-gamma))*pow((1-gamma)/gamma,gamma)*(
           exp(beta*gamma*(a50-double(j)))/(1+
           exp(beta*(a50-double(j)))));
  }
  // the selectivity is the same for the last two age classes

FUNCTION get_numbers_at_age
 // Get initial Age composition into N======================
  for (int j=2;j<=nages;j++)
  {
    N(styr,j)=exp(log_avg_init_age + init_age_dev(j));
  }
 // Get recruitment in subsequent years======================
  for (int i=styr;i<=endyr;i++)
  {
    N(i,1)=exp(log_avg_rec + log_rec_dev(i));
  }
 // Get fishable stock, compute rest of N matrix======================
  for (i=styr;i<endyr;i++)
  {
    N_f(i)=N(i)*sel;
    u(i)=new_totcat(i)/N_f(i);
    for (j=1; j<nages-1 ;j++)
    {
      N(i+1,j+1)=N(i,j)*(1-u(i)*sel(j))*surv;
    }
    N(i+1,nages)=N(i,nages-1)*(1-u(i)*sel(nages-1))*surv +
                 N(i,nages)*(1-u(i)*sel(nages))*surv;
  }
  N_f(endyr)=N(endyr)*sel;
  u(endyr)=new_totcat(endyr)/N_f(endyr);
  //cout<<N_f<<endl;
 // Now compute predicted RPN==================================
  pred_rpn=exp(log_q)*N_f;
 // Now compute predicted 96 biomass ==========================
  biom96=elem_prod(N(1996),sel)*wt_age;

FUNCTION get_catch_at_age

  for (int i=styr;i<=endyr;i++)
  {
    pred_N(i)(1,7) = elem_prod(N(i)(1,7),sel(1,7)); 
    pred_N(i,8)    = N(i)(8,9)*sel(8,9); 
    pred_N(i,9)    = N(i)(10,14)*sel(10,14); 
    pred_N(i,10)   = N(i,15)*sel(15);

    pred_p_age(i)  = pred_N(i)/sum(pred_N(i));
  }

  pred_p_len_f = pred_p_age * trans_f; 
  pred_p_len_m = pred_p_age * trans_m; 

FUNCTION evaluate_the_objective_function
  // penalty functions to ``regularize '' the solution

  like_rec  =.01*norm2(log_rec_dev);
  like_init =.01*norm2(init_age_dev);
  
  avg_hr=sum(u)/size_count(u);
  if (current_phase() > 1 )
    {
     like_hr=0.;  // 0.1*square(log(avg_hr/.1));
    }
    else
    {
     like_hr=10.*square(log(avg_hr/.1));
    }
  like_age=0.;
  for(int i=1;i<=nyrs_age;i++)
  {
    like_age -= 50.*sum(elem_prod(obs_p_age(i),
      log(elem_div(1.e-3+pred_p_age(yrs_age(i)),1.e-3+obs_p_age(i)))));
  }
  like_len= -100.*sum(elem_prod(obs_p_len_m,
             log(elem_div(1.e-3+pred_p_len_m,
             1.e-3+obs_p_len_m))));

  like_len -= 100.*sum(elem_prod(obs_p_len_f,
             log(elem_div(1.e-3+pred_p_len_f,
             1.e-3+obs_p_len_f))));

  like_rpn=lambda*(norm2(log(obs_rpn)-log(pred_rpn)));
  if (last_phase())
    {
      f+=1.*norm2(catch_dev);
     }

  f+=like_age+like_len+like_hr+like_rpn + like_rec + like_init;
  

RUNTIME_SECTION
  convergence_criteria 1.e-4 1.e-7  
  maximum_function_evaluations 200,1000

TOP_OF_MAIN_SECTION
  arrmblsize=160000;

REPORT_SECTION
  report << "Estimated numbers of fish " << endl;
  report << N << endl;
  report << "Estimated harvest rate " << endl;
  report << u<<endl;
  report << "Observed RPN" << endl;
  report << obs_rpn<<endl;
  report << "Predicted RPN" << endl;
  report << pred_rpn<<endl;
  report << "Selectivity " << endl;
  report << sel<<endl;
  report << "Pred P Age  " << endl;
  report << pred_p_age <<endl;
  report << "Obsr P Age  " << endl;
  report << obs_p_age <<endl;
  report << "Pred Len Males " << endl;
  report << pred_p_len_m <<endl;
  report << "Obsr len M  " << endl;
  report << obs_p_len_m <<endl;
  report << "Pred Len FeMales " << endl;
  report << pred_p_len_f <<endl;
  report << "Obsr len F  " << endl;
  report << obs_p_len_f <<endl;
  report << "obs, pred catch" << endl;
  report << obs_totcat <<endl;
  report << new_totcat <<endl;
  report << "Wt Maturity " << endl;
  report << wt_age <<endl;
  report << maturity <<endl;
  report << "Like Age, Len, rpn, hr, rec, init "<<endl;
  report << like_age<<" "<<like_len<<" "<<like_rpn<<" "<<like_hr <<" "<<like_rec<<" "<<like_init<<endl;