diff --git a/.gitattributes b/.gitattributes
new file mode 100644
index 0000000..e69de29
diff --git a/.gitignore b/.gitignore
index a75452f..f5d15ff 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,4 @@
 .DS_Store
 .*~
 __MACOSX
+#*
diff --git a/bsglmm/covarGPU.cu b/bsglmm/covarGPU.cu
index 8c46122..9bc4a44 100644
--- a/bsglmm/covarGPU.cu
+++ b/bsglmm/covarGPU.cu
@@ -213,7 +213,7 @@ __global__ void Spat_Coef_for_Spat_Covar(float *covar,float *alpha,float *Z,floa
 	float var[20*20];
 	float tmp=0;
 	float tmp2[20];
-//	float alph[20];
+	float alph[20];
 	
 	int voxel=0,IDX=0,IDXSC=0;
 	int current_covar=0;
@@ -302,8 +302,8 @@ __global__ void Spat_Coef_Probit(float *covar,float *covF,float *fix,float *alph
 	float var[20*20];
 	float tmp=0;
 	float tmp2[20];
-//	float alph[20];
-//	float fixed[20];
+	float alph[20];
+	float fixed[20];
 	
 	int voxel=0,IDX=0,IDXSC=0;
 	float bwhtmat=0;
@@ -407,7 +407,7 @@ __global__ void Spat_Coef(float *covar,float *covF,float *fix,float *alpha,float
 	float var[20*20];
 	float tmp=0;
 	float tmp2[20];
-//	float alph[20];
+	float alph[20];
 	float fixed[20];
 	
 	int voxel=0,IDX=0,IDXSC=0;
@@ -1018,7 +1018,7 @@ void updateSpatCoefPrecGPU_Laplace(float *SpatCoefPrec,unsigned long *seed)
 	double *var;
 	float *dbeta;
 	double betaSqr = 2.0;//0.2;
-//	double tau;
+	double tau;
 
 	var = (double *)calloc(NCOVAR*NCOVAR,sizeof(double));
 	hbeta = (float*)calloc(BPG,sizeof(float));
@@ -1372,7 +1372,7 @@ __global__ void Z_GPU1(float *dZ,unsigned char *data,float *covar,float *covarF,
 	int localsize = blockDim.x;  
 	int localid = threadIdx.x; 
 	
-	float mean[300];
+	float mean[2000];
 	float tmp=0;
 	float SC[20];
 	int voxel=0,IDX=0,IDXSC=0;
@@ -1421,7 +1421,7 @@ __global__ void Z_GPU1(float *dZ,unsigned char *data,float *covar,float *covarF,
 			int itmp = (int)data[isub*TOTVOX+IDX];
 			dZ[isub*TOTVOX+IDX] = truncNormGPU(mean[isub],1.0f,lim[itmp].x,lim[itmp].y,&localState);
 //			dZ[isub*TOTVOX+IDX] = ((float)data[isub*TOTVOX+IDX]-1.0f)*truncNormGPU(fabsf(mean[isub]),1.0f,0.0f,500.0f,&localState);
-			dR[isub*TOTVOX+IDX] += (fabsf(dZ[isub*TOTVOX+IDX] - mean[isub]) > 4.0f);
+//			dR[isub*TOTVOX+IDX] += (fabsf(dZ[isub*TOTVOX+IDX] - mean[isub]) > 4.0f);
 		}
 		state[idx] = localState;
 	}
@@ -1435,7 +1435,7 @@ __global__ void Z_GPU2(float *dZ,unsigned char *data,float *covar,float *covarF,
 	int localid = threadIdx.x; 
 	
 	float a = alp;
-	float mean[300];
+	float mean[2000];
 	float tmp=0;
 	float SC[20];
 	int voxel=0,IDX=0,IDXSC=0;
@@ -1489,7 +1489,7 @@ __global__ void Z_GPU2(float *dZ,unsigned char *data,float *covar,float *covarF,
 			int itmp = (int)data[isub*TOTVOX+IDX];
 			Z = truncNormGPU(M,P,lim[itmp].x,lim[itmp].y,&localState);
 //			Z = ((float)data[isub*TOTVOX+IDX]-1.0f)*truncNormGPU(fabsf(M),rsqrtf(P),0.0f,500.0f,&localState);
-			dR[isub*TOTVOX+IDX] += (fabsf(Z - M) > 4.0f);
+//			dR[isub*TOTVOX+IDX] += (fabsf(Z - M) > 4.0f);
 			
 			b = Z - M;
 			b = (df + b*b)/2.0f;
@@ -1511,7 +1511,7 @@ __global__ void Phi_GPU(float *dZ,unsigned char *data,float *covar,float *covarF
 	int localid = threadIdx.x; 
 	
 	const float a = alp;
-	float beta2[300];
+	float beta2[2000];
 	float tmp=0;
 	float SC[20];
 	int voxel=0,IDX=0,IDXSC=0;
diff --git a/bsglmm/main.cu b/bsglmm/main.cu
index 1d65db8..285a5bf 100644
--- a/bsglmm/main.cu
+++ b/bsglmm/main.cu
@@ -27,8 +27,8 @@ float *covar;
 int NSUBTYPES;
 int NCOVAR;
 int NCOV_FIX = 0;
-int MAXITER = 1000000;
-int BURNIN  = 500000;
+int MAXITER = 100000;
+int BURNIN  =  50000;
 int RESTART;
 float *deviceCovar;
 float *hostCovar;
@@ -69,12 +69,31 @@ int main (int argc, char * const argv[]) {
 	void write_empir_prb(unsigned char *,float *,unsigned char *,int *);
 	unsigned char *get_WM_mask(float *,unsigned char *);
 	
-	if (argc !=5) 
-		exit(0);
+	if (argc !=5 && argc !=7) {
+	   printf("%s: Usage\n",argv[0]);
+	   printf("%s  NTypes  NCov  GPU  Design  [MaxIter BurnIn]  \n",argv[0]);
+	   printf("  NTypes  - Number of groups\n",argv[0]);
+	   printf("  NCov    - Number of covariates (count must include groups)\n",argv[0]);
+	   printf("  GPU     - 1 use GPU; 0 use CPU\n",argv[0]);
+	   printf("  Design  - Text file, tab or space separated data file\n",argv[0]);
+	   printf("  MaxIter - Number of iterations (defaults to 1,000,000)\n",argv[0]);
+	   printf("  BurnIn  - Number of burn-in iterations (defaults to 500,000)\n",argv[0]);
+	   printf("For documentation see: http://warwick.ac.uk/tenichols/BSGLMM\n",argv[0]);
+		exit(1);
+		}
 
 	NSUBTYPES = atoi(argv[1]);
 	NCOVAR = atoi(argv[2]);
 	GPU = atoi(argv[3]);
+	if (argc >5)  {
+		MAXITER = atoi(argv[5]);
+		BURNIN  = atoi(argv[6]);
+		if (MAXITER<1000) 
+		   printf("WARNING: Silly-small number of iterations used; recommend abort and use more\n");
+		if (BURNIN>MAXITER)
+		   printf("WARNING: Burn-in exceeds MAXITER, no results will be saved\n");
+	}	
+
 	
 	RESTART = 0;
 	int deviceCount = 0;
@@ -113,6 +132,18 @@ int main (int argc, char * const argv[]) {
 	
 	seed = (unsigned long *)calloc(3,sizeof(unsigned long));
 	fseed = fopen("seed.dat","r");
+	if (fseed == NULL){
+		printf("No seed.dat found; using [ 1 2 3 ]\n");
+		seed[0]=1L;seed[1]=2L;seed[2]=2L;
+	} else {
+		rtn = fscanf(fseed,"%lu %lu %lu\n",&(seed[0]),&(seed[1]),&(seed[2]));
+		if (rtn==0)
+		  exit(1);
+		rtn = fclose(fseed);
+		if (rtn != 0)
+		  exit(1);
+	}
+
 	
 	logit_factor = 1.0f;
 	t_df = 1000.0f;
@@ -124,13 +155,6 @@ int main (int argc, char * const argv[]) {
 		t_df = 3.0f;
 
 	
-	rtn = fscanf(fseed,"%lu %lu %lu\n",&(seed[0]),&(seed[1]),&(seed[2]));
-	if (rtn==0)
-		exit(0);
-	rtn = fclose(fseed);
-	if (rtn != 0)
-		exit(0);
-
 		
 //	msk  = read_nifti1_mask("./images/avg152T1_highres_brain.hdr",
 //				"./images/avg152T1_highres_brain.img");
@@ -250,6 +274,9 @@ int main (int argc, char * const argv[]) {
 
 	}
 	else {
+	     fprintf(stderr,"WARNING!!!\n");
+	     fprintf(stderr,"WARNING!!! CPU code not tested!  Results might not be right!\n");
+	     fprintf(stderr,"WARNING!!!\n");
 		for (int i=0;i<2;i++) {
 			INDX[i].hostVox = (int *)calloc(INDX[i].hostN,sizeof(int));
 		}
diff --git a/bsglmm/mcmc.cpp b/bsglmm/mcmc.cpp
index 553d2f0..bf8e614 100644
--- a/bsglmm/mcmc.cpp
+++ b/bsglmm/mcmc.cpp
@@ -87,10 +87,10 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 	float Prec = 0.000001f;
 	SUB  *subj;
 	char *S;
-	
+
 //	float *hostWM;
 	FILE *out,*out2,*fWM,*fcoef,*fcoefsd,*fBF,*fresid;
-	
+
 	void itoa(int,char *);
 //	unsigned char *get_neighbors();
 	void initializeAlpha(unsigned char *,float *,float *,float *,unsigned long *);
@@ -99,7 +99,7 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 	void updateZGPU(unsigned char *,float,unsigned long *);
 	void updateZ(float *,float *,unsigned char *,unsigned char *,float, float *,float *,float *,unsigned long *);
 	void updatePhiGPU(float,float *);
-	void updatePhi(unsigned char *data,float *Z,float *Phi,unsigned char *msk,float beta,float *WM,float *spatCoef,float *covar,unsigned long *seed); 
+	void updatePhi(unsigned char *data,float *Z,float *Phi,unsigned char *msk,float beta,float *WM,float *spatCoef,float *covar,unsigned long *seed);
 	void updateAlpha(float *,float ,float ,float *,unsigned char *,float,float *,float *,float *,int,unsigned long *);
 	void updateAlphaGPU(float *,float *,float *,float *,unsigned char *,float,float *,float *,float *,unsigned long *);
 	void updateAlphaPrecGPU(float *,float *,float,float,float,unsigned long *);
@@ -124,9 +124,8 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 	double compute_prb_DIC(double *,double *,float *,unsigned char *,unsigned char *,int);
 
 //	void tmp(float *,float *);
+	fresid = fopen("resid.dat","w");
 
-//	fresid = fopen("resid.dat","w");
-	
  	NSIZE = (NROW+2)*(NCOL+2)*(NDEPTH+2);
 	int RSIZE = NROW*NCOL*NDEPTH;
 //	printf("%d %d \n",TOTVOX,NSIZE);
@@ -139,11 +138,11 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 	//		if (data[i] == 1)
 	//			data2[i] = 2;
 	//	}
-	//	CUDA_CALL( cudaMemcpy(deviceData,data2,NSUBS*TOTVOX*sizeof(unsigned char),cudaMemcpyHostToDevice) );				
-		CUDA_CALL( cudaMemcpy(deviceData,data,NSUBS*TOTVOX*sizeof(unsigned char),cudaMemcpyHostToDevice) );				
+	//	CUDA_CALL( cudaMemcpy(deviceData,data2,NSUBS*TOTVOX*sizeof(unsigned char),cudaMemcpyHostToDevice) );
+		CUDA_CALL( cudaMemcpy(deviceData,data,NSUBS*TOTVOX*sizeof(unsigned char),cudaMemcpyHostToDevice) );
 	//	free(data2);
 	}
-	
+
 	SpatCoefMean = (float *)calloc(NCOVAR,sizeof(float));
 	SpatCoefPrec = (float *)calloc(NCOVAR*NCOVAR,sizeof(float));
 	for (int i=0;i<NCOVAR;i++) {
@@ -158,16 +157,16 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 	alphaPrec = (float *)calloc(NSUBTYPES,sizeof(float));
 	for (int i = 0;i<NSUBTYPES;i++)
 		alphaPrec[i] = 1.0;
-	
+
 	if (GPU) {
 	 	CUDA_CALL( cudaMalloc((void **)&deviceWM,TOTVOX*sizeof(float)) );
 	 	CUDA_CALL( cudaMemcpy(deviceWM,WM,TOTVOX*sizeof(float),cudaMemcpyHostToDevice) );
 		CUDA_CALL( cudaDeviceSynchronize() );
-	}	
+	}
+
 
-		
 	subj = (SUB *)calloc(NSUBS,sizeof(SUB));
-	
+
 	Z = (float *)calloc(TOTVOX*NSUBS,sizeof(float));
 	Phi = (float *)calloc(TOTVOX*NSUBS,sizeof(float));
 	if (GPU) {
@@ -175,22 +174,22 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 		CUDA_CALL( cudaMalloc((void **)&deviceZ,TOTVOX*NSUBS*sizeof(float)) );
 		CUDA_CALL( cudaMemset(deviceZ,0,TOTVOX*NSUBS*sizeof(float)) );
 	//	CUDA_CALL( cudaHostAlloc((void **)&Phi, TOTVOX*NSUBS*sizeof(float),cudaHostAllocDefault) );
-		CUDA_CALL( cudaMalloc((void **)&devicePhi,TOTVOX*NSUBS*sizeof(float)) );
-		CUDA_CALL( cudaMemset(devicePhi,1,TOTVOX*NSUBS*sizeof(float)) );
-		CUDA_CALL( cudaMalloc((void **)&deviceResid,TOTVOX*NSUBS*sizeof(unsigned int)) );
-		CUDA_CALL( cudaMemset(deviceResid,0,TOTVOX*NSUBS*sizeof(unsigned int)) );
+//		CUDA_CALL( cudaMalloc((void **)&devicePhi,TOTVOX*NSUBS*sizeof(float)) );
+//		CUDA_CALL( cudaMemset(devicePhi,1,TOTVOX*NSUBS*sizeof(float)) );
+//		CUDA_CALL( cudaMalloc((void **)&deviceResid,TOTVOX*NSUBS*sizeof(unsigned int)) );
+//		CUDA_CALL( cudaMemset(deviceResid,0,TOTVOX*NSUBS*sizeof(unsigned int)) );
 	}
 //	else
 //		Z = (float *)calloc(TOTVOX*NSUBS,sizeof(float));
 
-	if (!RESTART) {
+/*	if (!RESTART) {
 		initializeZ(Z,data,msk,seed);
-		initializePhi(Phi,data,msk,t_df,seed);
+//		initializePhi(Phi,data,msk,t_df,seed);
 		if (GPU) {
-			CUDA_CALL( cudaMemcpy(deviceZ,Z,NSUBS*TOTVOX*sizeof(float),cudaMemcpyHostToDevice) );	
-			CUDA_CALL( cudaMemcpy(devicePhi,Phi,NSUBS*TOTVOX*sizeof(float),cudaMemcpyHostToDevice) );	
+			CUDA_CALL( cudaMemcpy(deviceZ,Z,NSUBS*TOTVOX*sizeof(float),cudaMemcpyHostToDevice) );
+//			CUDA_CALL( cudaMemcpy(devicePhi,Phi,NSUBS*TOTVOX*sizeof(float),cudaMemcpyHostToDevice) );
 		}
-	}	
+	}	*/
 	alpha = (float *)calloc(TOTVOX*NSUBTYPES,sizeof(float));
 	if (GPU) {
 //		cutilSafeCall(cudaMemcpy(deviceZ,Z,NSUBS*TOTVOX*sizeof(float),cudaMemcpyHostToDevice));
@@ -211,13 +210,13 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 		if (GPU)
 		       cudaMemcpy(deviceAlpha,alpha,NSUBTYPES*TOTVOX*sizeof(float),cudaMemcpyHostToDevice);
 	}*/
-	
+
 //	Malpha = (double *)calloc(NSUBTYPES*RSIZE,sizeof(double));
 	MWM = (double *)calloc(RSIZE,sizeof(double));
 	Mprb = (double *)calloc(NSUBTYPES*RSIZE,sizeof(double));
 	MCov = (double *)calloc(NCOVAR*RSIZE,sizeof(double));
 	SCov = (double *)calloc(NCOVAR*RSIZE,sizeof(double));
-	
+
 	Qhat = (double **)calloc(NSUBS,sizeof(double *));
 	for (int i=0;i<NSUBS;i++)
 		Qhat[i] = (double *)calloc(NSUBTYPES,sizeof(double));
@@ -227,7 +226,7 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 		for (int j=1;j<NCOL+1;j++) {
 			for (int i=1;i<NROW+1;i++) {
 				IDX = idx(i,j,k);
-				if (msk[IDX]) 
+				if (msk[IDX])
 					spatCoef[hostIdxSC[IDX]] = -25.0f;
 			}
 		}
@@ -252,8 +251,8 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 //	else {
 //		prb = (float *)calloc(NSUBTYPES*TOTVOX,sizeof(float));
 //	}
-	
-	if (RESTART) {
+
+/*	if (RESTART) {
 		restart_iter(&beta,fixMean,alphaMean,alphaPrec,SpatCoefMean,SpatCoefPrec,spatCoef,alpha,Z,Phi);
 		printf("***RESTART VALUES***\n\n");
 		double *V = (double *)calloc(NCOVAR*NCOVAR,sizeof(double));
@@ -264,83 +263,77 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 		cholesky_invert(NCOVAR, V);
 		printf("iter = %d\t",iter);
 		printf("%f \n",beta);
-//		for (int ii=0;ii<NCOV_FIX;ii++)
-//			printf("%f \n",fixMean[ii]);
-//		for (int ii=0;ii<NSUBTYPES;ii++)
-//			printf("\t %10.6lf %10.6f\n",alphaMean[ii],V[ii+ii*NCOVAR]);
-//		printf("\n");
-//		for (int ii=NSUBTYPES;ii<NCOVAR;ii++)
-//			printf("\t %10.6lf %10.6f\n",alphaMean[ii],V[ii+ii*NCOVAR]);
-//		printf("\n");
-//		for (int ii=0;ii<NCOVAR;ii++) {
-//			for (int jj=0;jj<=ii;jj++)
-//				printf("%6.3lf ",V[jj + ii*NCOVAR]/sqrt(V[ii + ii*NCOVAR]*V[jj + jj*NCOVAR]));
-//			printf("\n");
-//		}
+		for (int ii=0;ii<NCOV_FIX;ii++)
+			printf("%f \n",fixMean[ii]);
+		for (int ii=0;ii<NSUBTYPES;ii++)
+			printf("\t %10.6lf %10.6f\n",alphaMean[ii],V[ii+ii*NCOVAR]);
+		printf("\n");
+		for (int ii=NSUBTYPES;ii<NCOVAR;ii++)
+			printf("\t %10.6lf %10.6f\n",alphaMean[ii],V[ii+ii*NCOVAR]);
+		printf("\n");
+		for (int ii=0;ii<NCOVAR;ii++) {
+			for (int jj=0;jj<=ii;jj++)
+				printf("%6.3lf ",V[jj + ii*NCOVAR]/sqrt(V[ii + ii*NCOVAR]*V[jj + jj*NCOVAR]));
+			printf("\n");
+		}
 		printf("\n");
 		fflush(NULL);
 		free(V);
 
 
-/*		printf("%f \n",beta);
-		for (int ii=0;ii<NSUBTYPES;ii++)
-			printf("\t %10.6lf %10.6lf\t %10.6lf %10.6lf\n",alphaMean[ii],1./alphaPrec[ii],SpatCoefMean[ii],1./SpatCoefPrec[ii]);
-		for (int ii=NSUBTYPES;ii<NCOVAR;ii++)
-			printf("\t \t \t \t %10.6lf %10.6lf\n",alphaMean[ii],1./SpatCoefPrec[ii]);
-		printf("\n");
-		fflush(NULL);*/
-//		out  = fopen("parms.dat","a");
-//		out2  = fopen("select_parms.dat","w");
-	}
-	else {
-//		out = fopen("parms.dat","w");	
-//		out2 = fopen("select_parms.dat","w");	
-	}
-//	fBF = fopen("fBF.dat","w");
+		out  = fopen("parms.dat","a");
+		out2  = fopen("select_parms.dat","w");
+	}*/
+	//else {
+		out = fopen("parms.dat","w");
+		out2 = fopen("select_parms.dat","w");
+	//}
+	fBF = fopen("fBF.dat","w");
 
 	cudaEvent_t start, stop;
 	float time;
 	CUDA_CALL( cudaEventCreate(&start) );
 	CUDA_CALL( cudaEventCreate(&stop) );
 //cudaEventRecord(start,0);
-
+    
 	if ((MAXITER-BURNIN)<10000){
-		SAVE_ITER = MAXITER-BURNIN;
+	  SAVE_ITER = MAXITER-BURNIN;
 	}
 	else {
-		SAVE_ITER = (MAXITER - BURNIN)/10000;
+	  SAVE_ITER = (MAXITER - BURNIN)/10000;
 	}
 	PRNtITER = 100;
 
 	for (iter=0;iter<=MAXITER;iter++) {//printf("iter = %d\n",iter);fflush(NULL);
 
-		
+
 		if (GPU) {
 			if (!(iter%PRNtITER))
-				CUDA_CALL( cudaEventRecord(start, 0) ); 
+				CUDA_CALL( cudaEventRecord(start, 0) );
 			//for (int tj=0;tj<10;tj++)
 			updateZGPU(data,beta,seed);
+
 			if (!(iter%PRNtITER)) {
 				CUDA_CALL( cudaEventRecord(stop, 0) );
 				CUDA_CALL( cudaEventSynchronize(stop) );
 				CUDA_CALL( cudaEventElapsedTime(&time, start, stop) );
-//				printf ("Time to updateZGPU kernel: %f ms\n", time);
+				printf ("Time to updateZGPU kernel: %f ms\n", time);
 			}
-			
+
 			if (MODEL != 1) {
 				if (!(iter%PRNtITER))
-					CUDA_CALL( cudaEventRecord(start, 0) ); 
+					CUDA_CALL( cudaEventRecord(start, 0) );
 				//for (int tj=0;tj<10;tj++)
 				updatePhiGPU(beta,Phi);
-				//updatePhi(data,Z,Phi,msk,beta,WM,spatCoef,covar,seed);			
+				//updatePhi(data,Z,Phi,msk,beta,WM,spatCoef,covar,seed);
 				if (!(iter%PRNtITER)) {
 					CUDA_CALL( cudaEventRecord(stop, 0) );
 					CUDA_CALL( cudaEventSynchronize(stop) );
 					CUDA_CALL( cudaEventElapsedTime(&time, start, stop) );
-//					printf ("Time to updatePhiGPU kernel: %f ms\n", time);
+					printf ("Time to updatePhiGPU kernel: %f ms\n", time);
 				}
 			}
-			
+
 			if (!(iter%PRNtITER))
 				CUDA_CALL( cudaEventRecord(start, 0) );
 			//for (int tj=0;tj<10;tj++)
@@ -349,7 +342,7 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 				CUDA_CALL( cudaEventRecord(stop, 0) );
 				CUDA_CALL( cudaEventSynchronize(stop) );
 				CUDA_CALL( cudaEventElapsedTime(&time, start, stop) );
-//				printf ("Time for the updateBetaGPU kernel: %f ms\n", time);
+				printf ("Time for the updateBetaGPU kernel: %f ms\n", time);
 			}
 
 			if (!(iter%PRNtITER))
@@ -360,7 +353,7 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 				CUDA_CALL( cudaEventRecord(stop, 0) );
 				CUDA_CALL( cudaEventSynchronize(stop) );
 				CUDA_CALL( cudaEventElapsedTime(&time, start, stop) );
-//				printf ("Time to updateSpatCoefGPU kernel: %f ms\n", time);
+				printf ("Time to updateSpatCoefGPU kernel: %f ms\n", time);
 			}
 
 			if (!(iter%PRNtITER))
@@ -372,7 +365,7 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 				CUDA_CALL( cudaEventRecord(stop, 0) );
 				CUDA_CALL( cudaEventSynchronize(stop) );
 				CUDA_CALL( cudaEventElapsedTime(&time, start, stop) );
-//				printf ("Time to updateSpatCoefPrecGPU kernel: %f ms\n\n", time);
+				printf ("Time to updateSpatCoefPrecGPU kernel: %f ms\n\n", time);
 			}
 			//CUDA_CALL( cudaEventRecord(start, 0) );
 			//for (int tj=0;tj<10;tj++)
@@ -394,19 +387,19 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 				updatefixMeanGPU(fixMean,Prec,deviceCovar,deviceZ,deviceSpatCoef,beta,seed);
 		//	updateAlphaMeanGPU(alphaMean,Prec,deviceCovar,deviceZ,deviceSpatCoef,beta,seed);
 		}
-		else {	
+		else {
 			if (!(iter%PRNtITER))
 				CUDA_CALL( cudaEventRecord(start, 0) );
-	
+
 		//	for (int tj=0;tj<10;tj++)
 			updateZ(Z,alpha,data,msk,beta,WM,spatCoef,covar,seed);
-			if (!(iter%PRNtITER)) {			
+			if (!(iter%PRNtITER)) {
 				CUDA_CALL( cudaEventRecord(stop, 0) );
 				CUDA_CALL( cudaEventSynchronize(stop) );
 				CUDA_CALL( cudaEventElapsedTime(&time, start, stop) );
-//				printf ("Time for the updateZCPU kernel: %f ms\n", time);
+				printf ("Time for the updateZCPU kernel: %f ms\n", time);
 			}
-		//	CUDA_CALL( cudaEventRecord(start, 0) );	
+		//	CUDA_CALL( cudaEventRecord(start, 0) );
 		//	for (int tj=0;tj<10;tj++)
 			updateBeta(&beta,Z,alpha,WM,msk,prior_mean_beta,prior_prec_beta,spatCoef,covar,seed);
 		//	CUDA_CALL( cudaEventRecord(stop, 0) );
@@ -414,27 +407,27 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 		//	CUDA_CALL( cudaEventElapsedTime(&time, start, stop) );
 		//	printf ("Time for the updateBetaCPU kernel: %f ms\n", time);
 
-			if (!(iter%PRNtITER))		
-				CUDA_CALL( cudaEventRecord(start, 0) );	
+			if (!(iter%PRNtITER))
+				CUDA_CALL( cudaEventRecord(start, 0) );
 		//	for (int tj=0;tj<10;tj++)
 			updateSpatCoef(covar,spatCoef,SpatCoefPrec,alpha,alphaMean,Z,msk,beta,WM,seed);
-			if (!(iter%PRNtITER)) {				
+			if (!(iter%PRNtITER)) {
 				CUDA_CALL( cudaEventRecord(stop, 0) );
 				CUDA_CALL( cudaEventSynchronize(stop) );
 				CUDA_CALL( cudaEventElapsedTime(&time, start, stop) );
-//				printf ("Time for the updateSpatCoefCPU kernel: %f ms\n", time);
+				printf ("Time for the updateSpatCoefCPU kernel: %f ms\n", time);
 			}
 			if (!(iter%PRNtITER))
 				CUDA_CALL( cudaEventRecord(start, 0) );
 		//	for (int tj=0;tj<10;tj++)
 			updateSpatCoefPrec(SpatCoefPrec,spatCoef,msk,seed);
-			if (!(iter%PRNtITER)) {		
+			if (!(iter%PRNtITER)) {
 				CUDA_CALL( cudaEventRecord(stop, 0) );
 				CUDA_CALL( cudaEventSynchronize(stop) );
 				CUDA_CALL( cudaEventElapsedTime(&time, start, stop) );
-//				printf ("Time for the updateSpatCoefPrecCPU kernel: %f ms\n\n", time);
+				printf ("Time for the updateSpatCoefPrecCPU kernel: %f ms\n\n", time);
 			}
-	
+
 		//	CUDA_CALL( cudaEventRecord(start, 0) );
 		//	for (int tj=0;tj<10;tj++)
 		//		compute_prb(prb,alpha,spatCoef,beta,WM,msk);
@@ -450,7 +443,7 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 			int j = 77;//46;
 			int k = 19;//72;
 			IDX = idx(i,j,k);
-	
+
 			for (int ii=NSUBTYPES;ii<NCOVAR;ii++)
 				fprintf(out2,"%f ",spatCoef[ii*TOTVOXp+hostIdxSC[IDX]]/logit_factor);
 
@@ -458,7 +451,7 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 			j = 91;//46;
 			k = 49;//72;
 			IDX = idx(i,j,k);
-	
+
 			for (int ii=NSUBTYPES;ii<NCOVAR;ii++)
 				fprintf(out2,"%f ",spatCoef[ii*TOTVOXp+hostIdxSC[IDX]]/logit_factor);
 
@@ -466,15 +459,15 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 			j = 39;//46;
 			k =  5;//72;
 			IDX = idx(i,j,k);
-	
+
 			for (int ii=NSUBTYPES;ii<NCOVAR;ii++)
 				fprintf(out2,"%f ",spatCoef[ii*TOTVOXp+hostIdxSC[IDX]]/logit_factor);
-	
+
 		 	i = 17;//57;
 			j = 42;//46;
 			k = 12;//72;
 			IDX = idx(i,j,k);
-	
+
 			for (int ii=NSUBTYPES;ii<NCOVAR;ii++)
 				fprintf(out2,"%f ",spatCoef[ii*TOTVOXp+hostIdxSC[IDX]]/logit_factor);
 
@@ -482,7 +475,7 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 			j = 93;//46;
 			k = 30;//72;
 			IDX = idx(i,j,k);
-	
+
 			for (int ii=NSUBTYPES;ii<NCOVAR;ii++)
 				fprintf(out2,"%f ",spatCoef[ii*TOTVOXp+hostIdxSC[IDX]]/logit_factor);
 			fprintf(out2,"\n");fflush(NULL);
@@ -492,7 +485,7 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 //			ss[1] = 134;
 //			ss[2] = 141;
 //			ss[3] = 39;
-			
+
 //			float *zs;
 //			zs = (float *)malloc(sizeof(float));
 //			for (int i=0;i<4;i++) {
@@ -502,14 +495,14 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 //					fprintf(out2,"%f ",spatCoef[ii*TOTVOXp+hostIdxSC[IDX]]/logit_factor);
 //				}
 //				cudaMemcpy(zs,&(deviceZ[ss[i]*TOTVOX+hostIdx[IDX]]),sizeof(float),cudaMemcpyDeviceToHost);
-//			
+//
 			//	printf("%7.4f ",zs[0] - xb);
 //				fprintf(fresid,"%f ",zs[0] - xb);
 //			}
 //			fprintf(fresid,"\n");
 //			free(zs);
 //			free(ss);
-			
+
 	//		i = 59;
 	//		j = 25;
 	//		k = 48;
@@ -517,11 +510,11 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 //			j = 34;
 //			k = 7;
 //			IDX = idx(i,j,k);
-			
+
 //			for (int ii=0;ii<NCOVAR;ii++)
 //				fprintf(out2,"%f ",spatCoef[ii*TOTVOXp+hostIdxSC[IDX]]/logit_factor);
 //			fprintf(out2,"\n");
-		
+
 			fprintf(out,"%f ",beta);
 			for (int ii=0;ii<NCOV_FIX;ii++)
 				fprintf(out,"%f ",fixMean[ii]);
@@ -529,11 +522,11 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 				fprintf(out,"%f ",alphaMean[ii]/logit_factor);
 			for (int ii=0;ii<NCOVAR;ii++)
 				for (int jj=0;jj<=ii;jj++)
-					fprintf(out,"%f ",SpatCoefPrec[jj + ii*NCOVAR]);				
+					fprintf(out,"%f ",SpatCoefPrec[jj + ii*NCOVAR]);
 			fprintf(out,"\n");
 		fflush(NULL);
 		}*/
-		
+
 		if (!(iter%PRNtITER)) {
 /*			CUDA_CALL( cudaMemcpy(spatCoef,deviceSpatCoef,NCOVAR*TOTVOXp*sizeof(float),cudaMemcpyDeviceToHost) );
 			float max = -500;
@@ -541,14 +534,14 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 				for (int j=1;j<NCOL+1;j++) {
 					for (int i=1;i<NROW+1;i++) {
 						IDX = idx(i,j,k);
-						if (msk[IDX]) 
+						if (msk[IDX])
 							max = (max >	(double)(spatCoef[4*TOTVOXp+hostIdxSC[IDX]]/logit_factor)) ?
 									max:(double)(spatCoef[4*TOTVOXp+hostIdxSC[IDX]]/logit_factor);
 					}
 				}
 			}
-			printf("max = %lf\n",max);	*/			
-			
+			printf("max = %lf\n",max);	*/
+
 			double *V = (double *)calloc(NCOVAR*NCOVAR,sizeof(double));
 			for (int ii=0;ii<NCOVAR;ii++)
 				for (int jj=0;jj<NCOVAR;jj++)
@@ -556,21 +549,21 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 			cholesky_decomp(V, NCOVAR);
 			cholesky_invert(NCOVAR, V);
 			printf("iter = %d\t",iter);
-//			printf("%f \n",beta);
-//			for (int ii=0;ii<NCOV_FIX;ii++)
-//				printf("%f \n",fixMean[ii]);
-//			for (int ii=0;ii<NSUBTYPES;ii++)
-//				printf("\t %10.6lf %10.6f\n",alphaMean[ii]/logit_factor,sqrt(V[ii+ii*NCOVAR])/logit_factor);
-//			printf("\n");
-//			for (int ii=NSUBTYPES;ii<NCOVAR;ii++)
-//				printf("\t %10.6lf %10.6f\n",alphaMean[ii]/logit_factor,sqrt(V[ii+ii*NCOVAR])/logit_factor);
-//			printf("\n");
-//			for (int ii=0;ii<NCOVAR;ii++) {
-//				for (int jj=0;jj<=ii;jj++)
+			printf("%f \n",beta);
+			for (int ii=0;ii<NCOV_FIX;ii++)
+				printf("%f \n",fixMean[ii]);
+			for (int ii=0;ii<NSUBTYPES;ii++)
+				printf("\t %10.6lf %10.6f\n",alphaMean[ii]/logit_factor,sqrt(V[ii+ii*NCOVAR])/logit_factor);
+			printf("\n");
+			for (int ii=NSUBTYPES;ii<NCOVAR;ii++)
+				printf("\t %10.6lf %10.6f\n",alphaMean[ii]/logit_factor,sqrt(V[ii+ii*NCOVAR])/logit_factor);
+			printf("\n");
+			for (int ii=0;ii<NCOVAR;ii++) {
+				for (int jj=0;jj<=ii;jj++)
 //					printf("%6.3lf ",sqrt(V[jj + ii*NCOVAR]));
-//					printf("%6.3lf ",V[jj + ii*NCOVAR]/sqrt(V[ii + ii*NCOVAR]*V[jj + jj*NCOVAR]));
-//				printf("\n");
-//			}
+					printf("%6.3lf ",V[jj + ii*NCOVAR]/sqrt(V[ii + ii*NCOVAR]*V[jj + jj*NCOVAR]));
+				printf("\n");
+			}
 			printf("\n");
 			fflush(NULL);
 			free(V);
@@ -591,7 +584,7 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 					cudaEventRecord(stop, 0);
 					cudaEventSynchronize(stop);
 					cudaEventElapsedTime(&time, start, stop);
-//					printf ("Time to compute_predictGPU: %f ms\n", time);
+					printf ("Time to compute_predictGPU: %f ms\n", time);
 				}
 
 				if (!(iter%PRNtITER))
@@ -603,7 +596,7 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 					cudaEventRecord(stop, 0);
 					cudaEventSynchronize(stop);
 					cudaEventElapsedTime(&time, start, stop);
-//					printf ("Time to compute_prbGPU: %f ms\n\n", time);
+					printf ("Time to compute_prbGPU: %f ms\n\n", time);
 				}
 			}
 			else {
@@ -614,15 +607,15 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 					cudaEventRecord(stop, 0);
 					cudaEventSynchronize(stop);
 					cudaEventElapsedTime(&time, start, stop);
-//					printf ("Time to compute_prbGPU: %f ms\n\n", time);
+					printf ("Time to compute_prbGPU: %f ms\n\n", time);
 				}
 			}
-			
+
 			for (int k=1;k<NDEPTH+1;k++) {
 				for (int j=1;j<NCOL+1;j++) {
 					for (int i=1;i<NROW+1;i++) {
 						IDX = idx(i,j,k);
-						int i1=i-1; 
+						int i1=i-1;
 						int j1=j-1;
 						int k1=k-1;
 						int IDX2 = idx2(i1,j1,k1);
@@ -637,12 +630,12 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 						//		Malpha[ii*RSIZE+IDX2] += (double)alpha[ii*TOTVOX+hostIdx[IDX]];
 								Mprb[ii*RSIZE+IDX2] += (double)prb[ii*TOTVOX+hostIdx[IDX]];
 							}
-						}						
+						}
 					}
 				}
 			}
 			TOTSAV++;
-		}	
+		}
 /*		if (!(iter%10000) && iter > 1) {
 			if (GPU) {
 				CUDA_CALL( cudaMemcpy(Z,deviceZ,NSUBS*TOTVOX*sizeof(float),cudaMemcpyDeviceToHost) );
@@ -658,23 +651,22 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 //			CUDA_CALL( cudaEventElapsedTime(&time, start, stop) );
 //			printf ("Time for the updateZGPU kernel: %f sec.\n", time/1000.0);
 
-//	fclose(fBF);
-	
-	unsigned int *Resid;
-	double *ResidMap;
-	Resid = (unsigned int *)calloc(NSUBS*TOTVOX,sizeof(unsigned int)); 	
-	ResidMap = (double *)calloc(RSIZE,sizeof(double)); 	
+	fclose(fBF);
+
+//	unsigned int *Resid;
+//	double *ResidMap;
+//	Resid = (unsigned int *)calloc(NSUBS*TOTVOX,sizeof(unsigned int));
+//	ResidMap = (double *)calloc(RSIZE,sizeof(double));
 	if (GPU) {
 		CUDA_CALL( cudaMemcpy(Z,deviceZ,NSUBS*TOTVOX*sizeof(float),cudaMemcpyDeviceToHost) );
-		CUDA_CALL( cudaMemcpy(Phi,devicePhi,NSUBS*TOTVOX*sizeof(float),cudaMemcpyDeviceToHost) );
+//		CUDA_CALL( cudaMemcpy(Phi,devicePhi,NSUBS*TOTVOX*sizeof(float),cudaMemcpyDeviceToHost) );
 //		cudaMemcpy(alpha,deviceAlpha,NSUBTYPES*TOTVOX*sizeof(float),cudaMemcpyDeviceToHost);
 		CUDA_CALL( cudaMemcpy(spatCoef,deviceSpatCoef,NCOVAR*TOTVOXp*sizeof(float),cudaMemcpyDeviceToHost) );
-		
-		CUDA_CALL( cudaMemcpy(Resid,deviceResid,NSUBS*TOTVOX*sizeof(unsigned int),cudaMemcpyDeviceToHost) );
+
+//		CUDA_CALL( cudaMemcpy(Resid,deviceResid,NSUBS*TOTVOX*sizeof(unsigned int),cudaMemcpyDeviceToHost) );
 	}
 //	save_iter(beta,fixMean,alphaMean,alphaPrec,SpatCoefMean,SpatCoefPrec,spatCoef,alpha,Z,Phi);
-
-	for (int isub=0;isub<NSUBS;isub++) {
+/*	for (int isub=0;isub<NSUBS;isub++) {
 		for (int k=1;k<NDEPTH+1;k++) {
 			for (int j=1;j<NCOL+1;j++) {
 				for (int i=1;i<NROW+1;i++) {
@@ -686,7 +678,7 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 					if (msk[IDX]) {
 						double tmp = (double)Resid[isub*TOTVOX+hostIdx[IDX]]/(double)(MAXITER-BURNIN);
 						if (tmp > 0.05) {
-//							fprintf(fresid,"%d %d %d %d %lf\n",isub,i,j,k,tmp);
+							fprintf(fresid,"%d %d %d %d %lf\n",isub,i,j,k,tmp);
 							ResidMap[IDX2]++;
 						}
 					}
@@ -695,23 +687,23 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 		}
 	}
 	free(Resid);
-//	fclose(out);
-//	fclose(out2);
-//	fclose(fresid);
+	fclose(out);
+	fclose(out2);
+	fclose(fresid);*/
 
 	char *RR = (char *)calloc(500,sizeof(char));
 	S = (char *)calloc(100,sizeof(char));
 //	S = strcpy(S,"ResidMap.nii");
 //	int rtn = write_nifti_file(NROW,NCOL,NDEPTH,1,S,S,ResidMap);
-	free(ResidMap);
-
+//	free(ResidMap);
+	int rtn;
 //	RR = strcpy(RR,"gzip -f ");
 //	RR = strcat(RR,(const char *)S);
 //	rtn = system(RR);
-	
-//	printf("TOTSAV = %d\n",TOTSAV);
 
-	out = fopen("Qhat.dat","w");	
+	printf("TOTSAV = %d\n",TOTSAV);
+
+	out = fopen("Qhat.dat","w");
 	for (int isub=0;isub<NSUBS;isub++) {
 		fprintf(out,"%d ",isub);
 		double max = -DBL_MAX;
@@ -742,10 +734,10 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 			}
 		}
 		fprintf(out,"%d\n",truetype);fflush(NULL);
-	}			
+	}
 	fclose(out);
 
