DATA_SECTION
  init_ivector control_flags(1,6)
  init_number phase_M
  init_number phase_explogsel
  init_number phase_selcoff
  init_int lambda1  
  init_int lambda2    
  init_int lambda3  
  init_int lambda4  
  init_int nyrs    
  init_int nages
  int dim_sel;
  int group_num
  int last_age_sel_group
  int nng
  int c_size
  int s_size
  int snyrs
 LOCAL_CALCS
  snyrs=nyrs-control_flags(6);
  last_age_sel_group=control_flags(3);
  c_size=control_flags(4);
  s_size=control_flags(5);
  nng= nages-last_age_sel_group+1;
  group_num=3;
  int ii=0;
  for (int i=1;i<snyrs;i++)
  {
    if (!(i%group_num))
    {
      ii++;
    }
  }
  dim_sel=ii;
 END_CALCS
  init_matrix agerr(1,nages,1,nages)
  init_matrix sobs_catch_at_age(1,nyrs,1,nages)
  matrix obs_catch_at_age(1,snyrs,1,nages)
  matrix OP(1,snyrs,1,nages)
  vector OTC(1,snyrs)
 LOCAL_CALCS
  for (i=1;i<=snyrs;i++)
  {
    obs_catch_at_age(i)=sobs_catch_at_age(i);
    OTC(i)=sum(sobs_catch_at_age(i));
    OP(i)=sobs_catch_at_age(i)/OTC(i);
  }
 END_CALCS
  init_vector seffort(1,nyrs)
  vector effort(1,snyrs)
  init_matrix sobs_catchsd_at_age(1,nyrs,1,nages)
  matrix obs_catchsd_at_age(1,snyrs,1,nages)
  init_vector relwt(2,nages)
  init_vector mass(1,nages)
  init_matrix sobs_survey_at_age(1,nyrs,1,nages)
  matrix obs_survey_at_age(1,snyrs,1,nages)
  matrix OSP(1,snyrs,1,nages)
  vector OTSC(1,snyrs)
 LOCAL_CALCS
  for (i=1;i<=snyrs;i++)
  {
    effort(i)=seffort(i);
    obs_survey_at_age(i)=sobs_survey_at_age(i);
    obs_catchsd_at_age(i)=sobs_catchsd_at_age(i);
    OTSC(i)=sum(obs_survey_at_age(i));
    OSP(i)=obs_survey_at_age(i)/OTSC(i);
  }
 END_CALCS
  init_vector ssurvey_cpue(1,nyrs)
  vector survey_cpue(1,snyrs)
  init_vector maturity(1,nages)
  init_vector scatch_biomass(1,nyrs)
  vector catch_biomass(1,snyrs)
  matrix catchw(1,snyrs,1,nages)
 LOCAL_CALCS
  for (i=1;i<=snyrs;i++)
  {
    survey_cpue(i)=ssurvey_cpue(i);
    catch_biomass(i)=scatch_biomass(i);
  }

INITIALIZATION_SECTION
  log_q -1
  log_avginit 4.0
  log_avgrec 4.0
  M 0.225

