 	cudaSetDevice(device);
 	cudaDeviceProp deviceProp;
 	cudaGetDeviceProperties(&deviceProp, device);
-	int threads = deviceProp.maxThreadsPerBlock/2.;
+	int threads = deviceProp.maxThreadsPerBlock/4.;
 	int blocks = (N+threads-1)/threads;
 
 	double * d_d;

	double result21 = 0.;

	int i;
	// Function values
	double fval1[10],fval2[10];
	for (i=0;i<10;i++){
		fval1[i] = integrand(center+halflength*xgk[i],d,ep1,ep2,sigma,k1,k2,tau);
		fval2[i] = integrand(center-halflength*xgk[i],d,ep1,ep2,sigma,k1,k2,tau);
	}

	// Evaluate gauss-10 and kronrod-21
	for (i=0;i<5;i++){
		result10 += wg[i]*(fval1[2*i+1] + fval2[2*i+1]);
		result21 += wgk[2*i+1]*(fval1[2*i+1] + fval2[2*i+1]);
		result21 += wgk[2*i]*(fval1[2*i] + fval2[2*i]);
	}

	double fc = integrand(center,d,ep1,ep2,sigma,k1,k2,tau);
	result21 += wgk[10]*fc;

	// // Removed error calculation.....
	// errors
	// double resabs = result21; // in this case resabs = result21 b/c it's all positive
	// double reskh = result21*.5;
	// double resasc = wgk[10]*fabs(fc-reskh);
	// for (i=0;i<10;i++){
		// resasc += wgk[i]*(fabs(fval1[i] - reskh) + fabs(fval2[i] - reskh));
	// }

	out[0] = result21*halflength;
	// resabs *= halflength;
	// resasc *= halflength;
	// out[1] = fabs(result10-result21)*halflength;
	// if (resasc != 0. && out[1] != 0.){
	// 	out[1] = resasc*fmin(1.,pow(200.*out[1]/resasc,1.5));
	// }
	// if (resabs > DBL_MIN/(50.*DBL_EPSILON)){
	// 	out[1] = fmax(resabs*(DBL_EPSILON*50.),out[1]);
	// }
}

__device__ double integrate(double d, double ep1, double ep2, double sigma, double k1, double k2, double tau,double epsilon){

	double out[2],integral = 0.,error = 0.;

	int i;
	for (i=0;i<20;i++){
		qgauss(0.05*(double)i,0.05*(double)(i+1),d,ep1,ep2,sigma,k1,k2,tau,out);
		integral += out[0];
		error += out[1];
	}
	return integral;
}


__global__ void kernel_loglikelihood(int N, double * d, double ep1, double ep2, double sigma, double k1, double k2, double tau, double epsilon, double * ll) {

	int idx = threadIdx.x + blockIdx.x*blockDim.x;

	if (idx < N) {
		double out;
		out = k2/(k1+k2) * exp(-k1*tau - .5 * pow((d[idx]-ep1)/sigma,2.)); // state 1
		out += k1/(k1+k2) * exp(-k2*tau - .5 * pow((d[idx]-ep2)/sigma,2.)); // state 2

		out += 2.*k1*k2/(k1+k2)*tau *integrate(d[idx],ep1,ep2,sigma,k1,k2,tau,epsilon); // both

		out = log(out) - .5 * log(2.* M_PI) - log(sigma); // prefactor
		ll[idx] = out; // transfer
	}
}


	double d[50] = {0.87755042,  0.90101722,  0.88297422,  0.90225072,  0.91185969,        0.88479424,  0.64257305,  0.23650566,  0.17532272,  0.24785572,        0.77647345,  0.12143413,  0.04994399,  0.19918067,  0.09625039,
        0.14283554,  0.30052487,  0.8937437 ,  0.90544194,  0.87350816,        0.62315481,  0.48258872,  0.77018322,  0.42989469,  0.69183523,        0.35556625,  0.90622313,  0.12529433,  0.74309849,  0.8860914 ,        0.8335358 ,  0.56208782,  0.45287218,  0.79373139,  0.42808399,        0.86643919,  0.70459052,  0.09161765,  0.53514735,  0.06578612,        0.09050594,  0.14923124,  0.8579178 ,  0.884698  ,  0.8745358 ,        0.89191605,  0.57743238,  0.80656044,  0.9069933 ,  0.65817311};

	int device = 1;
	double sum = 0;
	printf("Starting Test...\n");
	sum = sum_log_likelihood(device,50,d,0.,1.,.05,3.,8.,.1,1e-10);
	double truth = -45.312935111611729;

	printf("Real : %.20f\n",truth);
	printf("Calcd: %.20f\n",sum);
	printf("%%Diff: %.20f\n",abs((truth-sum)/truth)*100.);
}
*/
