 #include <stdlib.h>
 //#include <math.h>
 #include <time.h>
-
 #include <cuda.h>
 
 /* Compile with:
 nvcc -arch compute_50 cuda_biasd_simpson.c -o cuda_biasd_simpson
 */
 
-#define maxiter 2000
+#define DBL_EPSILON 2.22045e-16
+#define DBL_MIN 4.94066e-324
 
 extern "C" void log_likelihood(int N, double * d, double ep1, double ep2, double sigma, double k1, double k2, double tau, double epsilon, double * ll);
 extern "C" double sum_log_likelihood(int N, double *d, double ep1, double ep2, double sigma, double k1, double k2, double tau, double epsilon);
 
 __global__ void kernel_loglikelihood(int N, double * d, double ep1, double ep2, double sigma, double k1, double k2, double tau, double epsilon, double * ll);
-__device__ double adaptive_quad(int idx, double * d, double ep1, double ep2, double sigma, double k1, double k2, double tau, double epsilon);
 __device__ double integrand(double f, double d, double ep1, double ep2, double sigma, double k1, double k2, double tau);
+__device__ void qgauss(double a, double b, double d, double ep1, double ep2, double sigma, double k1, double k2, double tau, double *out);
+__device__ double integrate(double d, double ep1, double ep2, double sigma, double k1, double k2, double tau,double epsilon);
+
 
 // ##########################################################
 // ##########################################################


}

__device__ double adaptive_quad(int idx, double * d, double ep1, double ep2, double sigma, double k1, double k2, double tau, double epsilon) {

	// Doing it this way to cache all of the non-midpoint integrand calculations in order to minimize the number of calculations that are performed.

	int i = 0;
	double ival = 0., ax = 0., ay = 0., mx = 0., my = 0., h = 0., s1 = 0., s2_left = 0., s2_right = 0.;
	double bx[maxiter], by[maxiter];
	bx[0] = 1.;

	// Calculate f(0), and f(1)
	ay = integrand(ax,d[idx],ep1,ep2,sigma,k1,k2,tau);
	by[0] = integrand(bx[0],d[idx],ep1,ep2,sigma,k1,k2,tau);

	while (i >= 0) {
		h = (bx[i] - ax)/2.;
		mx = ax + h;
		my = integrand(mx,d[idx],ep1,ep2,sigma,k1,k2,tau); // Calc f(mid-point)

		s1 = h/3. * (ay + 4.*my + by[i]);

		s2_left = h/6. * (ay + 4.*integrand(ax + h/2.,d[idx],ep1,ep2,sigma,k1,k2,tau) + my);
		s2_right = h/6. * (my + 4.*integrand(ax + h*3./2.,d[idx],ep1,ep2,sigma,k1,k2,tau) + by[i]);

		// Success
		// if (fabs(s1 - s2_left - s2_right) < 15.*epsilon){ // absolute wasn't working... had strange pattern at ln(L) < ~-20
		if (fabs(s1 - s2_left - s2_right)/(s2_left + s2_right) < 15.*epsilon){
			ival += s2_left + s2_right;
			ax = bx[i];
			ay = by[i];
			i -= 1;
		} else { // Failure, so add midpoint to list if
			i += 1;
			if (i < maxiter - 1) { // there is space left in buffer
				bx[i] = mx;
				by[i] = my;
			} else { // there's no space in the buffer
				bx[maxiter - 1] = mx;
				by[maxiter - 1] = my;
				i = maxiter - 1;
			}
		}
	}
	return ival;
}


__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) {

		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
	}
