 // Constants used for the model SED
 static __constant__ double c_nu0 = 2.3e11;  // reference frequency, nu0 = 230 GHz
 static const double nu0 = 2.3e11;
-static const int m = 5;  // the number of measured features per data point
+static const int mfeat = 5;  // the number of measured features per data point
 static __constant__ int c_m = 5;
     // observational frequencies, corresponding to the Herschel PACS and SPIRES bands
-static const double nu[m] = {4.286e12, 1.765e12, 1.200e12, 8.571e11, 6.000e11};
-static __constant__ double c_nu[m] = {4.286e12, 1.765e12, 1.200e12, 8.571e11, 6.000e11};
+static const double nu[mfeat] = {4.286e12, 1.765e12, 1.200e12, 8.571e11, 6.000e11};
+static __constant__ double c_nu[mfeat] = {4.286e12, 1.765e12, 1.200e12, 8.571e11, 6.000e11};
 static const double hplanck = 6.6260755e-27;  // Planck's constant, in cgs
 static const double clight = 2.997925e10;
 static const double kboltz = 1.380658e-16;
 static __constant__ double c_2h_over_csqr = 1.4745002e-47;  // 2 * hplanck / clight^2, used in the Planck function
 static __constant__ double c_h_over_k = 4.7992157e-11;  // hplanck / kboltz
+
+// Global random number generator and distributions for generating random numbers on the host. The random number generator used
+// is the Mersenne Twister mt19937 from the BOOST library.
+boost::random::mt19937 rng;
+boost::random::normal_distribution<> snorm(0.0, 1.0); // Standard normal distribution
+boost::random::uniform_real_distribution<> uniform(0.0, 1.0); // Uniform distribution from 0.0 to 1.0
 
 // Function to return the model SED at a input frequency, given the characteristics chi = (norm, beta, temp)
 __device__

// Cuda Includes
#include "cuda_runtime.h"
#include "device_launch_parameters.h"

            for (int k=0; k<(j+1); k++) {
                // transform the unit proposal to the centered proposal, drawn from a multivariate normal
                // or t-distribution.
                scaled_proposal_j += cholfactor[cholfact_index] * snorm_deviate[k];
                cholfact_index++;
            }
            // scaled_proposal_j = 0.0;
            new_chi[j] = chi[n * j + i] + scaled_proposal_j;
            scaled_proposal[j] = scaled_proposal_j;
        }

        // Compute the conditional log-posterior of the proposed parameter values for this data point
        double meas_i[mfeat];
        double meas_unc_i[mfeat];
        for (int k=0; k<c_m; k++) {
            // grab the measurements for this data point
            meas_i[k] = meas[n * k + i];
            meas_unc_i[k] = meas_unc[n * k + i];
        }

        double logdens_meas_prop = logdensity_meas(meas_i, meas_unc_i, new_chi);
        double logdens_pop_prop = logdensity_pop(new_chi, c_theta);
        double logdens_prop = logdens_meas_prop + logdens_pop_prop;

        bool finite_logdens = isfinite(logdens_prop);

        // Compute the Metropolis ratio
        double logdens_old = logdens_pop[i] + logdens_meas[i];
        double ratio = logdens_prop - logdens_old;

        if (i == 0) {
			debug_info[0] = logdens_meas_prop - logdens_meas[i];
			debug_info[1] = logdens_pop_prop - logdens_pop[i];
			debug_info[2] = ratio;
		}

        ratio = (ratio < 0.0) ? ratio : 0.0;
        ratio = exp(ratio);

        // Now randomly decide whether to accept or reject
        double unif_draw = curand_uniform_double(&localState);

        if ((unif_draw < ratio) && finite_logdens) {
            // Accept this proposal, so save this parameter value and conditional log-posterior
            for (int j=0; j<c_p; j++) {
                chi[n * j + i] = new_chi[j];
            }
            logdens_pop[i] = logdens_pop_prop;
            logdens_meas[i] = logdens_meas_prop;
            ++naccept[i];
        }

        // Copy state back to global memory
        state[i] = localState;

        if (finite_logdens) {
			// Finally, adapt the Cholesky factor of the proposal distribution
			double unit_norm = 0.0;
			for (int j=0; j<c_p; j++) {
				unit_norm += snorm_deviate[j] * snorm_deviate[j];
			}
			unit_norm = sqrt(unit_norm);
			double decay_sequence = 1.0 / pow((double) current_iter, c_decay_rate);
			double scaled_coef = sqrt(decay_sequence * fabs(ratio - c_target_rate)) / unit_norm;
			for (int j=0; j<c_p; j++) {
				scaled_proposal[j] *= scaled_coef;
			}
			bool downdate = (ratio < c_target_rate);
			CholUpdateR1(cholfactor, scaled_proposal, c_p, downdate);

			// copy cholesky factor for this data point back to global memory
			for (int j=0; j<dim_cholfactor; j++) {
				prop_cholfact[n * j + i] = cholfactor[j];
			}
        }
	}
}

