             // Grab the SED normalization and power-law index
             double normalization = exp(h_true_chi[i]);
             double beta = h_true_chi[n + i];
-            
+            // Compute the model SED
             double nu_over_nu0 = frequencies[k] / nu0
             double SED_ik = normalization * pow(nu_over_nu0, beta) * bbody;
-            
+            // Generate measured SEDs
 			h_meas_unc[k * n + i] = sigma_msmt[k];
-            // Just assume E(meas|chi) = chi for now
             h_meas[k * n + i] = SED_ik + sigma_msmt[k] * snorm(rng);
 		}
 	}

    
	for (int i=0; i<n; i++) {
        // Loop over the data indices
		for (int j=0; j<p; j++) {
            // Loop over the chi indices to generate true value of the characteristics
			h_true_chi[i + n * j] = h_theta[j] + sqrt(h_theta[p + j]) * snorm(rng);
            // Initialize the scale of the chi proposal distributions to just be some small constant value
            h_jump_sigma[i + n * j] = 0.01;
        }
        for (int k=0; k<m; k++) {
            // Loop over the feature indices to generate measurements, given the characteristics


// Perform the MHA update on the n parameters, done in parallel on the GPU. This is what the kernel does.
__global__
void update_chi(double* chi, double* meas, double* meas_unc, int n, double* logdens, double* prop_cholfact,
                curandState* state, int current_iter)
{
	int i = blockDim.x * blockIdx.x + threadIdx.x;
