// Copyright (c) 2008, 2009, 2010 Regents of the University of California.
//
// ADModelbuilder and associated libraries and documentations are
// provided under the general terms of the "BSD" license.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2.  Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3.  Neither the name of the  University of California, Otter Research,
// nor the ADMB Foundation nor the names of its contributors may be used
// to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

// MODIFIED BY STEVE MARTELL, March 15, 2012 to demonstrate the use of
// contributed functions from statsLib.h


DATA_SECTION
  init_int nyrs  
  init_int nages
  init_int sel_type
  init_int nodes
  vector age(1,nages);
  !! age.fill_seqadd(1,1);
  init_matrix obs_catch_at_age(1,nyrs,1,nages)
  init_vector effort(1,nyrs)
  init_number M
  vector relwt(2,nages);
INITIALIZATION_SECTION
  log_q -1
  log_popscale 5
PARAMETER_SECTION
  init_number log_q(1)
  init_number log_popscale(1)
  init_bounded_dev_vector log_sel_coff(1,nodes,-15.,15.,2)
  init_bounded_dev_vector log_relpop(1,nyrs+nages-1,-15.,15.,2)
  init_bounded_dev_vector effort_devs(1,nyrs,-5.,5.,3)
  vector sel_par(1,nodes);
  vector log_sel(1,nages)
  vector log_initpop(1,nyrs+nages-1);
  matrix F(1,nyrs,1,nages)
  matrix Z(1,nyrs,1,nages)
  matrix S(1,nyrs,1,nages)
  matrix N(1,nyrs,1,nages)
  matrix C(1,nyrs,1,nages)
  objective_function_value f
  number recsum
  number initsum
  sdreport_number avg_F
  sdreport_vector predicted_N(2,nages)
  sdreport_vector ratio_N(2,nages)
  // changed from the manual because adjusted likelihood routine doesn't
  // work
  likeprof_number pred_B
PRELIMINARY_CALCS_SECTION
  // this is just to ``invent'' some relative average
  // weight at age numbers
  relwt.fill_seqadd(1.,1.);
  relwt=pow(relwt,.5);
  relwt/=max(relwt);
PROCEDURE_SECTION
  // example of using FUNCTION to structure the procedure section
  get_mortality_and_survivial_rates();

  get_numbers_at_age();

  get_catch_at_age();

  //evaluate_the_objective_function();
  calc_objective_function();

FUNCTION get_mortality_and_survivial_rates
  int i, j;
  // calculate the selectivity from the sel_coffs
  // sel_type is a flag for which selectivity function to use
  // 1) age-specific sel_coffs
  // 2) cubic_spline interpolation estimating n nodes
  // 3) 2dcubic spline.
  switch(sel_type)
  {
  	case 1:
  		for (j=1;j<=nodes;j++)
	  {
	    log_sel(j)=log_sel_coff(j);
	  }
	  // the selectivity is the same for the last two age classes
	  log_sel(nodes,nages)=log_sel_coff(nodes);
  	break;
	case 2:
	  sel_par=log_sel_coff;
	  log_sel = cubic_spline(sel_par,age);
	break;
	
  }


  

  // This is the same as F(i,j)=exp(q)*effert(i)*exp(log_sel(j));
  F=outer_prod(mfexp(log_q)*effort,mfexp(log_sel));
  if (active(effort_devs))
  {
    for (i=1;i<=nyrs;i++)
    {
      F(i)=F(i)*exp(effort_devs(i));
    }
  }
  // get the total mortality
  Z=F+M;
  // get the survival rate
  S=mfexp(-1.0*Z);

FUNCTION get_numbers_at_age
  int i, j;
  log_initpop=log_relpop+log_popscale;
  for (i=1;i<=nyrs;i++)
  {
    N(i,1)=mfexp(log_initpop(i));
  }
  for (j=2;j<=nages;j++)
  {
    N(1,j)=mfexp(log_initpop(nyrs+j-1));
  }
  for (i=1;i<nyrs;i++)
  {
    for (j=1;j<nages;j++)
    {
      N(i+1,j+1)=N(i,j)*S(i,j);
    }
  }
  // calculated predicted numbers at age for next year
  for (j=1;j<nages;j++)
  {
    predicted_N(j+1)=N(nyrs,j)*S(nyrs,j);
    ratio_N(j+1)=predicted_N(j+1)/N(1,j+1);
  }
  // calculated predicted Biomass for next year for
  // adjusted profile likelihood
  pred_B=(predicted_N * relwt);

FUNCTION get_catch_at_age
  C=elem_prod(elem_div(F,Z),elem_prod(1.-S,N));

FUNCTION calc_objective_function
	/*
	Modified by Steve Martell to show the use of the functions
	in the statsLib.h
	*/
	dvar_vector fvec(1,4);
	fvec.initialize();
	avg_F=sum(F)/double(size_count(F));
	fvec(1) = 0.01*norm2(log_relpop);
	if (last_phase())
	{
		fvec(2) = 0.001*square(log(avg_F/0.2));
	}
	else
	{
		fvec(2) = 1000.*square(log(avg_F/0.2));
	}
	fvec(3) = dnorm(obs_catch_at_age,C);
	
	//the std in the following expression is based
	//on the wt of 0.1, which is std=sqrt(1/(2*wt))
	fvec(4) = dnorm(effort_devs,2.236068);
	f = sum(fvec);
	

FUNCTION evaluate_the_objective_function
  // penalty functions to ``regularize '' the solution
  f+=.01*norm2(log_relpop);
  avg_F=sum(F)/double(size_count(F));
  if (last_phase())
  {
    // a very small penalty on the average fishing mortality
    f+= .001*square(log(avg_F/.2));
  }
  else
  {
    f+= 1000.*square(log(avg_F/.2));
  }
  f+=0.5*double(size_count(C)+size_count(effort_devs))
    * log( sum(elem_div(square(C-obs_catch_at_age),.01+C))
    + 0.1*norm2(effort_devs));

GLOBALS_SECTION
  #include <statsLib.h>

REPORT_SECTION
  report << "Estimated numbers of fish " << endl;
  report << N << endl;
  report << "Estimated numbers in catch " << endl;
  report << C << endl;
  report << "Observed numbers in catch " << endl;
  report << obs_catch_at_age << endl; 
  report << "Estimated fishing mortality " << endl;
  report << F << endl; 
  report << "Selectivity"  << endl;
  report << exp(log_sel) <<endl;
