// $Id$

#ifndef _BOOST_UBLAS_COUNT_STRUCTURAL_ELEMENTS_
#define _BOOST_UBLAS_COUNT_STRUCTURAL_ELEMENTS_

// Boost numeric uBlas library headers
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/storage.hpp>

namespace boost { namespace numeric { namespace ublas {


// Implementation of simple free unary function count_structural_elements() to iterate over all
// structural elements in a vector_expression<E> or matrix_expression<E>.
// These function are not particularly useful, but demonstrate how to iterate over only the
// structural elements, leading to a performance benefit in packed and sparse storage containers.



// Unary functor returning size_type scalar 
template<class T>
struct scalar_size_unary_functor {
  typedef std::size_t size_type;
  typedef std::ptrdiff_t difference_type;
  typedef T value_type;
  typedef size_type result_type;
};

// Unary vector_count_structural_elements functor for vector expressions
template<class T>
struct vector_count_structural_elements:
  public scalar_size_unary_functor<T> {
  typedef typename scalar_size_unary_functor<T>::size_type size_type;
  typedef typename scalar_size_unary_functor<T>::difference_type difference_type;
  typedef typename scalar_size_unary_functor<T>::value_type value_type;
  typedef typename scalar_size_unary_functor<T>::result_type result_type;

  template<class E>
  static BOOST_UBLAS_INLINE
  result_type apply (const vector_expression<E> &e) {
		// Initialise counter
		size_type nnz = 0;
		// Iterate over vector elements and count number of elements containing non-zero value
    size_type size( e().size () );
    for (size_type i = 0; i < size; ++ i) {
 			// Update count of structural elements
			nnz++;
    }
		// Return number of elements containing non-zero value
		return nnz;
	}
  // Dense case
  template<class I>
  static BOOST_UBLAS_INLINE
  result_type apply( difference_type size, I it ) {
		// Initialise counter
    size_type nnz = 0;
		// Iterate over vector elements and count number of elements containing non-zero value
    while (--size >= 0) {
  		// Update count of structural elements
			nnz++;
			// Advance iterator
			++ it;
    }
		// Return number of elements containing non-zero value
		return nnz;
	}
	// Sparse case
	template<class I>
	static BOOST_UBLAS_INLINE
	result_type apply( I it, const I &it_end ) {
		// Initialise counter
    size_type nnz = 0;
		// Iterate over vector elements and count number of elements containing non-zero value
		while (it != it_end) {
  		// Update count of structural elements
			nnz++;
			// Advance iterator
			++ it;
		}
		// Return number of elements containing non-zero value
		return nnz;
	}
};

// Free function to count number of structural elements in vector
template<class E>
BOOST_UBLAS_INLINE
typename vector_scalar_unary_traits< E, vector_count_structural_elements< typename E::value_type > >::result_type
count_structural_elements (const vector_expression<E> &e)
{
	typedef typename vector_scalar_unary_traits<E, vector_count_structural_elements<typename E::value_type> >::expression_type expression_type;
   
	return expression_type(e());
}



template<class T>
struct matrix_count_structural_elements:
  public scalar_size_unary_functor<T> {
  typedef typename scalar_size_unary_functor<T>::size_type size_type;
  typedef typename scalar_size_unary_functor<T>::difference_type difference_type;
  typedef typename scalar_size_unary_functor<T>::value_type value_type;
  typedef typename scalar_size_unary_functor<T>::result_type result_type;

	// Specialisation for matrix_reference<C>
	// A matrix reference is associated with a matrix_container which may use any of the following
	// storage types (dense, dense_proxy, packed, packed_proxy, sparse, sparse_proxy).
	// We therefore use a dispatcher to select the optimised iteration code.
	// Dispatcher
  template<class C>
  static BOOST_UBLAS_INLINE
  result_type apply( const matrix_reference<C> &m) {
		// Call specialised function
		return apply( m, C::storage_category() );
  }
  // Dense storage
  template<class C>
  static BOOST_UBLAS_INLINE
	result_type apply( const matrix_reference<C> &m, dense_proxy_tag ) {
		// Initialise counter
		size_type nnz = 0;
		// Iterate over dense matrix and count number of elements containing non-zero value
		for (size_type i = 0; i < m.size1(); i++) {
			for (size_type j = 0; j < m.size2(); j++) {
				// Update count of structural elements
				nnz++;
			}
		}
		// Return number of elements containing non-zero value
		return nnz;
  }
	// Packed storage
	template <class C>
  static BOOST_UBLAS_INLINE
	result_type apply( const matrix_reference<C> &m, packed_proxy_tag ) {
		typedef typename matrix_reference<C>::const_iterator1 const_iterator1;
		typedef typename matrix_reference<C>::const_iterator2 const_iterator2;
		// Initialise counter
		size_type nnz = 0;
		// Iterate over packed matrix and count number of elements containing non-zero value
		// Initialise it1 iterator
		const_iterator1 it1( m.begin1() );
		// Compute data size
		difference_type size1( m.end1() - it1 );
		// Outer loop
		while (--size1 >= 0) {
	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
			// Initialise it2 iterator
			const_iterator2 it2( it1.begin() );
			// Compute data size
			difference_type size2( it1.end() - it2 );
	#else
			// Initialise it2 iterator
			const_iterator2 it2( begin( it1, iterator1_tag () ) );
			// Compute number of elements in this sequence
			difference_type size2( end( it1, iterator1_tag() ) - it2 );
	#endif
			// Inner loop
			while (--size2 >= 0) {
				// Update count of structural elements
				nnz++;
				// Advance inner iterator
				++it2;
			}
			++ it1; // Advance outer iterator
		}
		// Return number of elements containing non-zero value
		return nnz;
	}
	// Sparse storage
	template <class C>
  static BOOST_UBLAS_INLINE
	result_type apply( const matrix_reference<C> &m, sparse_proxy_tag ) {
		typedef typename matrix_reference<C>::const_iterator1 const_iterator1;
		typedef typename matrix_reference<C>::const_iterator2 const_iterator2;
		// Initialise counter
		size_type nnz = 0;
		// Iterate over packed matrix and count number of elements containing non-zero value
		// Initialise it1 iterators
		const_iterator1 it1( m.begin1() );
		const_iterator1 it1_end( m.end1() );

		// Iterate over non-zero elements and populate arrays of indices
		// Outer loop
		while (it1 != it1_end) {
#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
			// Initialise it2 iterators
			const_iterator2 it2( it1.begin() );
			const_iterator2 it2_end( it1.end() );
#else
			// Initialise it2 iterators
			const_iterator2 it2( begin( it1, iterator1_tag() ) );
			const_iterator2 it2_end( end( it1, iterator1_tag() ) );
#endif // BOOST_UBLAS_NO_NESTED_CLASS_RELATION
			// Inner loop
			while (it2 != it2_end) {
				// Update count of structural elements
				nnz++;
				// Advance inner iterator
				++it2;
			}
			++ it1; // Advance outer iterator
		}
		// Return number of elements containing non-zero value
		return nnz;
	}
	// Specialisation for matrix_expression<E>
	// If the matrix_reference<C> specialisations have not been selected by the compiler, then we
	// have an matrix_expression that needs to be evaluated. As such, there are no optimisations to
	// the iterator relevent here.
	template<class E>
  static BOOST_UBLAS_INLINE
  result_type apply (const matrix_expression<E> &e) {
		// Initialise counter
		size_type nnz = 0;
		// Get number of rows
    size_type size2 (e ().size2 ());
		// Iterate over rows
    for (size_type j = 0; j < size2; ++ j) {
			// Get data size
      size_type size1 (e ().size1 ());
			// Inner loop
      for (size_type i = 0; i < size1; ++ i) {
        // size_type NNZ2( type_traits<value_type>::equals( e()(i, j) ) );
 				nnz++;
     }
    }
    return nnz; 
  }
};



// Free function to count number of non-zero elements in matrix
// real: nnz_count2(m) = (sum(m(i,j) != 0))
// complex: ?
template<class E>
BOOST_UBLAS_INLINE
typename matrix_scalar_unary_traits< E, matrix_count_structural_elements< typename E::value_type > >::result_type
count_structural_elements (const matrix_expression<E> &e)
{
	typedef typename matrix_scalar_unary_traits<E, matrix_count_structural_elements<typename E::value_type> >::expression_type expression_type;
   
	return expression_type(e());
}


}}} // end of namespace scope

#endif // _BOOST_UBLAS_COUNT_STRUCTURAL_ELEMENTS_