PARAMETER_SECTION
  init_number log_q(1)
  init_bounded_vector log_q_devs(1,snyrs-1,-1.,1.,4)
  init_number log_surv_q(1)
  init_bounded_vector log_surv_q_devs(1,snyrs-1,-1.,1.,4)
  init_bounded_number M(.03,1.8,phase_M)
  init_bounded_vector sel_coffs(1,nng,-8.,8.,3)
  init_bounded_vector surv_sel_coffs(1,nng,-8.,8.,3)
  init_bounded_matrix sel_devs(1,dim_sel,1,nng,-5.,5.,5)
  init_bounded_matrix surv_sel_devs(1,dim_sel,1,nng,-5.,5.,5)
  init_number log_avginit(1)
  init_number log_avgrec(1)
  init_bounded_dev_vector log_recdevs(1,snyrs,-15.,15.,2)
  init_bounded_dev_vector log_initdevs(2,nages,-15.,15.,2)
  init_bounded_vector effort_devs(1,snyrs,-5.,5.,3)
  init_bounded_vector surv_effort_devs(1,snyrs,-5.,5.,3)
  init_bounded_vector M_devs(1,snyrs,-5.,5.,7)
  matrix log_fish_sel(1,snyrs,1,nages);
  vector lcatch(1,snyrs)
  vector lsurvcatch(1,snyrs)
  matrix log_surv_sel(1,snyrs,1,nages);
  vector log_recruits(1,snyrs);
  vector log_initage(2,nages);
  matrix F(1,snyrs,1,nages)
  matrix Z(1,snyrs,1,nages)
  matrix S(1,snyrs,1,nages)
  matrix N(1,snyrs,1,nages)
  matrix C(1,snyrs,1,nages)
  matrix EC(1,snyrs,1,nages)
  matrix ESP(1,snyrs,1,nages)
  matrix EP(1,snyrs,1,nages)
  vector ETC(1,snyrs)
  vector ETSC(1,snyrs)
  matrix E_surv_age(1,snyrs,1,nages)
  objective_function_value f
  number stock
  number recsum
  number fcatch
  number sel_tmp
  number sel_shape
  number feffort
  number initsum
  number avgsel;
  number avgsurvsel;
  sdreport_number avg_F
  sdreport_vector predicted_N(2,nages)
  sdreport_vector available_biomass(1,snyrs)
  sdreport_vector ratio_N(2,nages)
  sdreport_number pred_B
  sdreport_number log_avg_M
  sdreport_number available_biomass_ratio
  
PRELIMINARY_CALCS_SECTION
  catchw=square(1/(obs_catchsd_at_age+.0001));

RUNTIME_SECTION
  maximum_function_evaluations 100,100,200,300,500,4000

PROCEDURE_SECTION
  get_mortality_and_survivial_rates();
  get_numbers_at_age();
  get_catch_at_age();
  evaluate_the_objective_function();

