 #include <stdio.h>
 #include <fstream>
 #include <stdlib.h>
+#include <iostream>
 #include <math.h>
 #include <vector>
+
+// Boost includes
+#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/normal_distribution.hpp>
+#include <boost/random/uniform_real_distribution.hpp>
 
 //#include <thrust/sort.h>
 //#include <thrust/copy.h>

        for (int j=0; j<dim_theta; j++) {
            snorm_deviate[j] = snorm(rng);
        }
        
            // transform unit proposal so that is has a multivariate normal distribution
        thrust::host_vector<double> scaled_proposal(dim_theta);
        thrust::fill(scaled_proposal.begin(), scaled_proposal.end(), 0.0);
        int cholfact_index = 0;
        for (int j=0; j<dim_theta; j++) {
            for (int k=0; k<(j+1); k++) {
                // CHECK THIS: SHOULDN'T THIS BE SNORM_DEVIATE[K]?????
                scaled_proposal[j] += h_theta_cholfact[cholfact_index] * snorm_deviate[j];
                cholfact_index++;
            }
            scaled_proposal[j] = 0.0;
            proposed_theta[j] = h_theta[j] + scaled_proposal[j];
        }
        
        CUDA_CALL(cudaDeviceSynchronize());
        h_chi = d_chi;
        h_logdens = d_logdens;
        
        h_mha_ratios = d_mha_ratios;
        double min_mha_ratio = 1e10;
        double max_mha_ratio = -1.0;
        for (int idata=0; idata<n; idata++) {
            min_mha_ratio = std::min(min_mha_ratio, h_mha_ratios[idata]);
            max_mha_ratio = std::max(max_mha_ratio, h_mha_ratios[idata]);
        }
        std::cout << std::endl;
        std::cout << "Iteration: " << current_iter << std::endl;
        std::cout << "Minimum, Maximum MHA ratio: " << min_mha_ratio << ", " << max_mha_ratio << std::endl;
        
        /*
        std::cout << "Proposed theta: ";
        for (int j=0; j<dim_theta; ++j) {
            std::cout << proposed_theta[j] << " ";
        }
        std::cout << std::endl;
        std::cout << "Current theta: ";
        for (int j=0; j<dim_theta; ++j) {
            std::cout << h_theta[j] << " ";
        }
        std::cout << std::endl;
         */
        // Compute Metropolis ratio
        // right now we loop over the elements of chi because we assume that they are statistically independent, i.e.,
        // the theta array does not contain correlations among the element of chi. this is unrealistic, and we will need
        // to come up with a better way to do the transform + reduction: maybe use CUBLAS?
        thrust::device_vector<double>::iterator chi_iter_begin = d_chi.begin();
        
        double logdens_pop = 0.0;
        for (int j=0; j<dim_theta/2; j++) {
            double proposed_mu_j = proposed_theta[j];
            double proposed_var_j = exp(proposed_theta[j + p]);
            /*
            std::cout << "h_chi (begin,end): " << h_chi[j * n] << " " << h_chi[(j+1)*n-1] << std::endl;
            std::cout << "d_chi (begin,end): " << d_chi[j * n] << " " << d_chi[(j+1)*n-1] << std::endl;
            std::cout << "chi_iter (begin,end): " << *chi_iter_begin << " " << *(chi_iter_begin+n-1) << std::endl;
             */
            for (int idata=0; idata<n; idata++) {
                std::cout << "h_chi: " << h_chi[idata + j*n] << ", iterator: " << *(chi_iter_begin+idata) << std::endl;
            }
            // transform and reduction is over d_chi[j * n : (j+1) * n - 1]
            zsqr zsqrj(proposed_mu_j, proposed_var_j);
            logdens_pop = thrust::transform_reduce(chi_iter_begin, chi_iter_begin+n,
                                                zsqrj, logdens_pop, thrust::plus<double>());
            thrust::advance(chi_iter_begin, n);
        }
        std::cout << "logdens_pop from transform_reduce: " << logdens_pop << std::endl;
        logdens_pop = 0.0;
        zsqr zsqr0(h_theta[0], exp(h_theta[3]));
        zsqr zsqr1(h_theta[1], exp(h_theta[4]));
        zsqr zsqr2(h_theta[2], exp(h_theta[5]));
        for (int j=0; j<n; j++) {
            double chi_cent = h_chi[j] - h_theta[0];
            logdens_pop += zsqr0(h_chi[j]);
            //logdens_pop += -0.5 * h_theta[3] - 0.5 * chi_cent * chi_cent / exp(h_theta[3]);
            chi_cent = h_chi[j+n] - h_theta[1];
            logdens_pop += zsqr1(h_chi[j+n]);
            //logdens_pop += -0.5 * h_theta[4] - 0.5 * chi_cent * chi_cent / exp(h_theta[4]);
            chi_cent = h_chi[j+2*n] - h_theta[2];
            logdens_pop += zsqr2(h_chi[j+2*n]);
            //logdens_pop += -0.5 * h_theta[5] - 0.5 * chi_cent * chi_cent / exp(h_theta[5]);
        }
        std::cout << "logdens_pop: " << logdens_pop << std::endl;
        std::cout << "h_theta: ";
        for (int j=0; j<dim_theta; j++) {
            std::cout << " " << h_theta[j];
        }
        std::cout << std::endl;
        
        double logdens_pop2 = 0.0;
        for (int j=0; j<n; j++) {
            double logdens_pop_j = 0.0;
            double chi_cent = h_chi[j] - h_theta[0];
            logdens_pop_j += -0.5 * h_theta[3] - 0.5 * chi_cent * chi_cent / exp(h_theta[3]);
            chi_cent = h_chi[j+n] - h_theta[1];
            logdens_pop_j += -0.5 * h_theta[4] - 0.5 * chi_cent * chi_cent / exp(h_theta[4]);
            chi_cent = h_chi[j+2*n] - h_theta[2];
            logdens_pop_j += -0.5 * h_theta[5] - 0.5 * chi_cent * chi_cent / exp(h_theta[5]);
            double logdens_pop_j2;
            double chi_j[3];
            chi_j[0] = h_chi[j];
            chi_j[1] = h_chi[j+n];
            chi_j[2] = h_chi[j+2*n];
            std::cout << std::endl;
            std::cout << "iteration: " << j << std::endl;
            logdens_pop_j2 = h_logdensity_pop(chi_j, &h_theta[0]);
            std::cout << "chi: " << h_chi[j] << " " << h_chi[j+n] << " " << h_chi[j+2*n] << std::endl;
            std::cout << chi_j[0] << " " << chi_j[1] << " " << chi_j[2] << std::endl;
            std::cout << "logdensity: " << logdens_pop_j << " " << logdens_pop_j2 << std::endl;
            logdens_pop2 += logdens_pop_j;
        }
        
        
        double logdens_old = thrust::reduce(d_logdens.begin(), d_logdens.end());
        std::cout << "logdens_pop: " << logdens_pop << " logdens_old: " << logdens_old <<
        " logdens_pop2: " << logdens_pop2 << std::endl;
        double lograt = logdens_pop - logdens_old;
        lograt = std::min(lograt, 0.0);
        double ratio = exp(lograt);
        double unif = uniform(rng);
        
        //std::cout << "MHA ratio for theta: " << ratio << std::endl;
        
        if (unif < ratio) {
            // Accept the proposed theta
            h_theta = proposed_theta;
            // grab the pointer in case thrust changed the memory location on us
            double* p_theta = thrust::raw_pointer_cast(&h_theta[0]);
            // copy the new theta to constant memory on the device
            CUDA_CALL(cudaMemcpyToSymbol(c_theta, p_theta, h_theta.size() * sizeof(*p_theta)));
            naccept_theta++;
        }
        
        // finally, adapt the proposal covariance of theta by updating its cholesky factor
        double unit_norm = 0.0;
        for (int j=0; j<dim_theta; j++) {
            unit_norm += snorm_deviate[j] * snorm_deviate[j];
        }
        
        unit_norm = sqrt(unit_norm);
        double decay_sequence = 1.0 / pow(current_iter, decay_rate);
        double scaled_coef = sqrt(decay_sequence * fabs(ratio - theta_target_rate)) / unit_norm;
        for (int j=0; j<dim_theta; j++) {
            scaled_proposal[j] *= scaled_coef;
        }
        
        bool downdate = (ratio < theta_target_rate);
        double* p_theta_cholfact = thrust::raw_pointer_cast(&h_theta_cholfact[0]);
        double* p_scaled_proposal = thrust::raw_pointer_cast(&scaled_proposal[0]);
        CholUpdateR1(p_theta_cholfact, p_scaled_proposal, dim_theta, downdate);
        
        // Save the theta values
        for (int j=0; j<dim_theta; j++) {
            thetafile << h_theta[j] << " ";
        }
        thetafile << std::endl;
    }
    std::cout << "Number of accepted thetas: " << naccept_theta << std::endl;
    thetafile.close();
    
    cudaFree(devStates);  // Free up the memory on the GPU from the RNG states
    
	return 0;
}

/*
 * dust_seds.cu
 *
 *  Created on: Jun 2, 2013
 *      Author: brandonkelly
 */

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