AdmArwiki

Most of the models presented so far have followed one cohort in the first year under consideration, two in the second, three in the third, and so on. In addition, the models following ages 3 and up have relied on a preprocessing of the input data in order to ignore age groups zero through two. We now present a model which follows a fixed number of cohorts every year, without the need to preprocess the data.

The properties of the model are given in the following table:

NOTE that N is structured differently than for the models in the other sections. The same goes for catch and survey data in the input. Essentially, the diagonals of this N correspond to the N used in other models, and similarly for the catch matrix.

The model is implemented in the following code. An example of a .dat file can be found here.

cod.tpl

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