DATA_SECTION
  init_int styr
  init_int endyr
  init_int nages
  init_int nselages
  init_int nsamples
  init_vector catch_bio(styr,endyr)
  init_int nobs_fish
  init_ivector yrs_fish(1,nobs_fish)
  init_int nobs_srv1
  init_ivector yrs_srv1(1,nobs_srv1)
  init_vector obs_srv1(1,nobs_srv1)
  init_matrix obs_p_fish(1,nobs_fish,1,nages)
  init_matrix obs_p_srv1(1,nobs_srv1,1,nages)
  init_vector wt(1,nages)
  init_vector maturity(1,nages)
  int styr_rec
 LOCAL_CALCS
  styr_rec=styr-nages+1;
  cout<<nobs_srv1<<endl;
  cout<<wt<<endl;
 END_CALCS
  

INITIALIZATION_SECTION
  M .20
  mean_log_rec 7.2
  log_avg_fmort -3.
  q1 1.

PARAMETER_SECTION
  init_bounded_number q1(.5,2,8)
  init_bounded_number M(.02,.25,-1)
  init_number mean_log_rec(1)
  init_bounded_dev_vector rec_dev(styr_rec,endyr,-5,5,3)
 
  init_number log_avg_fmort(2)
  init_bounded_dev_vector fmort_dev(styr,endyr,-15,15,2)

  init_vector log_selcoffs_fish(1,nselages,4)
  init_vector log_selcoffs_srv1(1,nselages,4)

  vector log_sel_fish(1,nages)
  vector log_sel_srv1(1,nages)
  vector sel(1,nages)
  vector sel_srv1(1,nages)
  number avgsel_fish;
  number avgsel_srv1;
  vector popn(styr,endyr)
  vector totn_srv1(styr,endyr)
  vector explbiom(styr,endyr)
  vector pred_srv1(styr,endyr)
  matrix pred_p_fish(styr,endyr,1,nages)
  matrix pred_p_srv1(styr,endyr,1,nages)
  vector pred_catch(styr,endyr)
  matrix natage(styr,endyr,1,nages)
  matrix catage(styr,endyr,1,nages)
  //matrix u(styr,endyr,1,nages)
  matrix Z(styr,endyr,1,nages)
  matrix F(styr,endyr,1,nages)
  matrix S(styr,endyr,1,nages)
  number rbar
  number surv
  vector offset(1,2)
  number rec_like
  number catch_like
  vector age_like(1,2)
  vector sel_like(1,2)
  number fpen    
  number surv_like
  sdreport_number endbiom
  sdreport_number depletion
  objective_function_value f

PRELIMINARY_CALCS_SECTION
  //Compute offset for multinomial
      for (int i=1; i <= nobs_fish; i++)
        {
      obs_p_fish(i)=obs_p_fish(i)/sum(obs_p_fish(i));
      //cout<<obs_p_fish(i)<<endl;
          for (int j=1; j<=nages; j++)
          {
            if (obs_p_fish(i,j)>0.0)
            {
              offset(1)-=nsamples*obs_p_fish(i,j)*log(obs_p_fish(i,j));
            }
          }
        }
      for (i=1; i < nobs_srv1; i++)
        {
      obs_p_srv1(i)=obs_p_srv1(i)/sum(obs_p_srv1(i));
      //cout<<obs_p_srv1(i)<<endl;
          for (int j=1; j<=nages; j++)
          {
            if (obs_p_srv1(i,j)>0.0)
            {
              offset(2)-=nsamples*obs_p_srv1(i,j)*log(obs_p_srv1(i,j));
            }
      }
    }
     cout<<endl<<offset<<endl;


PROCEDURE_SECTION
  // example of using FUNCTION to structure the procedure section
  get_selectivity();
  get_mortality();
  surv=mfexp(-1.0* M);
  get_numbers_at_age();
  get_catch_at_age();
  evaluate_the_objective_function();
  //cout<<"rec_like   age_like  fpen       catch_like       surv_like"<<endl;
  //cout<<rec_like<<"  "<< age_like<< " "<<fpen<<" "<<catch_like<<" "<<surv_like<<endl;

FUNCTION get_selectivity
  // calculate the selectivity from alpha beta
  for (int j=1;j<=nselages;j++)
  {
    log_sel_fish(j)=log_selcoffs_fish(j);
    log_sel_srv1(j)=log_selcoffs_srv1(j);
  }
  for (j=nselages+1;j<=nages;j++)
  {
    log_sel_fish(j)=log_sel_fish(j-1);
    log_sel_srv1(j)=log_sel_srv1(j-1);
  }
  
  avgsel_fish=log(mean(mfexp(log_selcoffs_fish)));
  avgsel_srv1=log(mean(mfexp(log_selcoffs_srv1)));
  log_sel_fish-=log(mean(exp(log_sel_fish)));
  log_sel_srv1-=log(mean(exp(log_sel_srv1)));
  sel=mfexp(log_sel_fish);
  sel_srv1=mfexp(log_sel_srv1);


FUNCTION get_mortality
  for (int i=styr;i<=endyr;i++)
    {
    for (int j=1;j<=nages;j++)
       {
        F(i,j)=sel(j)*mfexp(log_avg_fmort + fmort_dev(i));
       }
    }
    Z=F+M;
    S=mfexp(-1.0*Z);

