         return 2;
     }
 
+    /*
+    Set up parallel random number generations on the GPU
+    */
     curandState* devStates;  // Create state object for random number generators on the GPU
     // Allocate memory on GPU for RNG states
     CUDA_CALL(cudaMalloc((void **)&devStates, nThreads.x * nBlocks.x *


	for (int i=0; i<ndata; i++) {
        // Loop over the data indices
        int diag_index = 0;
		for (int j=0; j<p; j++) {
            // Loop over the chi indices to generate true value of the characteristics
			chi[i + ndata * j] = theta[j] + sqrt(exp(theta[p + j])) * snorm(rng);
            // Initialize the covariance matrix of the chi proposal distribution to be 0.01 * identity(p)
            prop_cholfact[i + ndata * diag_index] = 0.01;
            diag_index += j + 2;
        }
        for (int k=0; k<mfeat; k++) {
            // Loop over the feature indices to generate measurements, given the characteristics

        double* p_meas = thrust::raw_pointer_cast(&d_meas[0]);
        double* p_meas_unc = thrust::raw_pointer_cast(&d_meas_unc[0]);
        double* p_chi = thrust::raw_pointer_cast(&d_chi[0]);
        double* p_prop_cholfact = thrust::raw_pointer_cast(&d_chi_cholfactor[0]);
        double* p_logdens = thrust::raw_pointer_cast(&d_logdens[0]);
        int current_iter = i + 1;
        update_chi<<<nBlocks,nThreads>>>(p_chi, p_meas, p_meas_unc, ndata, p_logdens, p_prop_cholfact, devStates, current_iter);