FUNCTION get_mortality_and_survivial_rates
//---Calculate the selectivity --------------------------
  avgsel=log(mean(exp(sel_coffs)));
  avgsurvsel=log(mean(exp(surv_sel_coffs)));

  for (int j=1;j<=nng;j++)
  {
    log_fish_sel(1,j)=sel_coffs(j);
    log_surv_sel(1,j)=surv_sel_coffs(j);
  }
  //--The selectivity is the same for the last last_age_sel_group age classes--------------
  for (j=nng+1;j<=nages;j++)
  {
    log_fish_sel(1,j)=sel_coffs(nng);
    log_surv_sel(1,j)=surv_sel_coffs(nng);
  }

  int ii;
  log_fish_sel(1)-=log(mean(exp(log_fish_sel(1))));

  if (active(sel_devs))
  {
    ii=1;
    for (int i=1;i<snyrs;i++)
    {
      if (!(i%group_num))
      {
        for (int j=1;j<=nng;j++)
        { 
          log_fish_sel(i+1,j)=log_fish_sel(i,j)+sel_devs(ii,j);
        }      
        for (j=nng+1;j<=nages;j++)
        {
          log_fish_sel(i+1,j)=log_fish_sel(i+1,nng);
        }      
        ii++;
        log_fish_sel(i+1)-=log(mean(exp(log_fish_sel(i+1))));
      }
      else
      {
        for (int j=1;j<=nng;j++)
        { 
          log_fish_sel(i+1,j)=log_fish_sel(i,j);
        }
        for (j=nng+1;j<=nages;j++)
        {
          log_fish_sel(i+1,j)=log_fish_sel(i+1,nng);
        }      
        log_fish_sel(i+1)-=log(mean(exp(log_fish_sel(i+1))));
      }
    }
  }
  else
  {
    for (int i=1;i<snyrs;i++)
    {
      for (int j=1;j<=nng;j++)
      { 
        log_fish_sel(i+1,j)=log_fish_sel(i,j);
      }
      for (j=nng+1;j<=nages;j++)
      {
        log_fish_sel(i+1,j)=log_fish_sel(i+1,nng);
      }      
    }
  }


  log_surv_sel(1)-=log(mean(exp(log_surv_sel(1))));
  if (active(surv_sel_devs))
  {
    ii=1;
    for (int i=1;i<snyrs;i++)
    {
      if (!(i%group_num))
      {
        for (int j=1;j<=nng;j++)
        { 
          log_surv_sel(i+1,j)=log_surv_sel(i,j)+surv_sel_devs(ii,j);
        }      
        for (j=nng+1;j<=nages;j++)
        {
          log_surv_sel(i+1,j)=log_surv_sel(i+1,nng);
        }      
        ii++;
        log_surv_sel(i+1)-=log(mean(exp(log_surv_sel(i+1))));
      }
      else
      {
        for (int j=1;j<=nng;j++)
        { 
          log_surv_sel(i+1,j)=log_surv_sel(i,j);
        }
        for (j=nng+1;j<=nages;j++)
        {
          log_surv_sel(i+1,j)=log_surv_sel(i+1,nng);
        }      
        log_surv_sel(i+1)-=log(mean(exp(log_surv_sel(i+1))));
      }
    }
  }
  else
  {
    for (int i=1;i<snyrs;i++)
    {
      for (int j=1;j<=nng;j++)
      { 
        log_surv_sel(i+1,j)=log_surv_sel(i,j);
      }
      for (j=nng+1;j<=nages;j++)
      {
        log_surv_sel(i+1,j)=log_surv_sel(i+1,nng);
      }      
    }
  }

  lcatch(1)=log_q;
  if (active(log_q_devs))
  {
    for (int i=1;i<snyrs;i++)
    {
      lcatch(i+1)=lcatch(i)+log_q_devs(i);
    }
  }
  else
  {
    for (int i=1;i<snyrs;i++)
    {
      lcatch(i+1)=lcatch(i);
    }
  }

  lsurvcatch(1)=log_surv_q;
  if (active(log_surv_q_devs))
  {
    for (int i=1;i<snyrs;i++)
    {
      lsurvcatch(i+1)=lsurvcatch(i)+log_surv_q_devs(i);
    }
  }
  else
  {
    for (int i=1;i<snyrs;i++)
    {
      lsurvcatch(i+1)=lsurvcatch(i);
    }
  }


//--This is the same as F(i,j)=exp(q)*effort(i)*exp(log_sel(i,j));----------
  for (int i=1;i<=snyrs;i++)
  {
    F(i)=mfexp(lcatch(i))*effort(i)*mfexp(log_fish_sel(i));
  }

  if (active(effort_devs))
  {
    for (int i=1;i<=snyrs;i++)
    {
      F(i)=F(i)*exp(effort_devs(i));
    }
  }

//--Get the total mortality-----------

  if (active(M_devs))
  {
    Z=F;
    for (int i=1;i<=snyrs;i++)
    {
      Z(i)=Z(i) + M*exp(M_devs(i));
    }
  }
  else
  {
    Z=F+M;
  }

//--Get the survival rate-----------
//
  S=mfexp(-1.0*Z);

FUNCTION get_numbers_at_age
//--Initial age composition here-------------------
  log_initage=log_initdevs+log_avginit;
  for (int j=2;j<=nages;j++)
  {
    N(1,j)=mfexp(log_initage(j));
  }
   log_recruits = log_avgrec + log_recdevs;
//--Subsequent Recruitment here-------------------
  for (int i=1;i<snyrs;i++)
  {
    N(i,1)= exp(log_recruits(i));
    for (j=1;j<nages;j++)
    {
      N(i+1,j+1)=N(i,j)*S(i,j);
    }
  //--Cumulate oldest age bin here-------------------
      N(i+1,nages)+=N(i,nages)*S(i,nages);
  }
//--Recruitment in last year -----------
  N(snyrs,1)= exp(log_recruits(snyrs));

