 #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;

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 compute the rank-1 Cholesky update/downdate. Note that this is done in place.
__device__ __host__

	int current_iter;
};

// Base class for an individual data point's characteristic, i.e., chi_i
class Characteristic {
public:
	// constructor
	__device__ __host__ Characteristic(curandState& localState, int iter, double rate) :
	state(localState), current_iter(iter), target_rate(rate), decay_rate(0.667);
	__device__ __host__ virtual ~Characteristic() {}

	// compute the conditional log-posterior density of the measurements given the characteristic
	__device__ __host__ virtual double logdensity_meas(double* chi, double* meas, double* meas_unc, int pchi, int mfeat)
	{
		return 0.0;
	}

	// compute the conditional log-posterior dentity of the characteristic given the population parameter
	__device__ __host__ virtual double logdensity_pop(double* chi, double* theta, int pchi, int dim_theta)
	{
		return 0.0;
	}

	// propose a new value for the characteristic
	__device__ __host__ virtual double* Propose(double* chi, double* cholfact, double* snorm_deviate, double* scaled_proposal, int pchi)
	{
		// get the unit proposal
		for (int j=0; j<pchi; j++) {
			snorm_deviate[j] = curand_normal_double(&state);
		}

		// propose a new chi value
		double proposed_chi[pchi];
		int cholfact_index = 0;
		for (int j=0; j<pchi; j++) {
			double scaled_proposal_j = 0.0;
			for (int k=0; k<(j+1); k++) {
				// transform the unit proposal to the centered proposal, drawn from a multivariate normal.
				scaled_proposal_j += cholfact[cholfact_index] * snorm_deviate[k];
				cholfact_index++;
			}
			proposed_chi[j] = chi[j] + scaled_proposal_j;
			scaled_proposal[j] = scaled_proposal_j;
		}
		return proposed_chi;
	}

	// adapt the covariance matrix of the proposals for the characteristics
	__device__ __host__ virtual void AdaptProp(double* cholfact, double* snorm_deviate, double* scaled_proposal,
			double metro_ratio, int pchi)
	{
		double unit_norm = 0.0;
		for (int j=0; j<pchi; j++) {
			unit_norm += snorm_deviate[j] * snorm_deviate[j];
		}
		unit_norm = sqrt(unit_norm);
		double decay_sequence = 1.0 / pow((double) current_iter, decay_rate);
		double scaled_coef = sqrt(decay_sequence * fabs(metro_ratio - target_rate)) / unit_norm;
		for (int j=0; j<pchi; j++) {
			scaled_proposal[j] *= scaled_coef;
		}
		bool downdate = (metro_ratio < target_rate);
		// do rank-1 cholesky update to update the proposal covariance matrix
		CholUpdateR1(cholfact, scaled_proposal, pchi, downdate);
	}

	// decide whether to accept or reject the proposal based on the metropolist-hasting ratio
	__device__ __host__ bool AcceptProp(double logdens_prop, double logdens_current, double forward_dens,
			double backward_dens)
	{
		double lograt = logdens_prop - forward_dens - (logdens_current - backward_dens);
		lograt = min(lograt, 0.0);
		double ratio = exp(lograt);
		double unif = curand_uniform_double(&state);
		bool accept = (unif < ratio) && std::isfinite(ratio);
		return accept;
	}

protected:
	// random number generator state for this characteristic
	curandState& state;
	// MCMC sampler parameters
	int current_iter;
	double decay_rate;
	double target_rate;
};

#endif /* PARAMETERS_CUH_ */
