DATA_SECTION
  init_int startyr
  init_int nyrs  
  init_int nages
  init_int nsamples
  init_vector catch_bio(1,nyrs)
  init_int nobs_srv1
  init_ivector yrs_srv1(1,nobs_srv1)
  init_vector obs_srv1(1,nobs_srv1)
  init_int nobs_srv2
  init_ivector yrs_srv2(1,nobs_srv2)
  init_vector obs_srv2(1,nobs_srv2)
  init_vector frac(1,nyrs)
  init_matrix wt(1,nyrs,1,nages)
  init_matrix obs_p(1,nyrs,1,nages)

INITIALIZATION_SECTION
  M .05
  rho .8
  mean_log_rec 3.2

PARAMETER_SECTION
  init_bounded_number alpha(1.0001,20,4)
  init_bounded_number beta1(.001,.5,4) 
  init_bounded_number q1(.01,1,1)
  init_bounded_number q2(.01,1,1)
  init_bounded_number M(.02,.15,-1)
  init_bounded_number rho(.1,1,-1)
  init_bounded_number mean_log_rec(0.1,10.,1)
  init_bounded_dev_vector rec_dev(-21,nyrs,-5,5,3)
  vector sel(1,nages)
  vector fmort(1,nyrs)
  vector pred_srv1(1,nyrs)
  vector pred_srv2(1,nyrs)
  vector popn(1,nyrs)
  vector catn(1,nyrs)
  vector explbiom(1,nyrs)
  matrix natage(1,nyrs,1,nages)
  matrix pred_p(1,nyrs,1,nages)
  matrix u(1,nyrs,1,nages)
  number rbar
  number tauhatsq
  number kappahatsq
  number surv
  number offset
  number rec_like
  number age_like
  number fpen    
  number stmp    
  number pentmp
  number negpen    
  number surv_like
  objective_function_value f

PRELIMINARY_CALCS_SECTION
  // this is just to ``invent'' some relative average
  // weight at age numbers
  //cout<<obs_p<<endl;
      for (int i=1; i <= nyrs; i++)
        {
          for (int j=1; j<=nages; j++)
          {
            if (obs_p(i,j)>0)
            {
              offset-=nsamples*obs_p(i,j)*log(obs_p(i,j));
            }
          }
        }


PROCEDURE_SECTION
  // example of using FUNCTION to structure the procedure section
  //cout<<alpha<<" "<<beta1<<"  "<<nages<<endl;
  get_selectivity();
  //cout<<"Sel1"<<endl;
  surv=mfexp(-1.0* M);
  get_numbers_at_age();
  //cout<<"numbers"<<endl;
  evaluate_the_objective_function();
  //cout<<"objfunc"<<endl;
  //cout<<rec_like<<"  "<< age_like<< " LikeFuns "<<surv_like<<endl;

FUNCTION get_selectivity
  // calculate the selectivity from alpha beta
  for (int j=1;j<nages;j++)
  {
    sel(j)=1. - (1.-beta1) * pow((double(nages-j)/ double(nages-1.)),alpha);
     //cout<<j<<" "<<sel(j)<<endl;
  }
  sel(nages)=sel(nages-1);


