// Cod population model written by Lennart Frimannslund, lennartfr@imr.no, November 2010.
//
// Models ALL cohorts "alive" at any given time, typically the 11 cohorts in age group 3-13.
//
//


DATA_SECTION

  init_int start_year
  init_int end_year

  init_int min_age
  init_int max_age

  // Catch
  init_matrix C(start_year,end_year,min_age,max_age)
  // Survey indices
  init_matrix I(start_year,end_year,min_age,max_age)

  number catch_scaling
  number survey_scaling
  number N_scaling
  number small_value

  !! catch_scaling = 1e-3;
  !! survey_scaling = 1e0;
  !! N_scaling = 1e3;
  !! small_value = 1e-5;
  
PARAMETER_SECTION
  
  init_bounded_vector N_firstcol(start_year,end_year,0.01,10000)
  init_bounded_vector N_firstrow(min_age+1,max_age,0.01,10000)

  init_bounded_vector q(1,1,0.05,1,2)
  init_bounded_number logs(-5,3)
  init_bounded_number M(0,0.5)

  sdreport_matrix N(start_year,end_year,min_age,max_age)
  number s  
  number temp

  sdreport_vector N46(start_year,end_year)
  sdreport_vector N7plus(start_year,end_year)

  objective_function_value l

PRELIMINARY_CALCS_SECTION

  cout << "N_firstcol = " << N_firstcol << endl;
  cout << "N_firsrow  = " << N_firstrow << endl;
  cout << "q = " << q << endl;
  cout << "logs = " << logs << endl;
  cout << "M = " << M << endl;


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));
    }
  }


PROCEDURE_SECTION

  // Debug
  // cout << "Hey!" << endl;

  // counters
  int i,j;
  temp = 0.0;

  // Initialize l
  l = 0;
   
  //
  s = exp(logs);

  // Initialize population trajectory
  for (i=min_age+1;i<=max_age;i++) {
     N(start_year,i) = N_scaling * N_firstrow(i);
     //cout << N(i,start_year) << " ";
  }
  //cout << endl;

  // Debug
  // cout << "Hey!" << endl;

  for (i=start_year;i<=end_year;i++) {
    N(i,min_age) = N_scaling*N_firstcol(i);
    //cout << N(min_age,i) << " ";
  }   

  //cout << endl;  
  // Debug
  // cout << "Hey!" << endl;

  
  // Compute population trajectory with recursive formula
  for (i=start_year+1;i<=end_year;i++) {
    for (j=min_age+1;j<=max_age;j++) {
     
      N(i,j) = posfun(  N(i-1,j-1) - C(i-1,j-1)*catch_scaling, 0.01 ) * exp(-M);
      //cout << N(i,j) << " ";
    }
    //cout << endl;
  }

  // Debug
  // cout << "Hey!" << endl;

  // Contribution from survey misfit
  for (i=start_year;i<=end_year;i++) {
    for (j=min_age;j<=max_age;j++) {

       if (I(i,j) > small_value) {
         temp = ( log(I(i,j)*survey_scaling) - log(q(1)*N(i,j)) ) / s;
         l += -0.5 * temp * temp  - logs;
       }

    }
  }

  // Debug
  // cout << "Hey!" << endl;

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

  // Calculate plus classes
  for (i=start_year;i<=end_year;i++) {
    N46(i)    = 0.0;
    N7plus(i) = 0.0;
    for (j=4;j<=6;j++) {
        N46(i) += N(i,j);
    }

    for (j=7;j<=max_age;j++) {
        N7plus(i) += N(i,j);
    }

 
  } 

  // Debug
  // cout << "Ho!" << endl;

REPORT_SECTION
  
  report << N << endl;
  // report << endl << "# C " << endl;
  // report << C << endl;
  // report << endl << "# I " << endl;
  // report << I << 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);