Merge pull request #8 from nicholst/BigN

Fixed critical bugs for samples larger than n=300; Improved useability

.gitignore

@@ -1,3 +1,4 @@
 .DS_Store
 .*~
 __MACOSX
+#*

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;

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));
 		}