             for (int k=0; k<(j+1); k++) {
                 // transform the unit proposal to the centered proposal, drawn from a multivariate normal
                 // or t-distribution.
-                
-                // DOUBLE CHECK THIS: DOES NOT LOOK RIGHT, WHY DOES K NOT APPEAR????
-                scaled_proposal_j += cholfactor[cholfact_index] * snorm_deviate[j];
+                scaled_proposal_j += cholfactor[cholfact_index] * snorm_deviate[k];
                 cholfact_index++;
             }
-            scaled_proposal_j = 0.0;
+            // scaled_proposal_j = 0.0;
             new_chi[j] = chi[n * j + i] + scaled_proposal_j;
             scaled_proposal[j] = scaled_proposal_j;
         }
 
         // Compute the conditional log-posterior of the proposed parameter values for this data point
-        double meas_i[m];
-        double meas_unc_i[m];
+        double meas_i[mfeat];
+        double meas_unc_i[mfeat];
         for (int k=0; k<c_m; k++) {
             // grab the measurements for this data point
             meas_i[k] = meas[n * k + i];
             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;
-        
+
+        double logdens_meas_prop = logdensity_meas(meas_i, meas_unc_i, new_chi);
+        double logdens_pop_prop = logdensity_pop(new_chi, c_theta);
+        double logdens_prop = logdens_meas_prop + logdens_pop_prop;
+
+        bool finite_logdens = isfinite(logdens_prop);
+
         // Compute the Metropolis ratio
-        double ratio = logdens_prop - logdens[i] - logdens_meas;
+        double logdens_old = logdens_pop[i] + logdens_meas[i];
+        double ratio = logdens_prop - logdens_old;
+
+        if (i == 0) {
+			debug_info[0] = logdens_meas_prop - logdens_meas[i];
+			debug_info[1] = logdens_pop_prop - logdens_pop[i];
+			debug_info[2] = ratio;
+		}
+
         ratio = (ratio < 0.0) ? ratio : 0.0;
         ratio = exp(ratio);
-        
-        p_mha_ratio[i] = ratio;
-        
+
         // Now randomly decide whether to accept or reject
         double unif_draw = curand_uniform_double(&localState);
-        
-        if (unif_draw < ratio) {
+
+        if ((unif_draw < ratio) && finite_logdens) {
             // 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;
+            logdens_pop[i] = logdens_pop_prop;
+            logdens_meas[i] = logdens_meas_prop;
+            ++naccept[i];
         }
 
         // 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];
-        }
-	}
-}
-
+        if (finite_logdens) {
+			// Finally, adapt the Cholesky factor of the proposal distribution
+			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
+			for (int j=0; j<dim_cholfactor; j++) {
+				prop_cholfact[n * j + i] = cholfactor[j];
+			}
+        }
+	}
+}
+
+// calculate the logdensity of chi|meas,theta for each chi on the device, needed for computing the
+// log-densities of the initial values
+__global__
+void g_logdens_meas(double* chi, double* meas, double* meas_unc, int ndata, double* logdens_meas)
+{
+	int i = blockDim.x * blockIdx.x + threadIdx.x;
+	if (i < ndata)
+	{
+		double this_chi[p];
+		double this_meas[mfeat];
+		double this_meas_unc[mfeat];
+		for (int j=0; j<p; j++) {
+			// grab the characteristics for this data point
+			this_chi[j] = chi[j * ndata + i];
+		}
+        for (int k=0; k<c_m; k++) {
+            // grab the measurements for this data point
+            this_meas[k] = meas[ndata * k + i];
+            this_meas_unc[k] = meas_unc[ndata * k + i];
+        }
+        logdens_meas[i] = logdensity_meas(this_meas, this_meas_unc, this_chi);
+	}
+}
+
+// calculate the logdensity of theta for each chi on the device, needed for updating theta
+__global__
+void g_logdens_pop(double* chi, int ndata, double* logdens_pop)
+{
+	int i = blockDim.x * blockIdx.x + threadIdx.x;
+	if (i < ndata)
+	{
+		double this_chi[p];
+		for (int j=0; j<p; j++) {
+			this_chi[j] = chi[j * ndata + i];
+		}
+		logdens_pop[i] = logdensity_pop(this_chi, c_theta);
+	}
+}
+
+/*
 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 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,
+		thrust::host_vector<double>& meas_sigma) {
+
+	double sigma_msmt[5] = {2.2e-4, 3.3e-4, 5.2e-4, 3.7e-4, 2.2e-4};  // Standard deviation for the measurement errors
+
+	for (int j = 0; j < mfeat; ++j) {
+		sigma_msmt[j] *= 1.0;
+	}
+
+	for (int i=0; i<ndata; 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
+			chi[i + ndata * j] = theta[j] + sqrt(exp(theta[p + j])) * snorm(rng);
+        }
+        for (int k=0; k<mfeat; k++) {
+            // Loop over the feature indices to generate measurements, given the characteristics
+            double bbody_numer = 2.0 * hplanck * nu[k] * nu[k] * nu[k] / (clight * clight);
+            double temperature = exp(chi[ndata * 2 + i]); // Grab the temperature for this data point
+            double bbody_denom = exp(hplanck * nu[k] / (kboltz * temperature)) - 1.0;
+            double bbody = bbody_numer / bbody_denom;
+            // Grab the SED normalization and power-law index
+            double normalization = exp(chi[i]);
+            double beta = chi[ndata + i];
+            // Compute the model SED
+            double nu_over_nu0 = nu[k] / nu0;
+            double SED_ik = normalization * pow(nu_over_nu0, beta) * bbody;
+            // Generate measured SEDs
+			meas_sigma[k * ndata + i] = sigma_msmt[k];
+            meas[k * ndata + i] = SED_ik + sigma_msmt[k] * snorm(rng);
+		}
+	}
+}
+
+// Initialize the chi and theta values. Right now this just initializes them to their true values. Also initialize the Cholesky factors
+// for the chi and theta proposals. The elements of the chi and theta proposals are initially independent.
+void initialize_parameters(int ndata, thrust::host_vector<double>& theta, thrust::host_vector<double>& chi,
+		thrust::host_vector<double>& theta_cholfactor, thrust::host_vector<double>& chi_cholfactor, thrust::host_vector<double>& true_chi) {
+
+	// For now, initialize chi values for drawing from their parent distribution.
+	for (int i = 0; i < ndata; ++i) {
+		for (int j = 0; j < p; ++j) {
+			double chi_mean = theta[j];
+			double chi_sig = sqrt(exp(theta[p + j]));
+			chi[j * ndata + i] = chi_mean + chi_sig * snorm(rng);
+		}
+	}
+
+	for (int j = 0; j < dim_theta; ++j) {
+		theta[j] = theta[j] + theta[j] * 0.2 * snorm(rng);
+	}
+
+	// Initialize the cholesky factor for the theta proposals
+    int diag_index = 0;
+    for (int j=0; j<dim_theta; j++) {
+    	// Initialize the covariance matrix of the theta proposal distribution to be 0.01 * identity(p)
+        theta_cholfactor[diag_index] = 0.01;
+        diag_index += j + 2;
+    }
+
+    // Initialize the cholesky factor for the chi proposals
+    for (int i=0; i<ndata; i++) {
+        int diag_index = 0;
+        for (int j=0; j<p; j++) {
+        	// Initialize the covariance matrix of the chi proposal distribution to be 0.01 * identity(p)
+        	chi_cholfactor[i + ndata * diag_index] = 0.01;
+        	diag_index += j + 2;
+        }
+    }
+}
+
+// Test the rank-1 cholesky update
+void test_CholUpdateR1() {
+	int dim_cf = 4;
+	int size_cf = 10; //dim_cf * dim_cf - ((dim_cf - 1) * dim_cf) / 2
+	double L[10] = {31.088, 7.759, 31.559, 7.632, 6.998, 29.739, 6.345, 7.080, 6.108, 29.015};
+	// true values of the updated cholesky factor for A = L * transpose(L) +/- v * transpose(v)
+	double Lup0[10] = {31.542, 4.726, 36.504, 6.898, 8.527, 29.748, 7.414, 3.262, 5.894, 30.270};
+	double Ldown0[10] = {30.628, 10.884, 25.278, 8.389, 4.944, 29.718, 5.246, 13.219, 6.581, 26.017};
+	double v[4] = {-5.331, 17.284, 3.691, -6.861};
+
+	double Ldown[10];
+	double Lup[10];
+	double vup[4];
+	double vdown[4];
+	for (int i=0; i<size_cf; i++) {
+		Lup[i] = L[i];
+		Ldown[i] = L[i];
+	}
+	for (int i=0; i<dim_cf; i++) {
+		vup[i] = v[i];
+		vdown[i] = v[i];
+	}
+
+	// test the update first
+	CholUpdateR1(Lup, vup, dim_cf, false);
+	CholUpdateR1(Ldown, vdown, dim_cf, true);
+	double fracdiff_up = 0.0;
+	double fracdiff_down = 0.0;
+	for (int i=0; i<size_cf; i++) {
+		std::cout << "Lup[i]: " << Lup[i] << ", Ldown[i]: " << Ldown[i] << std::endl;
+		double fracdiff = abs(Lup[i] - Lup0[i]) / Lup0[i];
+		fracdiff_up = max(fracdiff_up, fracdiff);
+		fracdiff = abs(Ldown[i] - Ldown0[i]) / Ldown0[i];
+		fracdiff_down = max(fracdiff_down, fracdiff);
+	}
+	std::cout << "test_CholUpdateR1:" << std::endl;
+	std::cout << "Maximum fractional difference for rank-1 update: " << fracdiff_up << std::endl;
+	std::cout << "Maximum fractional difference for rank-1 downdat: " << fracdiff_down << std::endl;
+}
 
 int main(void)
 {
-    /* TODO: Should use the boost random number generator libraries for doing this on the host */
-    
-    // Random number generator and distribution needed to simulate some data
-	// unsigned int seed = 98724732;
-    unsigned int seed = 123455;
-	static thrust::minstd_rand rng(seed);
-	thrust::random::experimental::normal_distribution<double> snorm(0., 1.0);
-    thrust::random::uniform_real_distribution<double> uniform(0.0, 1.0);
-    
+    /* TODO: Add command line options */
+
     /*
      Population level parameters: theta, where theta parameterizes the distribution of chi.
      For this application, chi = (log10(norm), beta, log10(temp)), where the model SED is a
      modified black body:
-    
+
         SED(nu) = norm * (nu / nu0)^beta * Bbody(nu,temp).
-     
+
      Here, Bbody(nu,temp) is the Planck function. For simplicity, right now I just assume that
      log10(norm), beta, and log10(temp) are independently drawn from different Gaussian distributions,
      although we should later include correlations among the parameters. So, for now, theta is just the
      collection of the mean and variances of the normal distributions.
     */
-    
-    int n = 3; // # of data points
-    
+
+    int ndata = 1000; // # of data points
+
     double mu_norm = 8.5 * log(10.0);  // Average value of the natural logarithm of the SED normalization
     double sig_norm = 0.5 * log(10.0);  // Standard deviation in the SED normalization
     double mu_beta = 2.0;  // Average value of the SED power-law index
     double sig_beta = 0.2;  // Standard deviation in the SED power-law index
     double mu_temp = log(15.0);  // Average value of the logarithm of the dust temperature
     double sig_temp = 0.2;  // Standard deviation in the logarithm of the SED temperature
-    
+
     thrust::host_vector<double> h_theta(dim_theta);  // Allocate memory on host
     h_theta[0] = mu_norm;
     h_theta[1] = mu_beta;

// Cuda Includes
#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <stdio.h>
#include <fstream>
#include <stdlib.h>
#include <iostream>
#include <math.h>
#include <vector>

// Boost includes
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/normal_distribution.hpp>
#include <boost/random/uniform_real_distribution.hpp>

//#include <thrust/sort.h>
//#include <thrust/copy.h>
