 static const double kboltz = 1.380658e-16;
 static __constant__ double c_2h_over_csqr = 1.4745002e-47;  // 2 * hplanck / clight^2, used in the Planck function
 static __constant__ double c_h_over_k = 4.7992157e-11;  // hplanck / kboltz
+
+// Global random number generator and distributions for generating random numbers on the host. The random number generator used
+// is the Mersenne Twister mt19937 from the BOOST library.
+
+boost::random::mt19937 rng;
+boost::random::normal_distribution<> snorm(0.0, 1.0); // Standard normal distribution
+boost::random::uniform_real_distribution<> uniform(0.0, 1.0); // Uniform distribution from 0.0 to 1.0
 
 // Function to return the model SED at a input frequency, given the characteristics chi = (norm, beta, temp)
 __device__

#include <iostream>
#include <math.h>
#include <vector>

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

    }
};

// 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 i=0; i<ndata; i++) {
        // Loop over the data indices
        int diag_index = 0;
		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);
            // Initialize the covariance matrix of the chi proposal distribution to be 0.01 * identity(p)
            prop_cholfact[i + ndata * diag_index] = 0.01;
            diag_index += j + 2;
        }
        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);
		}
	}
}

int main(void)
{
    /* TODO: Should use the boost random number generator libraries for doing this on the host */

    /*
     Population level parameters: theta, where theta parameterizes the distribution of chi.
