#include <iostream>
#include <string>
#include <ctime>

// #define USE_GENERALIZED_VECTOR_OF_VECTOR

// Boost
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/vector_of_vector.hpp>
#include <boost/numeric/ublas/vector_sparse.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/fwd.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/io.hpp>




using namespace boost::numeric::ublas;

typedef generalized_vector_of_vector< double, row_major, vector<compressed_vector<double> > > XType;

/**
 * Compute v * M
 */
// vector<double> ourProd(const vector<double>& v, const XType& M){
// 	vector<double> r(M.size2());
// 	const vector<compressed_vector<double> >& data = M.data();
// 	for(size_t i=0; i < M.size1(); ++i){
// 		r += v(i) * data(i);
// 	}
// 	return r;
// }

template <class RType, class MType, class VType>
RType ourProd(const VType& v, const MType& M){
	RType r(M.size2());
	
	std::cout << "M: " << M.size1() << "x" << M.size2() << std::endl;
	r = v(0) * matrix_row<const MType>(M,0);
	
	
	for(size_t i=1; i < M.size1(); ++i){
		matrix_row<const MType> mr(M,i);
		r += v(i) * mr;
	}
	return r;
}




template <class M> 
void fill_ublas(const size_t row, const size_t col, const char* data, M& m){
  for(size_t i = 0; i<row; ++i) {
    for(size_t j = 0; j<col; ++j) {
			char value = data[i*col+j];
    	if(value > 0){
        m(i,j) = value;
			}
    } // col
  }// row
}


template <class MType>
void test_ublas(const size_t row, const size_t col, const char* data, const unsigned int maxloop){
	MType X(row,col);
  fill_ublas(row, col, data, X);
  vector<double> w(X.size2());
	for(size_t i=0; i<w.size(); ++i) w(i) = 1.0;
	
  
  vector<double> result;
	
  const clock_t rightMultStartTime = clock();
  for(unsigned int loop=0; loop<maxloop; loop++) {
		result = prod(X, w);
  }
  const clock_t rightMultEndTime  = clock();
  const double rm_1norm = norm_1(result);
  
  vector<double> l(row);
  for(size_t i=0; i<l.size(); ++i){
    l[i] = 1.0;
  }
  
  vector<double> ourGrad;
  
  const clock_t leftMultStartTime = clock();
  for(unsigned int i=0; i<maxloop; i++) {
  		ourGrad = ourProd<vector<double>, MType, vector<double> >(l,X);
  }
	const clock_t leftMultEndTime = clock();
	const double our_l1norm = norm_1(ourGrad);
	assert(l.size() == row);
	assert(ourGrad.size() == col);
	assert(X.size1() == row);
	assert(X.size2() == col);

  

	vector<double> uBLASGrad;
  const clock_t uBlasLeftMultStart = clock();
  for(unsigned int i=0; i<maxloop; i++) {
    uBLASGrad = prod(l, X);
 	}

	const clock_t uBlasLeftMultStop = clock();
  const double ublas_l1norm = norm_1(uBLASGrad);

  const double ublasleftTime  = ((double)(uBlasLeftMultStop - uBlasLeftMultStart))/CLOCKS_PER_SEC;
  const double leftTime       = ((double)(leftMultEndTime - leftMultStartTime))/CLOCKS_PER_SEC;
  const double rightTime      = ((double)(rightMultEndTime - rightMultStartTime))/CLOCKS_PER_SEC;
  
  // display time elapsed
  printf("rightmult time            : %.4lf\n",rightTime);
  printf("our ublas leftmult time   : %.4lf\n",leftTime);
  printf("ublas leftmult time       : %.4lf\n",ublasleftTime);
  printf("rightmult 1-norm          : %.8lf\n",rm_1norm);
  printf("ublas leftmult 1-norm     : %.8lf\n",ublas_l1norm);
  printf("our ublas leftmult 1-norm : %.8lf\n",our_l1norm);
}



// zero-based indexing sparse vector
class sparsevec {
public:

  int len;
  long *fidx;
  double *fval;
  
  sparsevec():len(0),fidx(0),fval(0) {}

  sparsevec(int _len) {
    len = _len;
    fidx = new long[len];
    fval = new double[len];

    memset(fidx,0,sizeof(long)*len);
    memset(fval,0,sizeof(double)*len);
  }  

  sparsevec(int actual_len, const sparsevec& sp)
  {
    len = actual_len;
    fidx = new long[len];
    fval = new double[len];
    
    memcpy(fidx, sp.fidx, sizeof(long)*len);
    memcpy(fval, sp.fval, sizeof(double)*len);
  }

 
  ~sparsevec() {
    if(fidx) delete [] fidx;
    if(fval) delete [] fval;
    len = 0;
  }
};

