 }
 	  
 
-int main(void)
+int main(int argc, char *argv[])
 {
 	// measurements
-	int n = 2000000; // # of items
+	int n = 10000000; // # of items
 	int m = 1; // # of features
-	/*
-	wrapvec d_features(m,n);
-	wrapvec d_sigmas(m,n);
-	*/
+
+	int devId = 0;
+
+	if (argc > 1) n = atoi(argv[1]);
+	if (argc > 2) devId = atoi(argv[2]);
+
+	cudaError_t err = cudaSetDevice(devId);
+	if (err != CUDA_SUCCESS) { return 1; }
+
+
 	thrust::host_vector<double> h_features(n*m);
 	thrust::host_vector<double> h_sigmas(n*m);
 

	// 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);

//#include <cstdlib>

#include <vector>
// constants for integration
#define N_GH 10

