#ifndef _utl_SparseVector_h_
#define _utl_SparseVector_h_

#include <map>
#include <numeric>

#include <boost/lambda/lambda.hpp>


namespace utl {

  /**
    \class SparseVector SparseVector.h "utl/SparseVector.h"
    \brief Sparse container class for vectorial data

    Importatant notes:
    <ul>
    <li>
      the operator[] access is not creating entries. It should be used for
      reading of the whole vector.
    <li>
      the operator() access creates the corresponding entry on the fly and
      initializes it to GetInitializer() <b>before</b> the reference to the
      element is returned
    </ul>

    Please refer to the testSparseVector.cc and test SparseMatrixVectorOps.cc
    files for detailed usage cases.

    \author Darko Veberic
    \version $Id$
    \date 3 Jun 2008
    \ingroup math
  */

  template<typename T = double>
  class SparseVector {
  public:
    typedef int IndexType;

  private:
    typedef std::map<IndexType, T> MapType;
    typedef typename MapType::value_type MapValueType;
    typedef typename MapType::iterator MapIterator;
    typedef typename MapType::const_iterator MapConstIterator;

  public:
    template<class IteratorType, typename ValueType>
    class IteratorTransformer {
    public:
      IteratorTransformer(const IteratorType& it)
        : fIterator(it) { }

      IndexType GetIndex() const { return fIterator->first; }

      ValueType& operator()() const { return fIterator->second; }

      ValueType& operator*() const { return fIterator->second; }

      IteratorTransformer& operator++() { ++fIterator; return *this; }

      bool operator==(const IteratorTransformer& it) const { return fIterator == it.fIterator; }

      bool operator!=(const IteratorTransformer& it) const { return fIterator != it.fIterator; }

    private:
      IteratorType fIterator;
    };

    SparseVector(const T& init = T())
      : fInitializer(init) { }

    template<typename U>
    explicit
    SparseVector(const SparseVector<U>& v)
      : fInitializer(v.GetInitializer())
    {
      for (typename SparseVector<U>::ConstIterator it = v.SparseBegin();
           it != v.SparseEnd(); ++it)
        (*this)(it.GetIndex()) = it();
    }

    unsigned int GetNonEmptySize() const { return fVector.size(); }

    const T& GetInitializer() const { return fInitializer; }

    void Clear() { fVector.clear(); }

    /// noncreating element access
    const T&
    operator[](const IndexType index)
      const
    {
      const MapConstIterator it = fVector.find(index);
      return (it == fVector.end()) ? fInitializer : it->second;
    }

    /// creating element access
    T&
    operator()(const IndexType index)
    {
      return fVector.insert(MapValueType(index, fInitializer)).first->second;
    }

    const T&
    operator()(const IndexType index)
      const
    {
      const MapConstIterator it = fVector.find(index);
      return (it == fVector.end()) ? fInitializer : it->second;
    }

    bool IsEmpty() const { return fVector.empty(); }

    bool
    IsEmpty(const IndexType index)
      const
    {
      const MapConstIterator it = fVector.find(index);
      return it == fVector.end();
    }

    /// remove "empty" elements (that are equal to the initializer)
    unsigned int
    Reclaim()
    {
      unsigned int erased = 0;
      for (MapIterator it = fVector.begin();
           it != fVector.end(); )
        if (it->second != fInitializer)
          ++it;
        else {
          Erase(it++);
          ++erased;
        }
      return erased;
    }

    bool
    operator==(const SparseVector<T>& v)
      const
    {
      return fInitializer == v.fInitializer &&
             fVector == v.fVector;
    }

    bool operator!=(const SparseVector<T>& v) const { return !operator==(v); }

    template<typename U>
    bool
    operator==(const SparseVector<U>& v)
      const
    {
      if (fInitializer != v.GetInitializer())
        return false;
      ConstIterator it1 = SparseBegin();
      typename SparseVector<U>::ConstIterator it2 = v.SparseBegin();
      while (it1 != SparseEnd() && it2 != v.SparseEnd()) {
        if (it1.GetIndex() != it2.GetIndex() ||
            it1() != it2())
          return false;
        ++it1;
        ++it2;
      }

      if (it1 != SparseEnd() || it2 != v.SparseEnd())
        return false;

      return true;
    }

