DATA_SECTION
 init_int n_sta
 init_int ncol_X //10
 init_ivector n_per_sta(1,n_sta)
 int nobs
  !! nobs = sum(n_per_sta);
 init_vector n(1,nobs)
 init_vector big(1,nobs)
 init_vector offst(1,nobs)
 
 init_matrix S(1,ncol_X,1,ncol_X) //zeros added for unpenalized intercept
 init_matrix X(1,nobs,1,ncol_X)
 init_int betas_phase
 init_int smooth_reff_phase
 init_int lambda_par_phase
 matrix U(1,ncol_X,1,ncol_X)
 vector u(1,ncol_X)
 int n_F
 int n_R
 LOCAL_CALCS
  U = eigenvectors(S);
  u = eigenvalues(S);
  n_F = 0;
  n_R = 0;
  for(int i = 1; i <= ncol_X; i++){
   if(u(i) < 1e-10) n_F +=1;
   else n_R += 1;
  }
  cout << "n_F: " << n_F << endl; //n_F = 2, number of fixed effects, intercept and 1 for smoother
  cout << "n_R: " << n_R << endl; //n_R = 8, number of random effects
 END_CALCS
 matrix U_F(1,ncol_X,1,n_F)
 matrix U_R(1,ncol_X,1,n_R)
 matrix bad_U_F(1,ncol_X,1,n_F)
 matrix X_F(1,nobs,1,n_F)
 matrix X_R(1,nobs,1,n_R)
 vector D_plus(1,n_R)
 int jj
 LOCAL_CALCS
  jj = 0;
  for(int i = 1; i <= ncol_X; i++){
   if(u(i) < 1e-10) {
    jj += 1;
    U_F.colfill(jj,column(U,i));
   }
   else {
    U_R.colfill(i-jj,column(U,i));
    D_plus(i-jj) = u(i); //eigenvalues of S are "known" variance components for random effects
   }
  }
  X_R = X * U_R;
  X_F = X * U_F;
  bad_U_F = U_F;
  bad_U_F.colfill(2,-1.0*column(U_F,2)); //I think second column corresponds to centering constraint 

 END_CALCS
  
PARAMETER_SECTION
 
 init_vector betas(1,n_F,betas_phase) //fixed effects (intercept and fixed effects portion of smoother)
 init_number lambda_par(lambda_par_phase)  //log smoothing parameter
 //init_number dummy(1)
 random_effects_vector smooth_reff(1,n_R,smooth_reff_phase) // random effects portion of smoother
 sdreport_vector transformed_betas(1,ncol_X) //betas for X that match mgcv results
 sdreport_vector wrong_transformed_betas(1,ncol_X) //gets se's for betas that match mgcv results
 objective_function_value nll

PROCEDURE_SECTION
	
	//nll = dummy*dummy;

	//prior distribution of smoother random effects
	nll -= -0.5*(n_R*(log(2.0*M_PI) - lambda_par) - sum(log(D_plus)) + mfexp(lambda_par)* 
		sum(elem_prod(smooth_reff, elem_prod(D_plus, smooth_reff))));
	
	dvar_vector eta(1,nobs);
	dvar_vector mu(1,nobs);
	
	eta = offst + X_F * betas + X_R * smooth_reff;
	mu = mfexp(-1.0*log(1 + mfexp(-eta)));
	nll -= sum(gammln(n + 1.0) - gammln(big + 1.0) - gammln(n - big + 1.0)) + big*log(mu) + (n-big)*log(1-mu);
	
	if(sd_phase()) {
		transformed_betas = U_F * betas + U_R * smooth_reff;
		wrong_transformed_betas = bad_U_F * betas + U_R * smooth_reff;
	}		
