             std::cout << h_theta[j] << " ";
         }
         std::cout << std::endl;
+         */
         // Compute Metropolis ratio
         double logdens_pop = 0.0;
         // right now we loop over the elements of chi because we assume that they are statistically independent, i.e.,

        double* p_theta_cholfact = thrust::raw_pointer_cast(&h_theta_cholfact[0]);
        double* p_scaled_proposal = thrust::raw_pointer_cast(&scaled_proposal[0]);
        CholUpdateR1(p_theta_cholfact, p_scaled_proposal, dim_theta, downdate);
        */
        // Save the theta values
        for (int j=0; j<dim_theta; j++) {
            thetafile << h_theta[j] << " ";

        // 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();
        thrust::device_vector<double>::iterator chi_iter_end = d_chi.begin();
        /*
        for (int j=0; j<dim_theta/2; j++) {
            thrust::advance(chi_iter_end, n);
            double proposed_mu_j = proposed_theta[j];
