 
 // Decay rate for proposal scaling factor updates. The proposal scaling factors decay as 1 / niter^decay_rate.
 // This is gamma in the notation of Vihola (2012)
-static __constant__ double c_decay_rate = 2.0 / 3.0;
+static __constant__ double c_decay_rate = 0.66667;
 
 // Current iteration of the MCMC sampler, needed to calculate the decay sequence for the RAM algorithm. Is it
 // best to put this in constant memory on the GPU and update from the CPU every iteration?
 static __constant__ int c_current_iter = 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 main(void)
{
	// measurements
	int n = 2000; // # of items
	int m = 1; // # of features
	/*
     wrapvec d_features(m,n);
     wrapvec d_sigmas(m,n);
     */
	thrust::host_vector<double> h_features(n*m);
	thrust::host_vector<double> h_sigmas(n*m);
    curandState* devStates;

	unsigned int seed = 98724732;
	double sigma_msmt = 3.;
	static thrust::minstd_rand rng(seed);
	thrust::random::experimental::normal_distribution<double> dist(0., sigma_msmt);
    
	double mu_popn_true = 20.;
	double sigma_popn_true = 1.;
	// Create simulated data; copies to GPU dumbly!
	for (int i=0; i<n; i++) {
		for (int j=0; j<m; j++) {
			h_features[i*m+0] = mu_popn_true + dist(rng);
			h_sigmas[i*m+1] = sigma_popn_true;
		}
	}
	thrust::device_vector<double> d_features = h_features;
	thrust::device_vector<double> d_sigmas = h_sigmas;
    
	// Create array of hyperparameter values.
	// Here we condition on sigma_popn_true (for simple analytical result).
	// Currently ineffecient; should build host vector and copy over.
	int dim_theta = 2;
	int n_theta = 11;
	double mu_lo = mu_popn_true - 2*sigma_msmt/sqrt(n);
	double mu_hi = mu_popn_true + 2*sigma_msmt/sqrt(n);
	double dmu = (mu_hi - mu_lo)/(n_theta-1.);
	double mu;
	thrust::host_vector<double> h_theta(dim_theta*n_theta);
	for (int i=0; i<n_theta; i++) {
		mu = mu_lo + i*dmu;
		h_theta[i*dim_theta] = mu;
		h_theta[i*dim_theta+1] = sigma_popn_true;
	}
	thrust::device_vector<double> d_theta = h_theta;
    
	// To load a single set of thetas into constant memory, copy to a global:
	// cudaMemcpyToSymbol(c_theta, p_theta, d_theta.size() * sizeof(*p_theta));
    
	// Alloc mem for marginals for individuals, for each set of hyperparams.
	thrust::device_vector<double> d_marg(n*n_theta);

    // Cuda grid launch
    dim3 nThreads(256);
    dim3 nBlocks((n + nThreads.x-1) / nThreads.x);

		}
	}
    */
    cudaFree(devStates);  // Free up the memory on the GPU from the RNG states
    
	return 0;
}