//--Calculated predicted numbers at age for next year--------------------
  for (j=1;j<nages;j++)
  {
    predicted_N(j+1)=N(snyrs,j)*S(snyrs,j);
    ratio_N(j+1)=predicted_N(j+1)/N(1,j+1);
  }
//--Calculate predicted Biomass for next year for reporting confidence information--------
  pred_B=(predicted_N * relwt);


FUNCTION get_catch_at_age
  C=elem_prod(elem_div(F,Z),elem_prod(1.-S,N));
  for (int i=1;i<=snyrs;i++)
  {
    {
      EC(i)=C(i)*agerr;
      if (!active(surv_effort_devs))
      {
        E_surv_age(i)=mfexp(lsurvcatch(i))*
          elem_prod(mfexp(log_surv_sel(i)),N(i)*agerr);
      }
      else
      {
        E_surv_age(i)=mfexp(surv_effort_devs(i)+lsurvcatch(i))*
          elem_prod(mfexp(log_surv_sel(i)),N(i)*agerr);
      }
    }
  }
  for (i=1;i<=snyrs;i++)
  {
    ETC(i)=sum(EC(i));
    EP(i)=EC(i)/ETC(i);
  }
  for (i=1;i<=snyrs;i++)
  {
    ETSC(i)=sum(E_surv_age(i));
    ESP(i)=E_surv_age(i)/ETSC(i);
  }


FUNCTION evaluate_the_objective_function

  if (sd_phase())
  {
    dvar_matrix fish_sel=mfexp(log_fish_sel);
    for (int i=1;i<=snyrs;i++)
    {
      dvar_vector real_sel=fish_sel(i)/max(fish_sel(i));
      dvar_vector avbio= elem_prod(elem_prod(N(i),real_sel),mass);
      available_biomass(i)=sum(avbio);
    }
    available_biomass_ratio= available_biomass(snyrs)/available_biomass(1);
  }
