// Space time 
#include <TMB.hpp>


template<class Type>
Type objective_function<Type>::operator() ()
{
  // Data
  DATA_INTEGER(n_data);         // Total number of observations
  DATA_VECTOR(Y);       	// Count data
  DATA_FACTOR(NAind);		// 1 = Y is NA, 0 = is not NA
  DATA_FACTOR(NegBin);       // Use negative binomial? 0:No; 1:Yes
  DATA_INTEGER(n_stations)	// Number of stations 
  DATA_INTEGER(n_years)          // Number of years 
  DATA_INTEGER(n_p)          	// number of columns in covariate matrix X
  DATA_MATRIX(X);		// Covariate design matrix

  // Fixed effects
  PARAMETER_VECTOR(alpha);  // Gompertz-drift parameter
  PARAMETER(phi);           // Offset of beginning from equilibrium
  PARAMETER(log_tau);        // log-inverse SD of process variability
  PARAMETER(rho);                // Autocorrelation (i.e. density dependence)
  PARAMETER_VECTOR(ln_VarInfl);   // Overdispersion parameters
  
  // Random effects
  PARAMETER_VECTOR(Epsilon_input);      // Process variation

  using namespace density;
  int i,j;
  Type g = 0;
  
  // Probability of random effects
  g = AR1(rho)(Epsilon_input);

  // Likelihood contribution from observations
  vector<Type> eta(n_data); 
  vector<Type> nu(n_data);
  vector<Type> mean_abundance(n_years);
  matrix<Type> Dji(n_stations,n_years);
  vector<Type> Epsilon(n_years);
  eta = X*alpha.matrix();
  
  // Probability of data
  int ii = 0;
  for (int i=0;i<n_years;i++){
    for (int j=0;j<n_stations;j++){ 
      Epsilon(i) = Epsilon_input(i)/exp(log_tau);
      nu[ii] = phi * pow(rho, i);
      Dji(j,i) =  nu[ii] + ( eta(ii) + Epsilon_input(i)/exp(log_tau) )/(1-rho);
      mean_abundance(i) = mean_abundance(i) + exp( Dji(j,i) );      
      if(!NAind(ii)){               // nu(ii) + 
        Type mean_y = exp( Dji(j,i) );
        Type var_y = mean_y*(1.0+exp(ln_VarInfl(0)))+pow(mean_y,2.0)*exp(ln_VarInfl(1));
        if(!NegBin(0)){ g -= dpois( Y[ii], mean_y, 1 ); }
        if(NegBin(0)){ g -= dnbinom2( Y[ii], mean_y, var_y, 1 ); }
      }
      ii++;
    }
    mean_abundance[i] = mean_abundance[i] / n_stations;      
  }
  
  // Change in abundance
  Type Change = mean_abundance[n_years-1] / mean_abundance[0]; 
  Type Prop_Change = (mean_abundance[n_years-1] - mean_abundance[0]) / mean_abundance[0]; 
  REPORT( Change );
  REPORT( log(Change) );
  REPORT( Prop_Change );
  ADREPORT( Change );
  ADREPORT( log(Change) );
  ADREPORT( Prop_Change );
  // Parameters
  REPORT( Dji );
  REPORT( Epsilon );
  // Total abundance
  ADREPORT( log(mean_abundance) );
  ADREPORT(mean_abundance);
  return g;
}