DATA_SECTION

  init_int nobs			// Number of observations
  !! TRACE(nobs)
  init_int nroute
  !! TRACE(nroute)
  init_ivector route(1,nobs)    // vector of routes
  !! TRACE(route)
  init_ivector obs_per_route(1,nroute)
  !! TRACE(obs_per_route)
  init_int S			// Number of response categories
  int S1			// For variable allocation
  !! S1=S-1;			
  int S2			// For variable allocation
  !! S2=S-2;			
  init_int p;			// Number of fixed effects
  !! TRACE(p)
  init_imatrix y(1,nroute,1,obs_per_route) // Response vector
  !! TRACE(y)
  init_3darray X(1,nroute,1,obs_per_route,1,p)	// Fixed effects design matrix
  !! TRACE(X)
  init_matrix year(1,nroute,1,obs_per_route)  // Year for each observation
  !! TRACE(year)
  init_matrix dist(1,nroute,1,nroute)
  !! TRACE(dist)

PARAMETER_SECTION

  init_bounded_vector beta(1,p,-5,5,1) 		// Fixed effects
  init_bounded_vector tmpk(1,S2,-3,1,1)		// Variables underlying kappa
  init_bounded_number a(0.01,3.0,2)

  random_effects_vector b_int(1,nroute,2)
  random_effects_vector b_year(1,nroute,2)
  normal_prior M(b_int)

  objective_function_value g

  vector kappa(1,S1)

PRELIMINARY_CALCS_SECTION
  cout << setprecision(4);

GLOBALS_SECTION

  #include <df1b2fun.h>
  #include <fstream>
  #include <math.h>
  ofstream clogf("program.log");
  #define TRACE(obj) clogf<<"line "<<__LINE__<<", file "<<__FILE__<<", "<<#obj" =\n" \
                        <<obj<<endl<<endl;
  double pi = M_PI; 

PROCEDURE_SECTION

  g = 0;

  for (int route = 1; route <= nroute; route++)
    sepfun(route, beta, tmpk, kappa, b_int(route), b_year(route));

SEPARABLE_FUNCTION void sepfun(int route, const dvar_vector& beta, const dvar_vector& tmpk, const dvar_vector& kappa, const dvariable& b_int, const dvariable& b_year)

  g += 0.5*log(2*pi) + 0.5*square(b_int) + 0.5*square(b_year);

  // Sets up kappa from underlying variable tmpk
  kappa(1) = 0;  // Shift route intercepts by b_int
  for (int s = 2; s <= S1; s++)
    kappa(s) = kappa(s-1) + mfexp(tmpk(s-1));

  dvariable P,eta;                      	// Variable allocation 

  for (int obs = 1; obs <= y(route).indexmax(); obs++)
  {
    eta = X(route)(obs)*beta + b_int + year(route)(obs)*b_year;
    if(y(route)(obs) == S)				
      P = 1.0;
    else
      P = 1/(1+exp(-(kappa(y(route)(obs))-eta)));
    if(y(route)(obs) > 1)
      P -= 1/(1+exp(-(kappa(y(route)(obs)-1)-eta)));
    g -= log(1.e-20+P);                   		// log[f(y|u)] for i'th super individuals
  }

NORMAL_PRIOR_FUNCTION void get_M(const dvariable& _a)

  int i,j;
  dvar_matrix tmpM(1,nroute,1,nroute);

  for (i=1;i<=nroute;i++)
  {
    tmpM(i,i)=1.0;
    for ( j=1;j<i;j++)
    {
      tmpM(i,j)=exp(-_a*dist(i,j));			// Exponentially decaying correlation
      tmpM(j,i)=tmpM(i,j);
    }
  }
  M=tmpM;

FUNCTION void evaluate_M(void)

  get_M(a);

REPORT_SECTION

  report << "#beta " << endl << beta << endl;
  report << "#kappa " << endl << kappa << endl;

TOP_OF_MAIN_SECTION
  arrmblsize = 4000000L;
  gradient_structure::set_MAX_NVAR_OFFSET(40207);
  gradient_structure::set_GRADSTACK_BUFFER_SIZE(3000000);
  gradient_structure::set_CMPDIF_BUFFER_SIZE(200000);