Model following ALL cohorts

Description

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:

 Input data (example) start_year (1994), end_year (2001), minage (3), maxage (12). Corresponding catch (C) and survey data (I). Matrix N, indexed as N(i,j) = cohort aged j in year i. First row and column of N, catchability q, survey standard deviation s, mortality M. log(I) ~ N( log(qN), s). N(i,j) = [ N(i-1,j-1) - C(i-1,j-1) ] exp(-M). 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; 