-	out = fopen("Qhat2.dat","w");	
+	out = fopen("Qhat2.dat","w");
 	for (int isub=0;isub<NSUBS;isub++) {
 		fprintf(out,"%d ",isub);
 		double max = -DBL_MAX;
@@ -785,7 +777,7 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 	char *char_ii;
 	char_ii = (char *)calloc(3,sizeof(char));
 //	S = (char *)calloc(100,sizeof(char));
-	
+
 
 	for (int k=0;k<NDEPTH;k++) {
 		for (int j=0;j<NCOL;j++) {
@@ -794,14 +786,14 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 				MWM[IDX] /= (double)TOTSAV;
 //				fprintf(fWM,"%f ",MWM[IDX]);
 			}
-		}	
+		}
 	}
-	
+
 /*	for (int k=1;k<NDEPTH+1;k++) {
 		for (int j=1;j<NCOL+1;j++) {
 			for (int i=1;i<NROW+1;i++) {
 				IDX = idx(i,j,k);
-				int i1=i-1; 
+				int i1=i-1;
 				int j1=j-1;
 				int k1=k-1;
 				int IDX2 = idx2(i1,j1,k1);
@@ -815,20 +807,20 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 						Malpha[ii*RSIZE+IDX2] = NAN;
 						Mprb[ii*RSIZE+IDX2] = NAN;
 					}
-				}						
+				}
 			}
 		}
 	}*/
 
 //	fclose(fWM);
 	S = strcpy(S,"bWM.nii");
