         return 2;
     }
 
+    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 *
                          sizeof(curandState)));
     // Initialize the random number generator states on the GPU
     setup_kernel<<<nBlocks,nThreads>>>(devStates);
 
-    // Wait until everything is done running on the GPU
+    // Wait until everything is done running on the GPU, make sure everything went OK
     CUDA_CALL(cudaDeviceSynchronize());
 
+    // Now run the MCMC sampler
+    
     /*
 	{
 		// log marginal likelhoods in parallel on all threads independently

		}
	}
    */
    cudaFree(devStates);
    
	return 0;
    
}


int main(void)
{
	int n = 2000; // # of data points
	int m = 1; // # of features per data point (e.g., # of points in the SED for the i^th source)
    int p = 1; // * of characteristics per data point

    // Random number generator and distribution needed to simulate some data
	unsigned int seed = 98724732;
	static thrust::minstd_rand rng(seed);
	thrust::random::experimental::normal_distribution<double> snorm(0., 1.0);
    
    // Population level parameters: theta, where theta parameterizes the distribution of chi
    int dim_theta = 2;
	double mu_popn_true = 20.;  // Average value of the chi values
	double sigma_popn_true = 1.;  // Standard deviation of the chi values
    thrust::host_vector<double> h_theta(dim_theta);  // Allocate memory on host
    h_theta[0] = mu_popn_true;
    h_theta[1] = sigma_popn_true * sigma_popn_true;  // theta = (mu,var)
    
    // Allocate memory for arrays on host
    thrust::host_vector<double> h_meas(n * m);  // The measurements, m values for each of n data points
	thrust::host_vector<double> h_meas_unc(n * m);  // The measurement uncertainties
    thrust::host_vector<double> h_chi(n * p);  // Unknown characteristics, p values for each of n data points
    
    // Scale of proposal distribution for each data point, used by the Metropolis algorithm. When p > 1 this will be
    // the Cholesky factor of the proposal covariance matrix.
    thrust::host_vector<double> h_jump_sigma(n);
    
	// Create simulated characteristics and data
    std::vector<double> true_chi(n * p);
    double sigma_msmt = 3.;  // Standard deviation for the measurement errors
    
	for (int i=0; i<n; i++) {
		for (int j=0; j<p; j++) {
            // First generate true value of the characteristics
			true_chi[i * p + j] = mu_popn_true + sigma_popn_true * snorm(rng);
            // Initialize the scale of the chi proposal distributions to just be the measurement uncertainty
            h_jump_sigma[i * m + j] = sigma_msmt;
        }
        for (int k=0; k<m; k++) {
            // Now generate measurements given the true characteristics
			h_meas_unc[i * m + k] = sigma_msmt;
            // Just assume E(meas|chi) = chi for now
            h_meas[i * m + k] = true_chi[i * p + k] + sigma_msmt * snorm(rng);
		}
	}
    
    // Allocate memory for arrays on device and copy the values from the host
    thrust::device_vector<double> d_meas = h_meas;
	thrust::device_vector<double> d_meas_unc = h_meas_unc;
    thrust::device_vector<double> d_chi = h_chi;
    thrust::device_vector<double> d_jump_sigma = h_jump_sigma;
    thrust::device_vector<double> d_logdens;  // log posteriors for an individual data point
    
	// Load a single set of thetas into constant memory
    double* p_theta = thrust::raw_pointer_cast(&h_theta[0]);
	CUDA_CALL(cudaMemcpyToSymbol(c_theta, p_theta, h_theta.size() * sizeof(*p_theta)));
    
    // Cuda grid launch
    dim3 nThreads(256);
    dim3 nBlocks((n + nThreads.x-1) / nThreads.x);