//--Regularize the solution with penalties
  if (current_phase()<3)
  {
    f+=100*norm2(log_recdevs)+ 100*norm2(log_initdevs);
    f+=100*norm2(first_difference(first_difference(log_recdevs)))
     +100*norm2(first_difference(first_difference(log_initdevs)));
    f+=1000*norm2(first_difference(first_difference(log_fish_sel(1))))
     +1000*norm2(first_difference(first_difference(log_surv_sel(1))));
  }
  f+=lambda1*norm2(log_recdevs)+ lambda2*norm2(log_initdevs);
  f+=lambda1*norm2(first_difference(first_difference(log_recdevs)))
   +lambda2*norm2(first_difference(first_difference(log_initdevs)));

  if (control_flags(1))
  {
    for (int i=1;i<=snyrs;i++)
    {
      for (int j=1;j<nages;j++)
      {
        if (log_fish_sel(i,j)>log_fish_sel(i,j+1))
        {
          f+=control_flags(1)*square(log_fish_sel(i,j)-log_fish_sel(i,j+1));
        }
      }
    }
  }

  if (control_flags(2))
  {
    for (int i=1;i<=snyrs;i++)
    {
      for (int j=1;j<nages;j++)
      {
        if (log_surv_sel(i,j)>log_surv_sel(i,j+1))
        {
          f+=control_flags(2)*square(log_surv_sel(i,j)-log_surv_sel(i,j+1));
        }
      }
    }
  }

  if (active(surv_sel_devs))
  {
    f+=10./dim_sel*norm2(first_difference(first_difference(log_surv_sel(1))));
    for (int i=1;i<snyrs;i++)
    {
      if (!(i%group_num))
      {
        f+=10./dim_sel*
          norm2(first_difference(first_difference(log_surv_sel(i+1))));
      }
    }
  }
  else
  {
    f+=10*norm2(first_difference(first_difference(log_surv_sel(1))));
  }

  if (active(sel_devs))
  {
    f+=10./dim_sel*norm2(first_difference(first_difference(log_fish_sel(1))));
    for (int i=1;i<snyrs;i++)
    {
      if (!(i%group_num))
      {
        f+=10./dim_sel*
          norm2(first_difference(first_difference(log_fish_sel(i+1))));
      }
    }
  }
  else
  {
    f+=10*norm2(first_difference(first_difference(log_fish_sel(1))));
  }

  if (active(M_devs))
  {
    f+=25*norm2(M_devs);
  }
  avg_F=sum(F)/double(size_count(F));
  if (current_phase()>5)
  {
    // a very small penalty on the average fishing mortality
    f+= .001*square(log(avg_F/.2));
  }
  else if (current_phase()>3)
  {
    f+= 10.*square(log(avg_F/.2));
  }
  else 
  {
    f+= 1000.*square(log(avg_F/.2));
  }
  switch (current_phase())
  {
  case 1:
    f+= norm2(elem_div(square(EC-obs_catch_at_age),.01+EC))
      +norm2(elem_div(square(E_surv_age-obs_survey_at_age),.01+E_surv_age));
    break;
  case 2:
  case 3:
  case 4:
    f+=0.5*double(2*size_count(C))
      * log( sum(elem_div(square(EC-obs_catch_at_age),.02+EC))
      + sum(elem_div(square(E_surv_age-obs_survey_at_age),.02+E_surv_age)))
    + 10*norm2(effort_devs)
    + 10*norm2(surv_effort_devs)
    + 50*norm2(log_q_devs)
    + 250*norm2(log_surv_q_devs)
    +100./group_num*lambda4*norm2(sel_devs)
    +250./group_num*lambda4*norm2(surv_sel_devs);
    break;
  default:
    double rf=.1/nng;
    dvar_matrix vc=rf+2.*elem_prod(EP,1-EP);
    dvar_matrix vs=rf+2.*elem_prod(ESP,1-ESP);
    dvar_matrix lc=exp(-c_size*elem_div(square(EP-OP),vc));
    dvar_matrix ls=exp(-s_size*elem_div(square(ESP-OSP),vs));
    //cout << "lc = " << endl << lc << endl << endl;
    //dvar_matrix mtmp= c_size*lc+.01;
    //cout << "c_size*lc+.01= " << endl 
    //  << mtmp << endl << endl;
    //int rmin=mtmp.rowmin();
    //int rmax=mtmp.rowmax();
    //int cmin=mtmp.colmin();
    //int cmax=mtmp.colmax();
    //cout << "log(c_size*lc+.01) by element = " << endl;
    //cout << "log(c_size*lc+.01) = " << endl;
    //cout << log(mtmp) << endl << endl;
    f-=sum(log(lc+.01));
    f-=sum(log(ls+.01));
    f+=0.5*sum(log(vc));
    f+=0.5*sum(log(vs));
    f+= 100.*norm2(log(OTC+.01)-log(ETC+.01));
    f+= 100.*norm2(log(OTSC+.01)-log(ETSC+.01));
    f+= 10*norm2(effort_devs)
    + 10*norm2(surv_effort_devs)
    + 50*norm2(log_q_devs)
    + 250*norm2(log_surv_q_devs)
    +100./group_num*lambda4*norm2(sel_devs)
    +250./group_num*lambda4*norm2(surv_sel_devs);
    break;
  }

  f+=10.*square(avgsel);
  f+=10.*square(avgsurvsel);
//--Prior for M----------------------
   f+= 5.*square(log(M/0.20));

   if (sd_phase) 
   {
     log_avg_M=mean(M*exp(M_devs));
   }

//  if (active(effort_devs))
//  {
//    f+= 10*norm2(effort_devs);
//  }
//
//  if (active(sel_devs))
//  {
//    f+= 50*norm2(sel_devs);
//  }
//
//  if (active(log_q_devs))
//  {
//    f+= 50*norm2(log_q_devs);
//  }
//

