// This provides a way of constructing sparse matrices a la MTL via a hash_map
#ifndef HASH_MATRIX_H
#define HASH_MATRIX_H

#include <ext/hash_map>  // not part of STL
#include <vector>
#include <algorithm>
#include <mtl/mtl.h>

namespace hash_mat {

  template <class int_type>
  class myHash {
  public:
    size_t operator()(const std::pair<int_type,int_type> &p) const
    {  return p.first + 71437*p.second; }
  };

  template <class SpMat, 
            // class int_type = typename SpMat::int_type, 
            class int_type = int,
            class entry_loc = typename std::pair<int_type,int_type>,
            class value_type = typename SpMat::value_type>
  class hash_matrix : public __gnu_cxx::hash_map<entry_loc, value_type, myHash<int_type> > {
  public:
    typedef typename __gnu_cxx::hash_map<entry_loc, value_type, myHash<int_type> > hash_type;
    typedef typename hash_type::iterator ht_iterator;

    hash_matrix(int_type init_size =100) : hash_type(init_size)
    { }

    value_type &operator()(int_type i, int_type j)
    {
      return (this->operator[](entry_loc(i,j)));
    }

    // typedef typename hash_type::value_type triple;
    typedef typename std::pair<entry_loc, value_type> triple;
    typedef typename std::vector<triple> triple_list;
    class LexicoCompLocs {
    public:
      int operator()(const triple &t1, const triple &t2)
      {
        const entry_loc &e1 = t1.first, &e2 = t2.first;
        if ( e1.first < e2.first )
          return 1;
        else if ( e1.first == e2.first )
          return ( e1.second < e2.second );
        else
          return 0;
      }
    } lexico_cmp;

    SpMat &convert(SpMat &out)
    {
      triple_list triples;
      int idx = 0;

      triples.resize(this->size());
      checkpt();
      for ( ht_iterator i = this->begin(); i != this->end(); ++i )
        {
          assert(idx < this->size());
          triples[idx++] = *i;
        }
      checkpt();

      // Now sort triples
      std::sort(triples.begin(), triples.end(), this->lexico_cmp);

      checkpt();

      // std::cout << "length of triples vector = " << triples.size() << std::endl;
      for ( int i = 0; i+1 < triples.size(); ++i )
        {
          // std::cout << "checking triples: i = " << i << std::endl;
          if ( ! this->lexico_cmp(triples.at(i), triples.at(i+1)) )
            {
              checkpt();
              std::cout << "List not sorted" << std::endl;
              break;
            }
        }

      checkpt();

      // Now create sparse MTL matrix from external arrays
      // generated from the triples vector
      // First create the external arrays
      int nnz = triples.size();
      value_type *val_list   = new value_type[nnz];
      int_type   *col_list   = new int_type[nnz];
      int_type   *start_list = new int_type[out.nrows()+1];

      checkpt();

      // Fill val_list & col_list from triples
      int row_num;
      for ( row_num = 0; row_num <= out.nrows(); ++row_num )
        start_list[row_num] = -1;
      row_num = 0;
      start_list[row_num] = 0;
      checkpt();
      for ( idx = 0; idx < nnz; ++idx )
        {
          triple &t = triples[idx];
          val_list[idx] = t.second;
          col_list[idx] = t.first.second;
          if ( t.first.first != row_num )
            {
              row_num = t.first.first;
              start_list[row_num] = idx;
            }
        }
      checkpt();
      // start_list must end with nnz, and we must 
      // fill in any missing values
      start_list[out.nrows()] = nnz;
      for ( row_num = out.nrows(); row_num > 0; --row_num )
        if ( start_list[row_num-1] == -1 )
          start_list[row_num-1] = start_list[row_num];
      checkpt();

//       for ( row_num = 0; row_num < out.nrows(); ++row_num )
//         {
//           std::cout << "Row " << row_num << "  starts at index " << start_list[row_num] << std::endl;
//           for ( idx = start_list[row_num]; idx < start_list[row_num+1]; ++idx )
//             std::cout << col_list[idx] << ":" << val_list[idx] << "  ";
//           std::cout << std::endl;
//         }
//       std::cout << "col_list and val_list array values:" << std::endl;
//       for ( idx = 0; idx < nnz; ++idx )
//         {
//           std::cout << "idx = " << idx << ", col_list[idx] = " << col_list[idx] <<
//             ", val_list[idx] = " << val_list[idx] << std::endl;
//         }

      // construct the matrix using external arrays
      typedef typename mtl::matrix<value_type, mtl::rectangle<>, mtl::compressed<int, mtl::external, mtl::index_from_zero> , mtl::row_major>::type ExtSpMat;
      ExtSpMat ext(out.nrows(), out.ncols(), nnz, val_list, start_list, col_list);
      checkpt();
      // and copy this into out
      copy(ext, out);
      checkpt();
      // clean up
      delete[] val_list;
      delete[] col_list;
      delete[] start_list;
      checkpt();

      return out;
    }
      

  };

} // namespace hash_mat


#endif // HASH_MATRIX_H