-//	rtn = write_nifti_file(NROW,NCOL,NDEPTH,1,S,S,MWM);
+	rtn = write_nifti_file(NROW,NCOL,NDEPTH,1,S,S,MWM);
 
 	RR = strcpy(RR,"gzip -f ");
 	RR = strcat(RR,(const char *)S);
-//	rtn = system(RR);
+	rtn = system(RR);
 
-/*	for (int ii=0;ii<NSUBTYPES;ii++) {	
+/*	for (int ii=0;ii<NSUBTYPES;ii++) {
 		for (int k=0;k<NDEPTH;k++) {
 			for (int j=0;j<NCOL;j++) {
 				for (int i=0;i<NROW;i++) {
@@ -842,7 +834,7 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 			}
 		}
 
-		
+
 		S = strcpy(S,"spatRE");
 		itoa(ii,char_ii);
 		S = strcat(S,char_ii);
@@ -858,7 +850,7 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 	}*/
 
 
-	for (int ii=0;ii<NCOVAR;ii++) {	
+	for (int ii=0;ii<NCOVAR;ii++) {
 		double *standCoef = (double *)calloc(RSIZE,sizeof(double));
 		for (int k=0;k<NDEPTH;k++) {
 			for (int j=0;j<NCOL;j++) {
@@ -906,7 +898,7 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 					IDX = idx(ii,jj,kk);
 					int IDX2 = idx2(i,j,k);
 					if (msk[IDX])
-						SCov[ii*RSIZE + IDX2] = MCov[ii*RSIZE + IDX2]/sqrt(SCov[ii*RSIZE + IDX2]);					
+						SCov[ii*RSIZE + IDX2] = MCov[ii*RSIZE + IDX2]/sqrt(SCov[ii*RSIZE + IDX2]);
 				}
 			}
 		}*/
@@ -922,16 +914,15 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 
 		free(standCoef);
 	}
