 
 
 __global__
-void marginals(double *theta, int d, double **features, double **sigmas, int m, int n, double *marg)
+void marginals(double *theta, int dim_theta, int n_theta, double **features, double **sigmas, int m, int n, double *marg)
 {
 	int i = blockDim.x * blockIdx.x + threadIdx.x;
-	if (i<n)
+	int i_theta = blockDim.y * blockIdx.y + threadIdx.y;
+	if (i<n && i_theta<n_theta)
 	{
-		double mu = theta[0];
-		double var = theta[1];
+		double mu = theta[0+dim_theta*i_theta];
+		double var = theta[1+dim_theta*i_theta];
 		double x;
 		double rho;
 		marg[i] = 0;

			x = features[0][i] + c_rt2*sigmas[0][i]*c_absc[j];
			double t = x-mu;
			rho = exp(-t*t/(2.0*var))/c_rt2pi;
			marg[i] += c_wts[j]*rho;
		}
	}
}

	thrust::device_vector<double> d_marg(n);

	// theta (shd be broadcasted, too?)
	int dim_theta = 2;
	int n_theta = 100;
	thrust::device_vector<double> d_theta(dim_theta*n_theta);

	// init, etc
	d_theta[0] = 5;
	d_theta[1] = 3.;

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

	{
		// log marginal likelhoods in parallel on all threads independently
		double* p_marg = thrust::raw_pointer_cast(&d_marg[0]);
		double* p_theta = thrust::raw_pointer_cast(&d_theta[0]);
		double** p_features = d_features.ptrs();
