#include <windows.h>
#include <TMB.hpp>

template<class Type>
Type objective_function<Type>::operator() ()
{
  DATA_INTEGER(n); // Number of obs
  DATA_VECTOR(y);  // Tree circumference 
  DATA_VECTOR(t);  // Days since 1968
  DATA_INTEGER(M); // Number of trees
  DATA_FACTOR(ngroup); // Number of measurements on each tree
  PARAMETER_VECTOR(beta);
  PARAMETER(log_sigma);
  PARAMETER(log_sigma_u);
  PARAMETER(log_sigma_v);
  PARAMETER_VECTOR(u);
  PARAMETER_VECTOR(v);

  Type sigma=exp(log_sigma);
  Type sigma_sq=exp(2.0*log_sigma);
  Type sigma_u=exp(log_sigma_u);
  Type sigma_u_sq=exp(2.0*log_sigma_u);
  Type sigma_v=exp(log_sigma_v);
  Type sigma_v_sq=exp(2.0*log_sigma_v); 
  
  vector<Type> mu(n);
  vector<Type> a(3); // Parameters of curve, offset by constants
  a[0] = 200.0 + beta[0];
  a[1] = 700.0 + beta[1];
  a[2] = 300.0 + beta[2];

  int i,j,k;
  Type g=0;
  
  k=0;
  for(i=0;i< 7;i++) 
  { g -= dnorm(v[i],Type(0.0),sigma_v,TRUE); } // Random day effects contribution
  for(i=0;i< M;i++) 
  {
    g -= dnorm(u[i],Type(0.0),sigma_u,TRUE);  // Random tree effects contribution
    for(j=0;j<ngroup(i);j++)
    {
      mu[k] = (a[0]+u[i]+v[j])/(1+exp(-(t[k]-a[1])/a[2]));
      g -= dnorm(y[k],mu[k],sigma,TRUE); // Observation terms
      k++;
    }
  }
 
  freopen("out.txt","a",stdout);
  std::cout << "beta=" << beta << std::endl;
  std::cout << "log_sigma=" << log_sigma << std::endl;
  std::cout << "log_sigma_u=" << log_sigma_u << std::endl;  
  std::cout << "sigma_u=" << sigma_u << std::endl;  
  fclose(stdout);
  
  // These need to be here, not at top of calculations
  REPORT(sigma);
  REPORT(sigma_u);
  REPORT(sigma_v);
  
  ADREPORT(sigma);
  ADREPORT(sigma_sq);
  ADREPORT(sigma_u);
  ADREPORT(sigma_u_sq);  
  ADREPORT(sigma_v);
  ADREPORT(sigma_v_sq); 
  
  return g;
}