FUNCTION get_numbers_at_age
  int itmp;
  for (int j=1;j<nages;j++)
    {
      itmp=2-j;
      //cout<<itmp<<" "<<j<<endl;
      //Eq. A.9
      natage(1,j)=mfexp(mean_log_rec-M*double(j-1)+rec_dev(itmp));
    }
 // Last age group in first year
    itmp=2-nages;
  //Eq. A.10
    natage(1,nages)=mfexp(mean_log_rec+rec_dev(itmp)-M*(nages-1))/(1. - surv);
    //cout<<natage(1)<<endl;
    //cout<<natage(1,nages)<<" Yr1 N at maxage" <<endl;
   
 // Now do for next several years----------------------------------
  for (int i=2;i<=nyrs;i++)
    {
     // Eq. A.11
      natage(i,1)=mfexp(mean_log_rec+rec_dev(i));
    }
    popn(1)=sel*natage(1);
    u(1)=elem_prod(sel,natage(1))/popn(1);
    catn(1)=catch_bio(1)/(u(1)*wt(1));
    //cout<<1<<" "<<popn(1)<<endl;

   negpen=0.;
  for (i=2;i<=nyrs;i++)
    {
       for (int j=2;j<nages;j++)
       {
    // Eq. A.12
         natage(i,j)=surv*(natage(i-1,j-1)-u(i-1,j-1)*catn(i-1));
         //stmp = natage(i,j)/natage(i-1,j-1);
   //stmp = posfun(stmp,.1,negpen);
         //natage(i,j) = stmp * natage(i-1,j-1);
       }
    // Eq. A.13
       natage(i,nages)=surv*(natage(i-1,nages-1)+natage(i-1,nages)-(u(i-1,j-1)+u(i-1,nages))*catn(i-1));


    // Eq. A.3 POPN
       popn(i)=sel*natage(i);
       u(i)=elem_prod(sel,natage(i))/popn(i);
       catn(i)=catch_bio(i)/(u(i)*wt(i));
       stmp = 1.-catn(i)/popn(i);
    stmp = posfun(stmp,.2,negpen);
    //cout<<negpen<<endl;
       catn(i)=(1.-stmp)*popn(i);
       /* if(catn(i) > popn(i))          
        { 
          fpen+=10.*square(catn(i)-popn(i));
          catn(i)=popn(i)*.8;
          cout<<i<<" Ooops "<<popn(i)<<endl;
        }
       else
        {
          fpen=0.;
        }
  */
    }
       //cout<<catn<<" catn"<<endl;

  for (i=1;i<=nyrs;i++)
    {
       explbiom(i)=0.;
       for (int j=1;j<=nages;j++)
       {
       // Eq. A.5
         explbiom(i)+=sel(j)*wt(i,j)*natage(i,j);
         pred_p(i,j)=sel(j)*natage(i,j)/popn(i);
       }
     // Eq. A.8
     fmort(i)=log(popn(i)/(popn(i)-catn(i)));
     // cout<<natage(i,1)<<" ";
    }
       pred_srv1=q1*(explbiom-elem_prod(frac,catch_bio));
       pred_srv2=q2*(explbiom-elem_prod(frac,catch_bio));
       //pred_srv2=q2*(explbiom(i)-frac(i)*catch_bio(i));
       //cout<<"ExpB "<<explbiom<<endl;
       //cout<<frac<<endl;
       //cout<<pred_srv1<<endl;
       //cout<<pred_srv2<<endl;
     //cout<<endl;

FUNCTION evaluate_the_objective_function
  f=.0;
  age_like=0.;
  for (int i=1; i <= nyrs; i++)
    {
      for (int j=1; j<=nages; j++)
      {
        age_like-=nsamples*obs_p(i,j)*log(pred_p(i,j));
      }
    }
    age_like-=offset;
  // Fit to indices (lognormal)                       
  surv_like=norm2(log(obs_srv1+.01)-log(pred_srv1(yrs_srv1)+.01));
  surv_like+=norm2(log(obs_srv2+.01)-log(pred_srv2(yrs_srv2)+.01));
  rec_like=norm2(rec_dev);
  //cout<<surv_like<<" SurvHere2"<<endl;
  //cout<<yrs_srv1<<endl;
  //cout<<pred_srv1(yrs_srv1)<<endl;
  //cout<<obs_srv1<<endl;
  f+=age_like;
  f+=100.*surv_like;
  f+=10.*rec_like;
  // Phases less than 3, penalize High F's
  if (current_phase() < 3)
    {
   fpen+=10.*norm2(fmort-.1);
 }
  f+=fpen;
  f+=10.*negpen;



REPORT_SECTION
  report << "Estimated numbers of fish " << endl;
  report << natage << endl;
  report << "Observed Survey 1 " << endl;
  report << obs_srv1 << endl;
  report << "Predicted Survey 1 " << endl;
  report << pred_srv1 << endl;
  report << "Observed Survey 2 " << endl;
  report << obs_srv2 << endl;
  report << "Predicted Survey 2 " << endl;
  report << pred_srv2 << endl;
  report << "Observed Prop " << endl;
  report << obs_p << endl;
  report << "Predicted prop  " << endl;
  report << pred_p << endl;
  report << "Observed catch biomass " << endl;
  report << catch_bio << endl;
  report << "Estimated annual fishing mortality " << endl;
  report << (fmort) << endl;