FUNCTION get_numbers_at_age
  int itmp;
  for (int j=1;j<nages;j++)
    {
      itmp=styr+1-j;
      natage(styr,j)=mfexp(mean_log_rec-M*double(j-1)+rec_dev(itmp));
    }
    itmp=styr+1-nages;
    natage(styr,nages)=mfexp(mean_log_rec+rec_dev(itmp)-M*(nages-1))/(1. - surv);
   
 // Now do for next several years----------------------------------
  for (int i=styr+1;i<=endyr;i++)
    {
      natage(i,1)=mfexp(mean_log_rec+rec_dev(i));
    }

  for (i=styr;i< endyr;i++)
    {
       popn(i)= natage(i)*sel;
       natage(i+1)(2,nages)=elem_prod(natage(i)(1,nages-1),S(i)(1,nages-1))++;
       natage(i+1,nages)+=natage(i,nages)*S(i,nages);
    }
    popn(endyr)=natage(endyr)* sel;
  for (i=styr;i<=endyr;i++)
    {
     totn_srv1(i)=(natage(i)*sel_srv1);
     pred_srv1(i)=q1*(natage(i)*elem_prod(sel_srv1,wt));
       explbiom(i)=natage(i)*elem_prod(sel,wt);
       pred_p_srv1(i)=elem_prod(sel_srv1,natage(i))/totn_srv1(i);
    }
    depletion=explbiom(endyr)/explbiom(1);
    endbiom=explbiom(endyr);

FUNCTION get_catch_at_age
  for (int i=styr; i<=endyr; i++)
  {
     pred_catch(i)=0.;
   //--Baranov's equation here-----------------------------------
    for (int j = 1 ; j<= nages; j++)
      {
       catage(i,j) = natage(i,j)*F(i,j)*(1.-S(i,j))/Z(i,j);
       pred_catch(i)+=catage(i,j)*wt(j);
      }
      pred_p_fish(i)=catage(i)/sum(catage(i));
    }

FUNCTION evaluate_the_objective_function
  f=.0;
  if (active(rec_dev))
  {
    age_like=0.;
    int ii;
    rec_like=norm2(rec_dev);
    f+=2.*rec_like;
    for (int i=1; i <= nobs_fish; i++)
    {
      ii=yrs_fish(i);
      //cout<<pred_p_fish(ii)<<endl<<endl;
      for (int j=1; j<=nages; j++)
      {
        age_like(1)-=nsamples*obs_p_fish(i,j)*log(pred_p_fish(ii,j)+1.e-13);
       //cout << age_like << " " << i << " "<<j<< endl;
      }
    }
    age_like(1)-=offset(1);
    for (i=1; i < nobs_srv1; i++)
      {
        ii=yrs_srv1(i);
        for (int j=1; j<=nages; j++)
        {
          age_like(2)-=nsamples*obs_p_srv1(i,j)*log(pred_p_srv1(ii,j)+1.e-13);
         //cout << age_like << " " << i << " "<<j<< endl;
        }
      }
      age_like(2)-=offset(2);
      f+=sum(age_like);
  }

  // Fit to indices (lognormal)                       
    surv_like=norm2(log(obs_srv1+.01)-log(pred_srv1(yrs_srv1)+.01));

    catch_like=100.*norm2(log(catch_bio+.000001)-log(pred_catch+.000001));
 
    sel_like(1)=10.*norm2(first_difference(first_difference(log_sel_fish)));
    sel_like(2)=10.*norm2(first_difference(first_difference(log_sel_srv1)));
  
    f+=10.*sum(sel_like);
    f+=10.*square(avgsel_fish);
    f+=10.*square(avgsel_srv1);

    f+=10.*surv_like;

    f+=catch_like;

  // Phases less than 3, penalize High F's
    if (current_phase()<3)
    {
     fpen=10.*norm2(mfexp(fmort_dev+log_avg_fmort)-.2);
    }
     else
    {
     fpen=.01*norm2(mfexp(fmort_dev+log_avg_fmort)-.2);
    }

    if (active(fmort_dev))
	{
     fpen+=norm2(fmort_dev);
	}
    f+=fpen;



REPORT_SECTION
  report << "Estimated numbers of fish " << endl;
  report << natage << endl;
  report << "Estimated F mortality " << endl;
  report << F << endl;
  report << "Observed Survey 1 " << endl;
  report << obs_srv1 << endl;
  report << "Predicted Survey 1 " << endl;
  report << pred_srv1 << endl;
  report << "Observed Prop " << endl;
  report << obs_p_fish << endl;
  report << "Predicted prop  " << endl;
  report << pred_p_fish << endl;
  report << "Predicted prop  " << endl;
  report << pred_p_srv1<< endl;
  report << "Observed catch biomass " << endl;
  report << catch_bio << endl;
  report << "predicted catch biomass " << endl;
  report << pred_catch << endl;
  report << "Estimated annual fishing mortality " << endl;
  report << mfexp(log_avg_fmort+fmort_dev) << endl;
  report << "Selectivity" << endl;
  report << sel << endl;
  report << sel_srv1 << endl;
  report<<"rec_like   age_like  fpen sel_like  catch_like  surv_like"<<endl;
  report <<rec_like<<"  "<< age_like<< " "<<fpen<<" "<<sel_like<<" "<<catch_like<<" "<<surv_like<<endl;