// $Id $

// Simple application to test vector quad relational operator implementation

#include <iostream>
#include <complex>


// Boost numeric uBlas library
#include <boost/numeric/ublas/traits.hpp>
#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/triangular.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/storage.hpp>
#include <boost/numeric/ublas/functional.hpp>

namespace ublas = boost::numeric::ublas;



namespace boost { namespace numeric { namespace ublas {


// Vector expression used to store quad relational op expression
template<class E1, class E2, class R1, class R2, class F>
class vector_quad_relop:
    public vector_expression<vector_quad_relop<E1, E2, R1, R2, F> > {

    typedef E1 expression1_type;
    typedef E2 expression2_type;
		typedef R1 expression3_type;
		typedef R2 expression4_type;

    typedef F functor_type;
    typedef typename E1::const_closure_type expression1_closure_type;
    typedef typename E2::const_closure_type expression2_closure_type;
    typedef typename R1::const_closure_type expression3_closure_type;
    typedef typename R2::const_closure_type expression4_closure_type;
    typedef vector_quad_relop<E1, E2, R1, R2, F> self_type;
public:
#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
    using vector_expression<vector_quad_relop<E1, E2, R1, R2, F> >::operator ();
#endif
    typedef typename promote_traits<typename E1::size_type, typename E2::size_type>::promote_type size_type;
    typedef typename promote_traits<typename E1::difference_type, typename E2::difference_type>::promote_type difference_type;
    typedef typename F::result_type value_type;
    typedef value_type const_reference;
    typedef const_reference reference;
    typedef const self_type const_closure_type;
    typedef const_closure_type closure_type;
    typedef unknown_storage_tag storage_category;

    // Construction and destruction
    BOOST_UBLAS_INLINE
    vector_quad_relop (const expression1_type &e1, const expression2_type &e2, const expression3_type &r1, const expression4_type &r2):
        e1_ (e1), e2_ (e2), r1_ (r1), r2_ (r2) {}

    // Accessors
    BOOST_UBLAS_INLINE
    size_type size () const {
        return BOOST_UBLAS_SAME (e1_.size (), e2_.size ()); 
    }

private:
    // Accessors
    BOOST_UBLAS_INLINE
    const expression1_closure_type &expression1 () const {
        return e1_;
    }
    BOOST_UBLAS_INLINE
    const expression2_closure_type &expression2 () const {
        return e2_;
    }
    BOOST_UBLAS_INLINE
    const expression3_closure_type &expression3 () const {
        return r1_;
    }
    BOOST_UBLAS_INLINE
    const expression4_closure_type &expression4 () const {
        return r2_;
    }

public:
    // Element access - applies functor to elements
    BOOST_UBLAS_INLINE
    const_reference operator () (size_type i) const {
        return functor_type::apply (e1_ (i), e2_ (i), r1_(), r2_());
    }

    BOOST_UBLAS_INLINE
    const_reference operator [] (size_type i) const {
        return functor_type::apply (e1_ [i], e2_ [i], r1_(), r2_());
    }

    // Closure comparison
    BOOST_UBLAS_INLINE
    bool same_closure (const vector_quad_relop &vb) const {
        return (*this).expression1 ().same_closure (vb.expression1 ()) &&
                (*this).expression2 ().same_closure (vb.expression2 ());
    }

    // Iterator types
private:
    typedef typename E1::const_iterator const_subiterator1_type;
    typedef typename E2::const_iterator const_subiterator2_type;
    typedef const value_type *const_pointer;

public:
#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
    typedef typename iterator_restrict_traits<typename const_subiterator1_type::iterator_category,
                                              typename const_subiterator2_type::iterator_category>::iterator_category iterator_category;
    typedef indexed_const_iterator<const_closure_type, iterator_category> const_iterator;
    typedef const_iterator iterator;
#else
    class const_iterator;
    typedef const_iterator iterator;
#endif

    // Element lookup
    BOOST_UBLAS_INLINE
    const_iterator find (size_type i) const {
        const_subiterator1_type it1 (e1_.find (i));
        const_subiterator1_type it1_end (e1_.find (size ()));
        const_subiterator2_type it2 (e2_.find (i));
        const_subiterator2_type it2_end (e2_.find (size ()));
        i = (std::min) (it1 != it1_end ? it1.index () : size (),
                      it2 != it2_end ? it2.index () : size ());
#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
        return const_iterator (*this, i);
#else
        return const_iterator (*this, i, it1, it1_end, it2, it2_end);
#endif
    }

    // Iterator merges the iterators of the referenced expressions and
    // enhances them with the binary functor.

#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
    class const_iterator:
        public container_const_reference<vector_quad_relop>,
        public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
                                                                      typename E2::const_iterator::iterator_category>::iterator_category>::template
                    iterator_base<const_iterator, value_type>::type {
    public:
        typedef typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
                                                  typename E2::const_iterator::iterator_category>::iterator_category iterator_category;
        typedef typename vector_quad_relop::difference_type difference_type;
        typedef typename vector_quad_relop::value_type value_type;
        typedef typename vector_quad_relop::const_reference reference;
        typedef typename vector_quad_relop::const_pointer pointer;

        // Construction and destruction
        BOOST_UBLAS_INLINE
        const_iterator ():
            container_const_reference<self_type> (), i_ (), it1_ (), it1_end_ (), it2_ (), it2_end_ () {}
        BOOST_UBLAS_INLINE
        const_iterator (const self_type &vb, size_type i,
                        const const_subiterator1_type &it1, const const_subiterator1_type &it1_end,
                        const const_subiterator2_type &it2, const const_subiterator2_type &it2_end):
            container_const_reference<self_type> (vb), i_ (i), it1_ (it1), it1_end_ (it1_end), it2_ (it2), it2_end_ (it2_end) {}

    private: 
        // Dense specializations
        BOOST_UBLAS_INLINE
        void increment (dense_random_access_iterator_tag) {
            ++ i_; ++ it1_; ++ it2_;
        }
        BOOST_UBLAS_INLINE
        void decrement (dense_random_access_iterator_tag) {
            -- i_; -- it1_; -- it2_;
        }
        BOOST_UBLAS_INLINE
        void increment (dense_random_access_iterator_tag, difference_type n) {
            i_ += n; it1_ += n; it2_ += n;
        }
        BOOST_UBLAS_INLINE
        void decrement (dense_random_access_iterator_tag, difference_type n) {
            i_ -= n; it1_ -= n; it2_ -= n;
        }
        BOOST_UBLAS_INLINE
        value_type dereference (dense_random_access_iterator_tag) const {
            return functor_type::apply (*it1_, *it2_);
        }

        // Packed specializations
        BOOST_UBLAS_INLINE
        void increment (packed_random_access_iterator_tag) {
            if (it1_ != it1_end_)
                if (it1_.index () <= i_)
                    ++ it1_;
            if (it2_ != it2_end_)
                if (it2_.index () <= i_)
                    ++ it2_;
            ++ i_;
        }
        BOOST_UBLAS_INLINE
        void decrement (packed_random_access_iterator_tag) {
            if (it1_ != it1_end_)
                if (i_ <= it1_.index ())
                    -- it1_;
            if (it2_ != it2_end_)
                if (i_ <= it2_.index ())
                    -- it2_;
            -- i_;
        }
        BOOST_UBLAS_INLINE
        void increment (packed_random_access_iterator_tag, difference_type n) {
            while (n > 0) {
                increment (packed_random_access_iterator_tag ());
                --n;
            }
            while (n < 0) {
                decrement (packed_random_access_iterator_tag ());
                ++n;
            }
        }
        BOOST_UBLAS_INLINE
        void decrement (packed_random_access_iterator_tag, difference_type n) {
            while (n > 0) {
                decrement (packed_random_access_iterator_tag ());
                --n;
            }
            while (n < 0) {
                increment (packed_random_access_iterator_tag ());
                ++n;
            }
        }
        BOOST_UBLAS_INLINE
        value_type dereference (packed_random_access_iterator_tag) const {
            value_type t1 = value_type/*zero*/();
            if (it1_ != it1_end_)
                if (it1_.index () == i_)
                    t1 = *it1_;
            value_type t2 = value_type/*zero*/();
            if (it2_ != it2_end_)
                if (it2_.index () == i_)
                    t2 = *it2_;
            return functor_type::apply (t1, t2);
        }

        // Sparse specializations
        BOOST_UBLAS_INLINE
        void increment (sparse_bidirectional_iterator_tag) {
            size_type index1 = (*this) ().size ();
            if (it1_ != it1_end_) {
                if  (it1_.index () <= i_)
                    ++ it1_;
                if (it1_ != it1_end_)
                    index1 = it1_.index ();
            }
            size_type index2 = (*this) ().size ();
            if (it2_ != it2_end_) {
                if (it2_.index () <= i_)
                    ++ it2_;
                if (it2_ != it2_end_)
                    index2 = it2_.index ();
            }
            i_ = (std::min) (index1, index2);
        }
        BOOST_UBLAS_INLINE
        void decrement (sparse_bidirectional_iterator_tag) {
            size_type index1 = (*this) ().size ();
            if (it1_ != it1_end_) {
                if (i_ <= it1_.index ())
                    -- it1_;
                if (it1_ != it1_end_)
                    index1 = it1_.index ();
            }
            size_type index2 = (*this) ().size ();
            if (it2_ != it2_end_) {
                if (i_ <= it2_.index ())
                    -- it2_;
                if (it2_ != it2_end_)
                    index2 = it2_.index ();
            }
            i_ = (std::max) (index1, index2);
        }
        BOOST_UBLAS_INLINE
        void increment (sparse_bidirectional_iterator_tag, difference_type n) {
            while (n > 0) {
                increment (sparse_bidirectional_iterator_tag ());
                --n;
            }
            while (n < 0) {
                decrement (sparse_bidirectional_iterator_tag ());
                ++n;
            }
        }
        BOOST_UBLAS_INLINE
        void decrement (sparse_bidirectional_iterator_tag, difference_type n) {
            while (n > 0) {
                decrement (sparse_bidirectional_iterator_tag ());
                --n;
            }
            while (n < 0) {
                increment (sparse_bidirectional_iterator_tag ());
                ++n;
            }
        }
        BOOST_UBLAS_INLINE
        value_type dereference (sparse_bidirectional_iterator_tag) const {
            value_type t1 = value_type/*zero*/();
            if (it1_ != it1_end_)
                if (it1_.index () == i_)
                    t1 = *it1_;
            value_type t2 = value_type/*zero*/();
            if (it2_ != it2_end_)
                if (it2_.index () == i_)
                    t2 = *it2_;
            return functor_type::apply (t1, t2);
        }

    public: 
        // Arithmetic
        BOOST_UBLAS_INLINE
        const_iterator &operator ++ () {
            increment (iterator_category ());
            return *this;
        }
        BOOST_UBLAS_INLINE
        const_iterator &operator -- () {
            decrement (iterator_category ());
            return *this;
        }
        BOOST_UBLAS_INLINE
        const_iterator &operator += (difference_type n) {
            increment (iterator_category (), n);
            return *this;
        }
        BOOST_UBLAS_INLINE
        const_iterator &operator -= (difference_type n) {
            decrement (iterator_category (), n);
            return *this;
        }
        BOOST_UBLAS_INLINE
        difference_type operator - (const const_iterator &it) const {
            BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
            return index () - it.index ();
        }

        // Dereference
        BOOST_UBLAS_INLINE
        const_reference operator * () const {
            return dereference (iterator_category ());
        }
        BOOST_UBLAS_INLINE
        const_reference operator [] (difference_type n) const {
            return *(*this + n);
        }

        // Index
        BOOST_UBLAS_INLINE
        size_type index () const {
            return i_;
        }

        // Assignment
        BOOST_UBLAS_INLINE
        const_iterator &operator = (const const_iterator &it) {
            container_const_reference<self_type>::assign (&it ());
            i_ = it.i_;
            it1_ = it.it1_;
            it1_end_ = it.it1_end_;
            it2_ = it.it2_;
            it2_end_ = it.it2_end_;
            return *this;
        }

        // Comparison
        BOOST_UBLAS_INLINE
        bool operator == (const const_iterator &it) const {
            BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
            return index () == it.index ();
        }
        BOOST_UBLAS_INLINE
        bool operator < (const const_iterator &it) const {
            BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
            return index () < it.index ();
        }

    private:
        size_type i_;
        const_subiterator1_type it1_;
        const_subiterator1_type it1_end_;
        const_subiterator2_type it2_;
        const_subiterator2_type it2_end_;
    };
#endif

    BOOST_UBLAS_INLINE
    const_iterator begin () const {
        return find (0);
    }
    BOOST_UBLAS_INLINE
    const_iterator end () const {
        return find (size ());
    }

    // Reverse iterator
    typedef reverse_iterator_base<const_iterator> const_reverse_iterator;

    BOOST_UBLAS_INLINE
    const_reverse_iterator rbegin () const {
        return const_reverse_iterator (end ());
    }
    BOOST_UBLAS_INLINE
    const_reverse_iterator rend () const {
        return const_reverse_iterator (begin ());
    }

private:
    expression1_closure_type e1_;
    expression2_closure_type e2_;
    expression3_closure_type r1_;
    expression4_closure_type r2_;
};

template<class E1, class E2, class R1, class R2, class F>
struct vector_quad_relop_traits {
    typedef vector_quad_relop<E1, E2, R1, R2, F> expression_type;
#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
    typedef expression_type result_type; 
#else
    typedef typename E1::vector_temporary_type result_type;
#endif
};



// Scalar relational operator functor taking four arguments
// Arguments 3 and 4 must be the same type
// Result is the same type as arguments 3 and 4.
template<class E1, class E2, class R1, class R2>
struct scalar_rel_op_with_4args_functor
{
	typedef typename type_traits<E1>::const_reference argument1_type;
	typedef typename type_traits<E2>::const_reference argument2_type;
	typedef typename type_traits<R1>::const_reference argument3_type;
	typedef typename type_traits<R2>::const_reference argument4_type;
	typedef typename promote_traits<typename R1, typename R2>::promote_type result_type;
};



// Binary relational functors
// The relational operators are <, >, <=, >=, ==, and !=.
// Relational operators perform element-by-element comparisons between two vector expressions (or
// two matrix expressions). They return a boolean vector (or boolean matrix) of the same size, with
// elements set to true where the relation is true, and elements set to false where it is not.
// The operators <, >, <=, and >= use only the real part of their operands for the comparison.
// The operators == and != test real and imaginary parts.

// scalar_lt_quad_quad functor
template<class E1, class E2, class R1, class R2>
struct scalar_lt_quad
	: public scalar_rel_op_with_4args_functor<E1, E2, R1, R2>
{
	typedef typename scalar_rel_op_with_4args_functor<E1, E2, R1, R2>::argument1_type argument1_type;
  typedef typename scalar_rel_op_with_4args_functor<E1, E2, R1, R2>::argument2_type argument2_type;
  typedef typename scalar_rel_op_with_4args_functor<E1, E2, R1, R2>::argument3_type argument3_type;
  typedef typename scalar_rel_op_with_4args_functor<E1, E2, R1, R2>::argument4_type argument4_type;
  typedef typename scalar_rel_op_with_4args_functor<E1, E2, R1, R2>::result_type result_type;

  static BOOST_UBLAS_INLINE
  result_type apply (argument1_type t1, argument2_type t2, argument3_type r1, argument4_type r2) {
		return (ublas::type_traits< T1 >::real(t1)) < (ublas::type_traits< T2 >::real(t2)) ? r1 : r2;
  }
};

// --- Free quad functions for relational operators --- 

// Free quad function to model relational operator < for vector expressions
// The relational operators are <, >, <=, >=, ==, and !=.
// Performs element-by-element comparisons between two vector expressions e1 and e2
// Returns a vector expression of the same size as e1, with
// elements r1 where the relation is true, and elements set to r2 where it is not.
// Uses only the real part e1 and e2 for the comparison.
template<class E1, class E2, class R1, class R2>
BOOST_UBLAS_INLINE
typename vector_quad_relop_traits< E1, E2, R1, R2, scalar_lt_quad< typename E1::value_type, typename E2::value_type, typename R1::value_type, typename R2::value_type > >::result_type
lt_quad (const vector_expression<E1> &e1, const vector_expression<E2> &e2, const scalar_expression<R1> &r1, const scalar_expression<R2> &r2 )
{
	typedef typename vector_quad_relop_traits<E1, E2, R1, R2, scalar_lt_quad< typename E1::value_type, typename E2::value_type, typename R1::value_type, typename R2::value_type > >::expression_type expression_type;
   
	// Evaluate expression by calling scalar_lt_quad functor apply() method on elements in lock-step
	return expression_type(e1(), e2(), r1(), r2());
}


}}} // end namespace scope


// Declarations
int main( int argc, char **argv );

// Function definitions
int main( int argc, char **argv )
{
	ublas::vector<double> V1(2);
	V1(0) = 0.0;
	V1(1) = 1.0;
	ublas::vector<double> V2(2);
	V2(0) = 0.0;
	V2(1) = 1.0;

	// bool r1(true), r2(false);
	ublas::scalar_value<bool> r1(true), r2(false);

	std::cout << "lt_quad( V1, V2, r1, r2) = " << lt_quad( V1, V2, r1, r2) << std::endl;

	return 0;
}