class densevec {
public:
  int len;
  double *array;

  densevec():len(0),array(0) {}

  densevec(int _len) 
  {
    len = _len;
    array = new double[len];
    memset(array,0,sizeof(double)*len);
  }

  

  ~densevec() 
  {
    if(array) delete [] array;
    len = 0;
  }

  
  // dot product
  double dot(const sparsevec *svec) 
  {
    double res = 0.0;
    
    if(svec->len == 0) return res;


    assert(svec->len > 0);
    
    for(register int i=0; i<svec->len; i++)
      res += array[svec->fidx[i]] * svec->fval[i];
    
    return res;
  }


  // add with scaled sparsevec
  void add_and_scale(const sparsevec *svec, const double& weight) 
  {
    assert(len > 0);
    assert(svec->len > 0);
    assert(svec->fidx[svec->len-1] <= len);
    
    for(register int i=0; i<svec->len; i++)
      array[svec->fidx[i]] += weight * svec->fval[i];
  }


  // compute l1 norm
  double l1norm() 
  {
    double res = 0;
    for(int i=0; i<len; i++)
      res += fabs(array[i]);
    return res;
  }

	void zero(){
		memset(array, 0, sizeof(double)*len);
	}
};


void fill_fast(const size_t row, const size_t col, const char* data, std::vector<sparsevec*>& m)
{
  int actual_len = 0;
  sparsevec sp(col);
  
  for(size_t i=0; i<row; i++) {
    actual_len = 0;
    for(size_t j=0; j<col; j++) {
      if(fabs(data[i*col+j]) > 0) {
        sp.fidx[actual_len] = j;
        sp.fval[actual_len] = (double) data[i*col+j];
        actual_len++;
      }
    }
    m.push_back(new sparsevec(actual_len, sp));
  }
}


void test_fast(const size_t row, const size_t col, const char* data, const unsigned int maxloop)
{
  std::vector<sparsevec*> X;
  fill_fast(row, col, data, X);
  
  densevec w(col);
  for(size_t i=0; i<col; ++i){
    w.array[i] = 1.0;
  }
  
  vector<double> result(row);
  
  clock_t rightMultStartTime = clock();
  for(unsigned int loop=0; loop<maxloop; loop++) {
    for(size_t i=0; i<row; i++) {
      result[i] = w.dot(X[i]);
    }    
  }
  clock_t rightMultEndTime  = clock();
  double rm_1norm = 0;
  for(size_t i=0; i<row; i++) rm_1norm += fabs(result[i]);
  
  vector<double> l(row);
  for(size_t i=0; i<row; ++i){
    l[i] = 1.0;
  }
  
  densevec grad(col);
  
  clock_t leftMultStartTime = clock();
  for(unsigned int i=0; i<maxloop; i++) {  
		grad.zero();
		for(size_t i=0; i<row; i++) {
      grad.add_and_scale(X[i],l[i]);
    }
  }
  clock_t leftMultEndTime = clock();
  double lm_1norm = grad.l1norm();

  double leftTime  = ((double)(leftMultEndTime - leftMultStartTime))/CLOCKS_PER_SEC;
  double rightTime = ((double)(rightMultEndTime - rightMultStartTime))/CLOCKS_PER_SEC;
  
  // display time elapsed
  printf("rightmult time            : %.4lf\n",rightTime);
  printf("leftmult time             : %.4lf\n",leftTime);
  printf("rightmult 1-norm          : %.8lf\n",rm_1norm);
  printf("leftmult 1-norm           : %.8lf\n",lm_1norm);
  
}

int main(int argc, char** argv)
{	
  int row     = atoi(argv[1]);
  int col     = atoi(argv[2]);
  int maxLoop = atoi(argv[3]);
  double nnz_prob = atof(argv[4]);
  
  // generate random sparse matrix
  int nnz = 0;
  char *data = new char[row*col];
  for(int i=0; i<row; i++) {
    for(int j=0; j<col; j++) {
      if(drand48() < nnz_prob) {
        data[i*col + j] = 1;
        nnz++;
      }
    }
  }
  printf("-----------------------------------------------------------\n");

  printf("nnz in data: %d\n",nnz);

  // printf("\n\n--------------------------\n");
  // printf("ublas dense\n");
  // printf("--------------------------\n");
  // test_ublas<matrix<double> >(row, col, data, maxLoop);

  printf("\n\n--------------------------\n");
  printf("ublas compressed_matrix\n");
  printf("--------------------------\n");
  test_ublas<compressed_matrix<double> >(row, col, data, maxLoop);

  printf("\n\n--------------------------\n");
  printf("fast\n");
  printf("--------------------------\n");
  test_fast(row, col, data, maxLoop);
  
  printf("-----------------------------------------------------------\n\n\n");

  return 0;
}