// calculate the logdensity of chi|meas,theta for each chi on the device, needed for computing the
// log-densities of the initial values
__global__
void g_logdens_meas(double* chi, double* meas, double* meas_unc, int ndata, double* logdens_meas)
{
	int i = blockDim.x * blockIdx.x + threadIdx.x;
	if (i < ndata)
	{
		double this_chi[p];
		double this_meas[mfeat];
		double this_meas_unc[mfeat];
		for (int j=0; j<p; j++) {
			// grab the characteristics for this data point
			this_chi[j] = chi[j * ndata + i];
		}
        for (int k=0; k<c_m; k++) {
            // grab the measurements for this data point
            this_meas[k] = meas[ndata * k + i];
            this_meas_unc[k] = meas_unc[ndata * k + i];
        }
        logdens_meas[i] = logdensity_meas(this_meas, this_meas_unc, this_chi);
	}
}

// calculate the logdensity of theta for each chi on the device, needed for updating theta
__global__
void g_logdens_pop(double* chi, int ndata, double* logdens_pop)
{
	int i = blockDim.x * blockIdx.x + threadIdx.x;
	if (i < ndata)
	{
		double this_chi[p];
		for (int j=0; j<p; j++) {
			this_chi[j] = chi[j * ndata + i];
		}
		logdens_pop[i] = logdensity_pop(this_chi, c_theta);
	}
}

/*
struct zsqr : public thrust::unary_function<double,double> {
    double mu, var;
    zsqr(double m, double v) : mu(m), var(v) {}

    __device__ __host__
    double operator()(double* chi) {
        double chi_cent = chi - mu;
        double logdens_pop = -0.5 * log(var) - 0.5 * chi_cent * chi_cent / var;
        return logdens_pop;
    }
};
*/

// Generate a data set and fill the chi, meas, and meas_sigma host vector with the values
void generate_data(int ndata, thrust::host_vector<double>& theta, thrust::host_vector<double>& chi, thrust::host_vector<double>& meas,
		thrust::host_vector<double>& meas_sigma) {

	double sigma_msmt[5] = {2.2e-4, 3.3e-4, 5.2e-4, 3.7e-4, 2.2e-4};  // Standard deviation for the measurement errors

	for (int j = 0; j < mfeat; ++j) {
		sigma_msmt[j] *= 1.0;
	}

	for (int i=0; i<ndata; i++) {
        // Loop over the data indices
		for (int j=0; j<p; j++) {
            // Loop over the chi indices to generate true value of the characteristics
			chi[i + ndata * j] = theta[j] + sqrt(exp(theta[p + j])) * snorm(rng);
        }
        for (int k=0; k<mfeat; k++) {
            // Loop over the feature indices to generate measurements, given the characteristics
            double bbody_numer = 2.0 * hplanck * nu[k] * nu[k] * nu[k] / (clight * clight);
            double temperature = exp(chi[ndata * 2 + i]); // Grab the temperature for this data point
            double bbody_denom = exp(hplanck * nu[k] / (kboltz * temperature)) - 1.0;
            double bbody = bbody_numer / bbody_denom;
            // Grab the SED normalization and power-law index
            double normalization = exp(chi[i]);
            double beta = chi[ndata + i];
            // Compute the model SED
            double nu_over_nu0 = nu[k] / nu0;
            double SED_ik = normalization * pow(nu_over_nu0, beta) * bbody;
            // Generate measured SEDs
			meas_sigma[k * ndata + i] = sigma_msmt[k];
            meas[k * ndata + i] = SED_ik + sigma_msmt[k] * snorm(rng);
		}
	}
}

