#include <iostream>
#include <boost/timer.hpp>
#include <boost/random/linear_congruential.hpp>
#include <boost/random/uniform_int.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/random/variate_generator.hpp>
#include <boost/generator_iterator.hpp>
#include <boost/program_options.hpp>

// This is a basic application program which try
// to show the performance difference between an
// NRC storage and a classic one
//
// The computation result must be the same between
// nsr and classic access.
//
// Here is the program option
// Allowed options:
//   --help                produce help message
//   --width arg (=1024)   width of the array
//   --height arg (=768)   height of the array
//   --loop arg (=1000)    number of run loop
//   --seed arg (=42)      seed for random generator
//
// To compile this application, you'll need
// the boost program options library. (This is the
// libboost-program-options-dev under debian)
//
// compile this program with :
// g++ array.cpp -W -Wall -O3 -lboost_program_options

// a random generator using a static seed
// to reproduce result
typedef boost::minstd_rand base_generator_type ;
base_generator_type generator(42u) ;

// 2d array allocation using NRC storage
template <typename T>
T** alloc_array(int h, int w)
{
    typedef T* ptr_type ;
    typedef ptr_type* return_type ;

    return_type m = new ptr_type[h] ;
    if (!m) return 0 ;

    m[0] = new T[w*h] ;
    if (!m[0]) { delete[] m; return 0 ; }

    for (int i=1 ; i < h ; ++i) m[i] = m[i-1] + w ;
    return m ;
}

// helper function to release memory
template <typename T>
void release_array(T** array)
{
    delete[] array[0] ;
    delete[] array ;
}

// fill array with value between bv & uv
template <typename T>
void fill_array(T** array, int h, int w, T bv=-10, T uv=10)
{
    boost::uniform_real<T> uni_dist(bv, uv) ;
    boost::variate_generator<base_generator_type&, boost::uniform_real<T> > uni(generator, uni_dist) ;

    for (int i = 0 ; i < h ; ++i)
        for (int j = 0 ; j < w ; ++j) {
            array[i][j] = uni() ;
        }
}

// function using to simulate work and array access
// if ncr_access is true, then ncr access is used,
// otherwise, "classic" unidimensionnal access is used.
template <typename T>
T simulate_work(T** array, int h, int w, int loop, bool ncr_access=true)
{
    T res = 0 ;

    // we generate some amount of work here
    if (ncr_access) { // using ncr access
        for (int i = 0 ; i < loop ; ++i) {

            res = 0 ;
            // we loop over the whole array
            for (int ti = 0 ; ti < h ; ++ti)
            {
                for (int tj = 0 ; tj < w ; ++tj)
                {
                    res += array[ti][tj] ;
                }
            }

        }
        return res ;
    }
    else { // use naive access
        for (int i = 0 ; i < loop ; ++i) {

            // begin will help us to mimic the naive access
            // using the NRC allocated array
            // It's just like begin = new T[w*h]
            T* begin = array[0] ;

            res = 0 ;
            // we loop over the whole array
            for (int ti = 0 ; ti < h ; ++ti)
            {
                for (int tj = 0 ; tj < w ; ++tj)
                {
                  //classic element access v1
                  // the least effective as it requires arithmetic
                  //res += begin[w*ti + tj];

                  //classic element access v2
                  // seems better as it only increments a pointer
                  res += *begin++;
                }
            }

        }
        return res ;
     }
}


int main(int argc, char* argv[])
{
    namespace po = boost::program_options ;

    // we will use a 2d float array
    typedef float value_type ;
    typedef value_type** Array ;

    int width, height, bench_loop;
    unsigned int seed ;

    // options parsing
    po::options_description desc("Allowed options") ;
    desc.add_options()
        ("help", "produce help message")
        ("width",  po::value<int>(&width)->default_value(512), "width of the array")
        ("height", po::value<int>(&height)->default_value(512), "height of the array")
        ("loop",   po::value<int>(&bench_loop)->default_value(10000), "number of run loop")
        ("seed",   po::value<unsigned int>(&seed)->default_value(42u), "seed for random generator") ;

    po::variables_map vm ;

    try {
        po::store(po::parse_command_line(argc, argv, desc), vm) ;
        po::notify(vm) ;
    }
    catch (po::unknown_option& e) {
        std::cout << desc << '\n' ;
        return 1 ;
    }

    if (vm.count("help")) {
        std::cout << desc << '\n' ;
        return 1;
    }

    // array allocation
    Array ncr_array = alloc_array<value_type>(height, width) ;
    if (!ncr_array) {
        std::cerr << "Error during array allocation\n" ;
        return -1 ;
    }

    // we fill the array with some random value
    generator.seed(seed) ;
    fill_array(ncr_array, height, width) ;

    value_type res_nsr, res_classic ;
    double tt,ttc;

    boost::timer timer ;
    {
        res_nsr = simulate_work(ncr_array, height, width, bench_loop, true) ;
        tt = timer.elapsed() ;
        std::cout.precision(10);
        std::cout << "NRC access: " << tt << " s\n" ;
    }

    // cleanup
    release_array(ncr_array) ;

    // reallocation
    ncr_array = alloc_array<value_type>(height, width) ;
    if (!ncr_array) {
        std::cerr << "Error during array allocation\n" ;
        return -1 ;
    }

    // again, we will the array with some random value
    generator.seed(seed) ;
    fill_array(ncr_array, height, width) ;

    {
        timer.restart() ;
        res_classic = simulate_work(ncr_array, height, width, bench_loop, false) ;
        ttc = timer.elapsed() ;
        std::cout.precision(10);
        std::cout << "Naive access: " << ttc << " s\n" ;
    }

    std::cout.precision(3);
    std::cout << "Performance ration: " << ((tt-ttc)/tt)*100 << "%\n" ;
    if (res_classic == res_nsr) {
        std::cout << "Computation is ok\n" ;
    }
    else {
        std::cout << "Result between nsr access and classic access are different. Not normal.\n" ;
    }

    release_array(ncr_array) ;
}