 // Function to compute the conditional log-density of the chi values given theta
 __device__ __host__
 double logdensity_pop(double* chi, double* theta) {
+#ifdef __CUDA_ARCH__
+    dim_chi = c_p;  // function is called from device, so need to access constant memory on the GPU
+#else
+    dim_chi = p;
+#endif
     // right now this is just an independent p-dimensional gaussian distribution, for simplicity
     double logdensity = 0.0;
-    for (int j=0; j<c_p; j++) {
+    for (int j=0; j<dim_chi; j++) {
         double mu_j = theta[j];
-        double var_j = theta[j + c_p];
+        double var_j = theta[j + dim_chi];
         double centered_chi_j = chi[j] - mu_j;
         logdensity += -0.5 * log(var_j) - 0.5 * centered_chi_j * centered_chi_j / var_j;
     }
     return logdensity;
+}
+
+// Function to compute the rank-1 Cholesky update/downdate
+__device__ __host__
+void CholUpdateR1(double* cholfactor, double* v, bool downdate) {
+#ifdef __CUDA_ARCH__
+    dim_chi = c_p;  // function is called from device, so need to access constant memory on the GPU
+#else
+    dim_chi = p;
+#endif
+    
 }
 
 // Initialize the random number generator state

    beta = chi[1];
    temperature = exp(chi[2]);
    
    double logdensity = 0.0
    for (int k=0; k<c_m; k++) {
        model_sed = modified_bbody(c_nu[k], normalization, beta, temperature);
        double stnd_meas = (meas[k] - model_sed) / meas_unc[k]
        logdensity += -0.5 * stnd_meas * stnd_meas;
    }
    return logdensity;


        // 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(unitnorm);
        double decay_sequence = 1.0 / pow(current_iter, c_decay_rate);
        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, 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];
        }
	}
}

