 
 int main(void)
 {
-    /* TODO: Should use the boost random number generator libraries for doing this on the host */
+    /* TODO: Add command line options */
 
     /*
      Population level parameters: theta, where theta parameterizes the distribution of chi.


        CUDA_CALL(cudaDeviceSynchronize());
        h_chi = d_chi;
        h_logdens = d_logdens;

        // Compute Metropolis ratio
        //
        // right now we loop over the elements of chi because we assume that they are statistically independent, i.e.,
        // the theta array does not contain correlations among the element of chi. this is unrealistic, and we will need
        // to come up with a better way to do the transform + reduction: maybe use CUBLAS?
        thrust::device_vector<double>::iterator chi_iter_begin = d_chi.begin();

        double logdens_pop = 0.0;


    //std::ofstream chifile("chis.dat");
    std::ofstream thetafile("thetas.dat");
    std::ofstream chifile("chi0.dat");
    int mcmc_iter = 25000;
    int naccept_theta = 0;
    std::cout << "Running MCMC Sampler...." << std::endl;

    thrust::host_vector<double> h_logdens_meas;
    thrust::host_vector<double> h_logdens_pop;
    thrust::host_vector<int> h_naccept(ndata);
    thrust::fill(h_naccept.begin(), h_naccept.end(), 0);
    thrust::device_vector<int> d_naccept = h_naccept;

    thrust::device_vector<double> d_debug_info(3);
    thrust::host_vector<double> h_debug_info = d_debug_info;

    // Initialize the log-posteriors
    double* p_meas = thrust::raw_pointer_cast(&d_meas[0]); // CUDA kernels need to have the raw pointers
    double* p_meas_unc = thrust::raw_pointer_cast(&d_meas_unc[0]);
    double* p_chi = thrust::raw_pointer_cast(&d_chi[0]);
    double* p_logdens_meas = thrust::raw_pointer_cast(&d_logdens_meas[0]);
    double* p_logdens_pop = thrust::raw_pointer_cast(&d_logdens_pop[0]);
    g_logdens_pop<<<nBlocks,nThreads>>>(p_chi, ndata, p_logdens_pop);
    g_logdens_meas<<<nBlocks,nThreads>>>(p_chi, p_meas, p_meas_unc, ndata, p_logdens_meas);
    CUDA_CALL(cudaDeviceSynchronize());

    for (int i=0; i<mcmc_iter; i++) {
        // Now grab the pointers to the vectors, needed to run the kernel since it doesn't understand Thrust
        // We do this here because the thrust vectors are smart, and we want to make sure they don't reassign
        // memory for whatever reason. This is very cheap to do, so better safe than sorry.
        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_meas = thrust::raw_pointer_cast(&d_logdens_meas[0]);
        double* p_logdens_pop = thrust::raw_pointer_cast(&d_logdens_pop[0]);
        int* p_naccept = thrust::raw_pointer_cast(&d_naccept[0]);
        int current_iter = i + 1;

        /*
        std::cout << "current h_chi: " << std::endl;
        for (int j = 0; j < p; ++j) {
        	std::cout << h_chi[j*ndata] << " ";
		}
        std::cout << std::endl;
		*/

        // Update the characteristics on the GPU
        double* p_debug_info = thrust::raw_pointer_cast(&d_debug_info[0]);
        update_chi<<<nBlocks,nThreads>>>(p_chi, p_meas, p_meas_unc, ndata, p_logdens_meas, p_logdens_pop,
        		p_prop_cholfact, devStates, current_iter, p_naccept, p_debug_info);

        // Generate new theta in parallel with GPU calculation above
        thrust::host_vector<double> proposed_theta(dim_theta);