REPORT_SECTION
  report << "Estimated numbers of fish " << endl;
  report << setfixed << setprecision(1) << setw(6) << N << endl;
  report << "Fishery Catchability" << endl;
  report << setscientific << setprecision(3) << setw(8) << 
            mfexp(lcatch) << endl;
  report << "Fishery Selectivities" << endl;

  dvar_matrix fish_sel=mfexp(log_fish_sel);
  report << setfixed << setprecision(3) << setw(6) << log_fish_sel << endl;

  double qmax=max(value(lcatch));

  // init_vector mass(1,nages)
  report << "available biomass by year" << endl;

  for (int i=1;i<=snyrs;i++)
  {
    dvar_vector real_sel=fish_sel(i)/max(fish_sel(i));
    dvar_vector avbio= elem_prod(elem_prod(N(i),real_sel),mass);
    available_biomass(i)=sum(avbio);
    report << setscientific << setprecision(3) << available_biomass(i) << "  "
           <<  avbio  << endl;
  }
  available_biomass_ratio= available_biomass(snyrs)/available_biomass(1);
  report << "available_biomass_ratio last_year/first_year = " 
       << available_biomass_ratio << endl;


  report << "Fishery Catchability by age" << endl;
  for (i=1;i<=snyrs;i++)
  {
    report << setfixed << setprecision(3) << setw(6)  
           << mfexp(lcatch(i)+log_fish_sel(i)-qmax) << endl;
  }
  report << "Survey Catchability by age" << endl;
  double sqmax=max(value(lsurvcatch));
  for (i=1;i<=snyrs;i++)
  {
    report << setfixed << setprecision(3) << setw(6)  
           << mfexp(lsurvcatch(i)+log_surv_sel(i)-sqmax) << endl;
  }
  report << "Survey Catchability" << endl;
  report << setfixed << setprecision(3) << setw(6) << 
            mfexp(lsurvcatch) << endl;
  report << "Survey Selectivities" << endl;
  report << setfixed << setprecision(3) << setw(6) << 
            mfexp(log_surv_sel) << endl;
  report << "Estimated numbers in catch " << endl;
  report << EC << endl;
  report << "Observed numbers in catch " << endl;
  report << obs_catch_at_age << endl; 
  report << "Estimated average natural mortality " << endl;
  report << mean(M*exp(M_devs)) << endl;
  report << "Estimated natural mortality " << endl;
  report << M*exp(M_devs) << endl;
  report << "Estimated fishing mortality " << endl;
  report << F << endl; 
  report << "Estimated survey selectivity " << endl;
  report << mfexp(log_surv_sel) << endl; 
  report << "Estimated Survey at age " << endl;
  report << E_surv_age << endl;
  report << "Observed survey at age" << endl;
  report << obs_survey_at_age << endl; 
  report << "Maturity" << endl;
  report << maturity << endl;
  report << "ObsEffort predicted SurvBiom CatchBiom" << endl;
      for (i=1;i<=snyrs;i++)
	      {
  report << effort(i)<<" " ;
  report << effort(i)*mfexp(effort_devs(i)) << " " ;
  report << survey_cpue(i) << " ";
  report << catch_biomass(i) << " "<<endl;
           }
  report << "Predicted N Yr30 (age2+)" << endl;
  report << predicted_N << endl;
  report << "Fishery Catch_at_age residuals" << endl;
  report << elem_div(EC-obs_catch_at_age,sqrt(.02+EC)) << endl << endl;
  report << "Survey catch_at_age residuals" << endl;
  report << elem_div(E_surv_age-obs_survey_at_age,sqrt(.02+E_surv_age))
         << endl << endl;
  report << "fishery total catch residuals " << endl
         << log(OTC+.01)-log(ETC+.01) << endl;
  report << "survey total catch residuals " << endl
         << log(OTSC+.01)-log(ETSC+.01) << endl;
TOP_OF_MAIN_SECTION
  gradient_structure::set_MAX_NVAR_OFFSET(1000);
  gradient_structure::set_GRADSTACK_BUFFER_SIZE(100000);
  gradient_structure::set_CMPDIF_BUFFER_SIZE(1000000);