-	
-	
-	for (int ii=0;ii<NSUBTYPES;ii++) {	
+
+	for (int ii=0;ii<NSUBTYPES;ii++) {
 		for (int k=0;k<NDEPTH;k++) {
 			for (int j=0;j<NCOL;j++) {
 				for (int i=0;i<NROW;i++) {
 					IDX = idx2(i,j,k);
 					double tmp = Mprb[ii*RSIZE+IDX]/(double)TOTSAV;
 					Mprb[ii*RSIZE+IDX] = tmp;
-					
+
 //					tmp = Malpha[ii*RSIZE+IDX]/(double)TOTSAV;
 //					Malpha[ii*RSIZE+IDX] = tmp;
 				}
@@ -954,27 +945,26 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 		rtn = write_nifti_file(NROW,NCOL,NDEPTH,1,S,S,&Malpha[ii*RSIZE]);*/
 	}
 
-
-//	FILE *fdic;
-//	fdic = fopen("DIC.dat","w");
+	FILE *fdic;
+	fdic = fopen("DIC.dat","w");
 	ED /= (double)TOTSAV;
 	ED *= -2.0;
 	double DE = compute_prb_DIC(MCov,MWM,covar,data,msk,RSIZE);
 	double PD = ED - DE;
-//	printf("DE = %lf ED = %lf PD = %lf DIC = %lf\n",DE,ED,PD,DE + 2*PD);
-//	fprintf(fdic,"DE = %lf ED = %lf PD = %lf DIC = %lf\n",DE,ED,PD,DE + 2*PD);
-	
-//	fclose(fdic);
+	printf("DE = %lf ED = %lf PD = %lf DIC = %lf\n",DE,ED,PD,DE + 2*PD);
+	fprintf(fdic,"DE = %lf ED = %lf PD = %lf DIC = %lf\n",DE,ED,PD,DE + 2*PD);
+
+	fclose(fdic);
 	free(S);
 	free(RR);
 	free(char_ii);
-	
+
 	free(alphaMean);
 	free(alphaPrec);
 
 	free(MCov);
 	free(SCov);
-	
+
 	free(MWM);
 
 
@@ -983,11 +973,11 @@ void mcmc(float *covar,float *covar_fix,unsigned char *data,float *WM,unsigned c
 		CUDA_CALL( cudaFree(deviceZ) );
 		CUDA_CALL( cudaFree(devicePrb) );
 		CUDA_CALL( cudaFree(devicePredict) );
-		CUDA_CALL( cudaFree(devicePhi) );
+//		CUDA_CALL( cudaFree(devicePhi) );
 	}
 	free(prb);
-	free(predict);		
-	free(spatCoef);	
+	free(predict);
+	free(spatCoef);
 	free(SpatCoefMean);
 	free(SpatCoefPrec);
 	free(Mprb);
@@ -1029,15 +1019,15 @@ void updateAlphaMean(float *alphaMean,float *Z,float *spatCoef,float beta,float
 					if (msk[IDX]) {
 						tmp = (double)Z[isub*TOTVOX+hostIdx[IDX]] - (double)beta*(double)WM[hostIdx[IDX]];
 						for (int ii=0;ii<NCOVAR;ii++)
-							tmp -= (double)spatCoef[ii*TOTVOXp+hostIdxSC[IDX]]*(double)covar[ii*NSUBS+isub]; 
+							tmp -= (double)spatCoef[ii*TOTVOXp+hostIdxSC[IDX]]*(double)covar[ii*NSUBS+isub];
 						sum += tmp;
 					}
 				}
 			}
 		}
-			
+
 		for (int ii=0;ii<NCOVAR;ii++)
-			tmpmean[ii] += sum*(double)covar[ii*NSUBS+isub];		
+			tmpmean[ii] += sum*(double)covar[ii*NSUBS+isub];
 	}
 
 	for (int i=0;i<NCOVAR;i++) {
@@ -1046,10 +1036,10 @@ void updateAlphaMean(float *alphaMean,float *Z,float *spatCoef,float beta,float
 	}
 
 	rmvnorm(tmpmean,V,NCOVAR,mean,seed,0);
-		
+
 	for (int i=0;i<NCOVAR;i++)
 		alphaMean[i] = (float)tmpmean[i];
-	
+
 	CUDA_CALL( cudaMemcpy(deviceAlphaMean,alphaMean,NCOVAR*sizeof(float),cudaMemcpyHostToDevice) );
 
 	free(tmpmean);
@@ -1063,7 +1053,7 @@ void updateAlphaPrec(float *alphaPrec,float *alpha,float alphaMean,float prior_a
 	float c = 0;
 	float y,t;
 //	float temp = 0;
-	
+
 	Alpha = 0.5*TOTVOX + prior_alphaPrec_A;
 	beta = 0;
 
@@ -1073,18 +1063,18 @@ void updateAlphaPrec(float *alphaPrec,float *alpha,float alphaMean,float prior_a
 				int IDX = idx(i,j,k);
 				if (msk[IDX]) {
 					tmp = (alpha[hostIdx[IDX]] - alphaMean);
-					y = tmp*tmp - c;	
+					y = tmp*tmp - c;
 					t = beta + y;
 					c = (t - beta) - y;
 					beta = t;
-//					beta += tmp*tmp;	
+//					beta += tmp*tmp;
 				}
 			}
 		}
 	}
 	beta = 0.5*beta + prior_alphaPrec_B;
 	*alphaPrec = (float)rgamma((double)Alpha,(double)beta,seed);
