                                                 zsqrj, logdens_pop, thrust::plus<double>());
             thrust::advance(chi_iter_begin, ndata);
         }
-
-        double logdens_old = thrust::reduce(d_logdens.begin(), d_logdens.end());
-        double lograt = logdens_pop - logdens_old;
+        */
+
+        // get log-density of current value of theta
+        double logdens_old = thrust::reduce(d_logdens_pop.begin(), d_logdens_pop.end());
+
+        // get value of population log-density for proposed theta value for each chi value. note that
+        // this overwrites the value of d_logdens_pop.
+        double* p_theta = thrust::raw_pointer_cast(&proposed_theta[0]);
+        // copy the proposed theta to constant memory on the device
+        CUDA_CALL(cudaMemcpyToSymbol(c_theta, p_theta, proposed_theta.size() * sizeof(*p_theta)));
+        p_chi = thrust::raw_pointer_cast(&d_chi[0]);  // grab the pointers in case thrust changed the location
+        p_logdens_pop = thrust::raw_pointer_cast(&d_logdens_pop[0]);
+
+        g_logdens_pop<<<nBlocks,nThreads>>>(p_chi, ndata, p_logdens_pop);
+        CUDA_CALL(cudaDeviceSynchronize());
+
+        double logdens_prop = thrust::reduce(d_logdens_pop.begin(), d_logdens_pop.end());
+
+        double lograt = logdens_prop - logdens_old;
         lograt = std::min(lograt, 0.0);
         double ratio = exp(lograt);
         double unif = uniform(rng);

            meas_unc_i[k] = meas_unc[n * k + i];
        }

        double logdens_meas = logdensity_meas(meas_i, meas_unc_i, new_chi);
        double logdens_pop = logdensity_pop(new_chi, c_theta);
        double logdens_prop = logdens_meas + logdens_pop;

        // Compute the Metropolis ratio
        double ratio = logdens_prop - logdens[i] - logdens_meas;
        ratio = (ratio < 0.0) ? ratio : 0.0;
        ratio = exp(ratio);

        // Now randomly decide whether to accept or reject
        double unif_draw = curand_uniform_double(&localState);

        if (unif_draw < ratio) {
            // Accept this proposal, so save this parameter value and conditional log-posterior
            for (int j=0; j<c_p; j++) {
                chi[n * j + i] = new_chi[j];
            }
            logdens[i] = logdens_pop;
        }

        // Copy state back to global memory
        state[i] = localState;

        // Finally, adapt the Cholesky factor of the proposal distribution
        // TODO: check that the ratio is finite before doing this step
        double unit_norm = 0.0;
        for (int j=0; j<c_p; j++) {
            unit_norm += snorm_deviate[j] * snorm_deviate[j];
        }
        unit_norm = sqrt(unit_norm);
        double decay_sequence = 1.0 / pow((double) current_iter, c_decay_rate);
        double scaled_coef = sqrt(decay_sequence * fabs(ratio - c_target_rate)) / unit_norm;
        for (int j=0; j<c_p; j++) {
            scaled_proposal[j] *= scaled_coef;
        }
        bool downdate = (ratio < c_target_rate);
        CholUpdateR1(cholfactor, scaled_proposal, c_p, downdate);

        // copy cholesky factor for this data point back to global memory
        cholfact_index = 0;
        for (int j=0; j<dim_cholfactor; j++) {
            prop_cholfact[n * j + i] = cholfactor[j];
        }
	}
}

struct zsqr : public thrust::unary_function<double,double> {
    double mu, var;
    zsqr(double m, double v) : mu(m), var(v) {}

    __device__ __host__
    double operator()(double chi) {
        double chi_cent = chi - mu;
        double logdens_pop = -0.5 * log(var) - 0.5 * chi_cent * chi_cent / var;
        return logdens_pop;
    }
};

// Generate a data set and fill the chi, meas, and meas_sigma host vector with the values
void generate_data(int ndata, thrust::host_vector<double>& theta, thrust::host_vector<double>& chi, thrust::host_vector<double>& meas,


        CUDA_CALL(cudaDeviceSynchronize());
        h_chi = d_chi;
        h_logdens_pop = d_logdens_pop;
        h_logdens_meas = d_logdens_meas;

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

        h_debug_info = d_debug_info;
        std::cout << "Difference in logdens_pop: " << h_debug_info[0] << std::endl;
        std::cout << "Difference in logdens_meas: " << h_debug_info[1] << std::endl;
        std::cout << "Difference in logdens: " << h_debug_info[2] << std::endl;
		*/

        // Update the theta values

        /*
        thrust::device_vector<double>::iterator chi_iter_begin = d_chi.begin();

        double logdens_pop = 0.0;
