 //#include <algorithm>
 //#include <cstdlib>
 
-#include <vector>
-/*
- typedef thrust::device_vector<double> dvec;
- typedef std::vector<dvec> vdvec;
- 
- struct wrapvec
- {
- vdvec v;
- 
- wrapvec(int m, int n) : v(m)
- {
- for (int i=0; i<m; i++)	{
- v[i].reserve(n);
- }
- }
- 
- double** ptrs()
- {
- thrust::host_vector<double*> h_ptr(v.size());
- for (unsigned int i=0; i<v.size(); i++)
- h_ptr[i] = (double*) thrust::raw_pointer_cast(&(v[i][0]));
- 
- thrust::device_vector<double*> d_ptr = h_ptr;
- return (double**) thrust::raw_pointer_cast(&d_ptr[0]);
- }
- };
- */
 #define CUDA_CALL(x) do { if((x) != cudaSuccess) { \
 printf("Error at %s:%d\n",__FILE__,__LINE__); \
 return EXIT_FAILURE;}} while(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 id = threadIdx.x + blockIdx.x * blockDim.x;
    /* Each thread gets same seed, a different sequence


int main(void)
{
    // Random number generator and distribution needed to simulate some data
	// unsigned int seed = 98724732;
    unsigned int seed =123455;
