 //#include <algorithm>
 //#include <cstdlib>
 
-#include <vector>
-/*
- typedef thrust::device_vector<double> dvec;
- typedef std::vector<dvec> vdvec;
- 
- struct wrapvec
- {
- vdvec v;
- 
- wrapvec(int m, int n) : v(m)
- {
- for (int i=0; i<m; i++)	{
- v[i].reserve(n);
- }
- }
- 
- double** ptrs()
- {
- thrust::host_vector<double*> h_ptr(v.size());
- for (unsigned int i=0; i<v.size(); i++)
- h_ptr[i] = (double*) thrust::raw_pointer_cast(&(v[i][0]));
- 
- thrust::device_vector<double*> d_ptr = h_ptr;
- return (double**) thrust::raw_pointer_cast(&d_ptr[0]);
- }
- };
- */
 #define CUDA_CALL(x) do { if((x) != cudaSuccess) { \
 printf("Error at %s:%d\n",__FILE__,__LINE__); \
 return EXIT_FAILURE;}} while(0)


// Pointer to the current value of theta. This is stored in constant memory so that all the threads on the GPU
// can access the same theta quickly.
static __constant__ double c_theta[2];

// Initialize the random number generator state
__global__ void setup_kernel(curandState *state)
{
    int id = threadIdx.x + blockIdx.x * blockDim.x;
    /* Each thread gets same seed, a different sequence

	thrust::random::experimental::normal_distribution<double> snorm(0., 1.0);
    thrust::random::uniform_real_distribution<double> uniform(0.0, 1.0);
    
    /*
     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 = 2000; // # of data points
	int m = 5; // # of features per data point (i.e., # of points in the SED for the i^th source)
    int p = 3; // # of characteristics per data point.
    int dim_theta = 2 * p;
    
    double mu_norm = 8.5 * log(10);  // Average value of the natural logarithm of the SED normalization
    double sig_norm = 0.5 * log(10);  // 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;
    h_theta[2] = mu_temp;
    h_theta[3] = sig_norm * sig_norm;
    h_theta[4] = sig_beta * sig_beta;
    h_theta[5] = sig_temp * sig_temp;
    
    // Allocate memory for arrays on host
    thrust::host_vector<double> h_meas(n * m);  // The measurements, m values for each of n data points
	thrust::host_vector<double> h_meas_unc(n * m);  // The measurement uncertainties
    thrust::host_vector<double> h_chi(n * p);  // Unknown characteristics, p values for each of n data points
    
    // Cholesky factor of proposal distribution for each data point, used by the Metropolis algorithm. The proposal
    // covariance matrix for each data point is Covar_i = PropChol_i * transpose(PropChol_i).
    dim_cholfactor = p * p - ((p - 1) * p) / 2  // only need one of the off-diagonal terms
    thrust::host_vector<double> h_prop_cholfact(n * dim_chol_factor);
    
	// Create simulated characteristics and data
    thrust::host_vector<double> h_true_chi(n * p);
    double sigma_msmt = [2.2e-4, 3.3e-4, 5.2e-4, 3.7e-4, 2.2e-4];  // Standard deviation for the measurement errors
    double lambda = [70.0, 170.0, 250.0, 350.0, 500.0];  // frequency bands correspond to Herschel PACS and SPIRES wavelengths
    double frequencies[5];
    for (int k=0; k<m; k++) {
        frequencies[k] = 1e6 / lambda[i];
    }
    
	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
            double bbody_numer = 2.0 * hplanck * frequencies[k] * frequencies[k] * frequencies[k] / (clight * clight);
            double temperature = exp(h_true_chi[n * 2 + i]); // Grap the temperature for this data point
            double bbody_denom = exp(hplanck * frequencies[k] / (kboltz * temperature)) - 1.0;
            double bbody = bbody_numer / bbody_denom;
            // Grab the SED normalization and power-law index
            double normalization = exp(h_true_chi[i]);
            double beta = h_true_chi[n + i];
            
            double nu_over_nu0 = frequencies[k] / nu0
            double SED_ik = normalization * pow(nu_over_nu0, beta) * bbody;
            
			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);
		}
	}
    
    // Initialize the chis to their true values for now
    thrust::copy(h_true_chi.begin(), h_true_chi.end(), h_chi.begin());
    
    // Allocate memory for arrays on device and copy the values from the host
    thrust::device_vector<double> d_meas = h_meas;
    thrust::device_vector<double> d_meas_unc = h_meas_unc;
    thrust::device_vector<double> d_chi = h_chi;
    thrust::device_vector<double> d_prop_cholfact = h_prop_cholfact;
    thrust::device_vector<double> d_logdens(n);  // log posteriors for an individual data point
    
    thrust::fill(d_logdens.begin(), d_logdens.end(), -1.0e300);