    template<typename U>
    bool operator!=(const SparseVector<U>& v) const
    { return !operator==(v); }

    typedef IteratorTransformer<MapIterator, T> Iterator;
    typedef IteratorTransformer<MapConstIterator, const T> ConstIterator;

    Iterator SparseBegin() { return Iterator(fVector.begin()); }

    ConstIterator SparseBegin() const { return ConstIterator(fVector.begin()); }

    Iterator SparseEnd() { return Iterator(fVector.end()); }

    ConstIterator SparseEnd() const { return ConstIterator(fVector.end()); }

    Iterator Find(const IndexType index)
    { return Iterator(fVector.find(index)); }

    ConstIterator Find(const IndexType index) const
    { return ConstIterator(fVector.find(index)); }

    /// multiplication with scalar
    template<typename U>
    SparseVector&
    operator*=(const U& value)
    {
      for (MapIterator it = fVector.begin(); it != fVector.end(); ++it)
        it->second *= value;
      return *this;
    }

    /// vector addition
    template<typename U>
    SparseVector&
    operator+=(const SparseVector<U>& m)
    {
      for (typename SparseVector<U>::ConstIterator it = m.SparseBegin();
           it != m.SparseEnd(); ++it)
        (*this)(it.GetIndex()) += it();
      return *this;
    }

    /// vector subtraction
    template<typename U>
    SparseVector&
    operator-=(const SparseVector<U>& m)
    {
      for (typename SparseVector<U>::ConstIterator it = m.SparseBegin();
           it != m.SparseEnd(); ++it)
        (*this)(it.GetIndex()) -= it();
      return *this;
    }

  private:

    void Erase(const MapIterator it) { fVector.erase(it); }

    const T fInitializer;
    MapType fVector;
  };


  /// dot product (scalar product of vectors)
  template<typename T, typename U>
  typename boost::lambda::return_type_2<
    boost::lambda::arithmetic_action<boost::lambda::multiply_action>,
    SparseVector<T>,
    SparseVector<U>
  >::type
  operator*(const SparseVector<T>& a, const SparseVector<U>& b)
  {
    typedef typename boost::lambda::return_type_2<
      boost::lambda::arithmetic_action<boost::lambda::multiply_action>,
      SparseVector<T>,
      SparseVector<U>
    >::type Result;

    Result r = Result();

    typename SparseVector<T>::ConstIterator aIt = a.SparseBegin();
    typename SparseVector<U>::ConstIterator bIt = b.SparseBegin();
    if (aIt == a.SparseEnd() || bIt == b.SparseEnd())
      return r;
    typename SparseVector<T>::IndexType ai;
    typename SparseVector<U>::IndexType bi;
    do {
      ai = aIt.GetIndex();
      bi = bIt.GetIndex();
      if (ai < bi) {
        ++aIt;
        continue;
      }
      if (bi < ai) {
        ++bIt;
        continue;
      }
      r += aIt() * bIt();
      ++aIt;
      ++bIt;
    } while (aIt != a.SparseEnd() && bIt != b.SparseEnd());

    return r;
  }


  // output stuff

  template<class Stream, typename T>
  Stream&
  Output(Stream& s, const SparseVector<T>& v,
         const unsigned int size, const char separator = ' ')
  {
    if (!size)
      return s;

    if (v.IsEmpty(0))
      s << '.';
    else
      s << v[0];

    for (unsigned int i = 1; i < size; ++i)
      if (v.IsEmpty(i))
        s << separator << '.';
      else
        s << separator << v[i];

    return s;
  }


}


namespace boost {

  namespace lambda {

    template<typename T, typename U>
    class plain_return_type_2<
      arithmetic_action<multiply_action>,
      utl::SparseVector<T>,
      utl::SparseVector<U>
    > {
    private:
      // delegate to the action result of (T, Action, U)
      typedef typename
        return_type_2<
          arithmetic_action<multiply_action>,
          T,
          U
        >::type res_type;

    public:
      typedef res_type type;
    };

  }

}


#include <utl/SparseMatrixVectorOp.h>


#endif