// File: herring.tpl. 
// Nedstrippet versjon til bruk ifm torsk på HI, 2010.
//
// Opprinnelig:
// Date: 29.09.02
// Author: Hans Julius Skaug (skaug@imr.no)
//
// Description: ADMB implementation of Sigurd's herring assessment. 
//


DATA_SECTION
  !! cout << "***************************************************************" << endl;
  !! cout << "* VERSION WITH SCALING, INPUT N0 WILL BE MULTIPLIED WITH 1000 *" << endl;
  !! cout << "***************************************************************" << endl;
  !! cout << endl;

  // Fixed parameters in the model
  init_int 	m			// Number of cohorts

  !! cout << "Number of cohorts: " << m << endl;

  // Starting and ending year of cohort observation. Cohort is born in the
  // starting year.
  init_vector n_start_year(1,m)
  init_vector n_end_year(1,m)

  // Keep track of largest index of each cohort (it is not 
  // a given that cohorts are followed for the same amount of years)
  // note that this index starts at zero
  vector n_largest_year_index(1,m)
  !! n_largest_year_index = n_end_year - n_start_year;
  !! cout << "n_largest_year_index " << n_largest_year_index << endl;

  // MATRIX WITH CATCH DATA, one cohort per row, one year per column
  // NOTE: The size of C is such that catch in the last year 
  // (which won't be used for any computation) MUST be present!!
  // Setting it to zero is fine
  init_matrix C(1,m,0,n_largest_year_index)

  // Survey data
  init_int 	n_surv			// Number surveys

  !! cout << "Number of surveys: " << n_surv << endl;

  init_int 	n_I			// Number survey indices
  init_matrix I(1,n_I,1,4)		// "cohort","age","survey","I"
  // end of survey data

  // Number of calendar years our model spans. Computed from start and end years in data.
  int n_span
  !! n_span = max(n_end_year) - min(n_start_year) + 1;

  !! cout << "n_span = max(n_end_year) - min(n_start_year) + 1 = " << n_span << endl; 
  !! cout << "Ignoring survey data for which there is no catch data!" << endl;


PARAMETER_SECTION

  // Recruitement 
  init_bounded_vector   N0(1,m,0.01,10000) 

  // Population trajectory, one cohort per row, one year per column.
  // (Entry i,j is year j in the life span of cohort i, which is calendar year
  // n_start_year(i)+j.) The width corresponds to the largest
  // life span we follow, so not all columns will be used unless all spans are the same.
  matrix 		N(1,m,0,n_span-1)

  // Catchability (survey specific)
  init_bounded_number q0(.05,1.0)
 
  init_number	q_alpha
  init_number	q_beta

  init_bounded_number 	logs(-5,3)		// Log standard deviation of survey errors
  number 	s				// Standard deviation of survey errors
  init_bounded_number   M(0,.5)                 // Mortality
 
  number tmp  
  
  objective_function_value l

PRELIMINARY_CALCS_SECTION
	
  


GLOBALS_SECTION
  
  // Helper function to keep population nonnegative
  #include <fvar.hpp>
  dvariable posfun(_CONST dvariable&x,const double eps)
  {
    if (x>=eps) 
    {
      return x;
    }
    else
    {
      dvariable y=1.0-x/eps;
      return eps*(1./(1.+y+y*y));
    }
  }

// Removed by LF 12032010
// #include <df1b2fun.h>

  //#include <df1b2loc.h>


PROCEDURE_SECTION

  int i,j;

  s= exp(logs);
  
  // initialize likelihood value
  l = 0;

  // Initialize N
  for(i=1;i<=(int) m;i++)
    N(i,0) = 1000*N0(i);


  // Compute population trajectory, note use of helper function posfun
  for(i=1;i<=m;i++) {
    for(j=1;j<=(int) n_end_year(i)-n_start_year(i);j++) {

         N(i,j) = posfun(N(i,j-1)-C(i,j-1),.01)*exp(-M);
     
    }
  }


  // Likelihood contributions from surveys
  for(i=1;i<=(int) n_I;i++)
  {
    // NEW --- LF, 19 Sept 2009 
    // Only include survey index if it comes from a year where N actually has a meaningful value
    // That is, that it has been recursively computed
    if (  I(i,2) <= n_largest_year_index(I(i,1))  ) {
      // New for this model - sigmoid function for q
      //tmp = 0.0;
      tmp = (log(I(i,4))  - log(q0*(exp(q_alpha*I(i,2)+q_beta)/(1+exp(q_alpha*I(i,2)+q_beta))) *N(I(i,1),I(i,2)))) /s;
      l += -.5*tmp*tmp - logs;
    }
  }

  // Reverse objective function sign since ADMB does minimization
  l *= -1;

  cout << setprecision(4);


REPORT_SECTION
  report << setprecision(4);
  //cout << "l="<< l << " N0=" << N0 << " q="<< q << " M=" << M  << " s=" << s << endl;

  report << N << endl;
TOP_OF_MAIN_SECTION
  arrmblsize = 50000000L;
  gradient_structure::set_GRADSTACK_BUFFER_SIZE(1000000);
  gradient_structure::set_CMPDIF_BUFFER_SIZE(2000000);
  gradient_structure::set_MAX_NVAR_OFFSET(40205);