-	
+
 }
 
 
@@ -1094,13 +1084,13 @@ void updateBeta(float *beta,float *Z,float *alpha,float *WM,unsigned char *msk,f
 	float c = 0,cm=0;
 	float y,t,ym,tm;
 //	float temp = 0,tempm=0;
-		
+
 	for (int k=1;k<NDEPTH+1;k++) {
 		for (int j=1;j<NCOL+1;j++) {
 			for (int i=1;i<NROW+1;i++) {
 				int IDX = idx(i,j,k);
 				if (msk[IDX]) {
-					y = WM[hostIdx[IDX]]*WM[hostIdx[IDX]] - c;	
+					y = WM[hostIdx[IDX]]*WM[hostIdx[IDX]] - c;
 					t = var + y;
 					c = (t - var) - y;
 					var = t;
@@ -1111,7 +1101,7 @@ void updateBeta(float *beta,float *Z,float *alpha,float *WM,unsigned char *msk,f
 						for (int ii=0;ii<NSUBTYPES;ii++)
 							tmp -= (spatCoef[ii*TOTVOXp+hostIdxSC[IDX]])*covar[ii*NSUBS+isub];
 						for (int ii=NSUBTYPES;ii<NCOVAR;ii++)
-							tmp -= spatCoef[ii*TOTVOXp+hostIdxSC[IDX]]*covar[ii*NSUBS+isub]; 
+							tmp -= spatCoef[ii*TOTVOXp+hostIdxSC[IDX]]*covar[ii*NSUBS+isub];
 					}
 					tmp *= WM[hostIdx[IDX]];
 					ym = tmp - cm;
@@ -1134,7 +1124,7 @@ void updateBeta(float *beta,float *Z,float *alpha,float *WM,unsigned char *msk,f
 void updateAlpha(float *alpha,float alphaMean,float alphaPrec,float *Z,unsigned char *msk,float beta,float *WM,float *spatCoef,float *covar,int grp,unsigned long *seed)
 {
 	float mean=0,var=0,N=0;
-	
+
 	for (int k=1;k<NDEPTH+1;k++) {
 		for (int j=1;j<NCOL+1;j++) {
 			for (int i=1;i<NROW+1;i++) {
@@ -1146,14 +1136,14 @@ void updateAlpha(float *alpha,float alphaMean,float alphaPrec,float *Z,unsigned
 						if (covar[grp*NSUBS+isub]) {
 							N++;
 							mean += (Z[isub*TOTVOX+hostIdx[IDX]] - beta*WM[hostIdx[IDX]]);
-//							mean -= spatCoef[grp][i][j][k]*covar[grp*NSUBS+isub]; 
+//							mean -= spatCoef[grp][i][j][k]*covar[grp*NSUBS+isub];
 							for (int ii=0;ii<NCOVAR;ii++)
-								mean -= spatCoef[ii*TOTVOXp+hostIdxSC[IDX]]*covar[ii*NSUBS+isub]; 
+								mean -= spatCoef[ii*TOTVOXp+hostIdxSC[IDX]]*covar[ii*NSUBS+isub];
 						}
 					}
 					var = 1./(N + alphaPrec);
 					mean *= var;
-					double tmp = rnorm((double)mean,sqrt((double)var),seed);			
+					double tmp = rnorm((double)mean,sqrt((double)var),seed);
 					alpha[hostIdx[IDX]] = (float)tmp;
 				}
 			}
@@ -1164,7 +1154,7 @@ void updateAlpha(float *alpha,float alphaMean,float alphaPrec,float *Z,unsigned
 
 
 void initializeAlpha(unsigned char *msk,float *alphaMean,float *alphaPrec,float *alpha,unsigned long *seed){
-	
+
 	for (int ii=0;ii<NSUBTYPES;ii++) {
 
 /*		if (RESTART) {
@@ -1203,8 +1193,8 @@ void initializeAlpha(unsigned char *msk,float *alphaMean,float *alphaPrec,float
 }
 
 void initializeZ(float *Z,unsigned char *data,unsigned char *msk,unsigned long *seed){
-	
-		
+
+
 		for (int k=1;k<NDEPTH+1;k++) {
 			for (int j=1;j<NCOL+1;j++) {
 				for (int i=1;i<NROW+1;i++) {
@@ -1225,8 +1215,8 @@ void initializeZ(float *Z,unsigned char *data,unsigned char *msk,unsigned long *
 }
 
 void initializePhi(float *Phi,unsigned char *data,unsigned char *msk,float df,unsigned long *seed){
-	
-		
+
+
 		for (int k=1;k<NDEPTH+1;k++) {
 			for (int j=1;j<NCOL+1;j++) {
 				for (int i=1;i<NROW+1;i++) {
@@ -1246,9 +1236,9 @@ void initializePhi(float *Phi,unsigned char *data,unsigned char *msk,float df,un
 unsigned char *get_neighbors()
 {
 	unsigned char *nbrs;
-	
+
 	nbrs = (unsigned char *)calloc(NSIZE,sizeof(unsigned char));
-		
+
 	return nbrs;
 }
 
@@ -1258,7 +1248,7 @@ void save_iter(float beta,float *fixMean,float *alphaMean,float *alphaPrec,float
 {
 	FILE *out;
 	size_t rtn;
-	
+
 	out  = fopen("last_parms.dat","w");
 
 	fprintf(out,"%f ",beta);
@@ -1296,13 +1286,13 @@ void restart_iter(float *beta,float *fixMean,float *alphaMean,float *alphaPrec,f
 
 	rtn = fread(spatCoef,sizeof(float),NCOVAR*TOTVOXp,in);
 	fclose(in);
-	
+
 	in = fopen("last_Z.dat","r");
 	rtn = fread(Z,sizeof(float),NSUBS*TOTVOX,in);
-	
+
 	in = fopen("last_Phi.dat","r");
 	rtn = fread(Phi,sizeof(float),NSUBS*TOTVOX,in);
-	
+
 	fclose(in);
 
 	if (GPU) {
@@ -1310,10 +1300,10 @@ void restart_iter(float *beta,float *fixMean,float *alphaMean,float *alphaPrec,f
 		CUDA_CALL( cudaMemcpy(deviceZ,Z,NSUBTYPES*TOTVOX*sizeof(float),cudaMemcpyHostToDevice) );
 		CUDA_CALL( cudaMemcpy(devicePhi,Phi,NSUBTYPES*TOTVOX*sizeof(float),cudaMemcpyHostToDevice) );
 	}
-		
+
 	printf("RESTART = %d BURNIN = %d\n",RESTART,BURNIN);
 
-	
+
 	in = fopen("last_parms.dat","r");
 	rtn = fscanf(in,"%f ",&dta);
 	*beta = dta;