 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
+
+// constant memory contains array sizes and MCMC sampler parameters
+__constant__ int c_ndata; // # of data points
+__constant__ int c_mfeat; // # of measured features per data point
+__constant__ int c_pchi; // # of characteristics per data point
+__constant__ int c_dim_theta; // # dimension of population parameters
+__constant__ double c_target_rate; // MCMC sampler target acceptance rate
+__constant__ double c_decay_rate; // decay rate of robust adaptive metropolis algorithm
 
 // Function to compute the rank-1 Cholesky update/downdate. Note that this is done in place.
 __device__ __host__

	int pchi = 2;

public:
	// Constructor when storing the measurements in std::vector
	DataAugmentation(vecvec& meas, vecvec& meas_unc, dim3& nB, dim3& nT) : nBlocks(nB), nThreads(nT)
	{
		int size1 = meas.size();
		int size2 = meas[0].size();

		ndata = max(size1,size2); // assume that ndata < mfeat and figure out the values from the size
		mfeat = min(size1,size2); // of the arrays.

		_SetArraySizes();

		// copy input data to data members
		for (int j = 0; j < mfeat; ++j) {
			for (int i = 0; i < ndata; ++i) {
				if (size1 < size2) {
					h_meas[ndata * j + i] = meas[j][i];
					h_meas_unc[ndata * j + i] = meas_unc[j][i];
				} else {
					h_meas[ndata * j + i] = meas[i][j];
					h_meas_unc[ndata * j + i] = meas_unc[i][j];
				}
			}
		}
		// copy data from host to device
		d_meas = h_meas;
		d_meas_unc = h_meas_unc;

		thrust::fill(h_cholfact.begin(), h_cholfact.end(), 0.0);
		d_cholfact = h_cholfact;

#include <thrust/device_vector.h>
#include <thrust/host_vector.h>

#define CUDA_CHECK_RETURN(value) {											\
	cudaError_t _m_cudaStat = value;										\
	if (_m_cudaStat != cudaSuccess) {										\
		fprintf(stderr, "Error %s at line %d in file %s\n",					\
				cudaGetErrorString(_m_cudaStat), __LINE__, __FILE__);		\
		exit(1);															\
	} }

typedef std::vector<std::vector<double> > vecvec;
typedef hvector hvector;
typedef dvector dvector;