// Initialize the chi and theta values. Right now this just initializes them to their true values. Also initialize the Cholesky factors
// for the chi and theta proposals. The elements of the chi and theta proposals are initially independent.
void initialize_parameters(int ndata, thrust::host_vector<double>& theta, thrust::host_vector<double>& chi,
		thrust::host_vector<double>& theta_cholfactor, thrust::host_vector<double>& chi_cholfactor, thrust::host_vector<double>& true_chi) {

	// For now, initialize chi values for drawing from their parent distribution.
	for (int i = 0; i < ndata; ++i) {
		for (int j = 0; j < p; ++j) {
			double chi_mean = theta[j];
			double chi_sig = sqrt(exp(theta[p + j]));
			chi[j * ndata + i] = chi_mean + chi_sig * snorm(rng);
		}
	}

	for (int j = 0; j < dim_theta; ++j) {
		theta[j] = theta[j] + theta[j] * 0.2 * snorm(rng);
	}

	// Initialize the cholesky factor for the theta proposals
    int diag_index = 0;
    for (int j=0; j<dim_theta; j++) {
    	// Initialize the covariance matrix of the theta proposal distribution to be 0.01 * identity(p)
        theta_cholfactor[diag_index] = 0.01;
        diag_index += j + 2;
    }

    // Initialize the cholesky factor for the chi proposals
    for (int i=0; i<ndata; i++) {
        int diag_index = 0;
        for (int j=0; j<p; j++) {
        	// Initialize the covariance matrix of the chi proposal distribution to be 0.01 * identity(p)
        	chi_cholfactor[i + ndata * diag_index] = 0.01;
        	diag_index += j + 2;
        }
    }
}

// Test the rank-1 cholesky update
void test_CholUpdateR1() {
	int dim_cf = 4;
	int size_cf = 10; //dim_cf * dim_cf - ((dim_cf - 1) * dim_cf) / 2
	double L[10] = {31.088, 7.759, 31.559, 7.632, 6.998, 29.739, 6.345, 7.080, 6.108, 29.015};
	// true values of the updated cholesky factor for A = L * transpose(L) +/- v * transpose(v)
	double Lup0[10] = {31.542, 4.726, 36.504, 6.898, 8.527, 29.748, 7.414, 3.262, 5.894, 30.270};
	double Ldown0[10] = {30.628, 10.884, 25.278, 8.389, 4.944, 29.718, 5.246, 13.219, 6.581, 26.017};
	double v[4] = {-5.331, 17.284, 3.691, -6.861};

	double Ldown[10];
	double Lup[10];
	double vup[4];
	double vdown[4];
	for (int i=0; i<size_cf; i++) {
		Lup[i] = L[i];
		Ldown[i] = L[i];
	}
	for (int i=0; i<dim_cf; i++) {
		vup[i] = v[i];
		vdown[i] = v[i];
	}

	// test the update first
	CholUpdateR1(Lup, vup, dim_cf, false);
	CholUpdateR1(Ldown, vdown, dim_cf, true);
	double fracdiff_up = 0.0;
	double fracdiff_down = 0.0;
	for (int i=0; i<size_cf; i++) {
		std::cout << "Lup[i]: " << Lup[i] << ", Ldown[i]: " << Ldown[i] << std::endl;
		double fracdiff = abs(Lup[i] - Lup0[i]) / Lup0[i];
		fracdiff_up = max(fracdiff_up, fracdiff);
		fracdiff = abs(Ldown[i] - Ldown0[i]) / Ldown0[i];
		fracdiff_down = max(fracdiff_down, fracdiff);
	}
	std::cout << "test_CholUpdateR1:" << std::endl;
	std::cout << "Maximum fractional difference for rank-1 update: " << fracdiff_up << std::endl;
	std::cout << "Maximum fractional difference for rank-1 downdat: " << fracdiff_down << std::endl;
}

int main(void)
{
    /* TODO: Add command line options */

    /*
     Population level parameters: theta, where theta parameterizes the distribution of chi.
     For this application, chi = (log10(norm), beta, log10(temp)), where the model SED is a
     modified black body:

        SED(nu) = norm * (nu / nu0)^beta * Bbody(nu,temp).

     Here, Bbody(nu,temp) is the Planck function. For simplicity, right now I just assume that
     log10(norm), beta, and log10(temp) are independently drawn from different Gaussian distributions,
     although we should later include correlations among the parameters. So, for now, theta is just the
     collection of the mean and variances of the normal distributions.
    */

    int ndata = 1000; // # of data points

    double mu_norm = 8.5 * log(10.0);  // Average value of the natural logarithm of the SED normalization
    double sig_norm = 0.5 * log(10.0);  // Standard deviation in the SED normalization
    double mu_beta = 2.0;  // Average value of the SED power-law index
    double sig_beta = 0.2;  // Standard deviation in the SED power-law index
    double mu_temp = log(15.0);  // Average value of the logarithm of the dust temperature
    double sig_temp = 0.2;  // Standard deviation in the logarithm of the SED temperature

    thrust::host_vector<double> h_theta(dim_theta);  // Allocate memory on host
    h_theta[0] = mu_norm;
    h_theta[1] = mu_beta;
