#ifndef MORT_OBJECT_NUMVEC_HPP
#define MORT_OBJECT_NUMVEC_HPP

#define MORT_DONT_USE_VALARRAY indeed

#include <cassert>

#include <vector>
#include <valarray>

#ifdef MORT_DONT_USE_VALARRAY
#  define MORT_NUMVEC_PARENT std::vector<double>
#  include <numeric>
#else
#  define MORT_NUMVEC_PARENT std::valarray<double>
#endif // MORT_DONT_USE_VALARRAY


namespace mort
{
    class numvec : public MORT_NUMVEC_PARENT
    {

    public:

        typedef MORT_NUMVEC_PARENT parent_type;

        numvec() {}

        numvec(int n) : parent_type(n) {}

#       ifdef MORT_DONT_USE_VALARRAY
        numvec( const double& x, int n ) : parent_type(n, x) {}

        numvec( const double* x, int n )
        : parent_type(n)
        {
            std::copy(x, x + n, begin());
        }

        numvec(const numvec& v) { *this = v; }
#       else
        numvec( double n, int size ) : parent_type(n, size) {}
        numvec( const double* n, int size ) : parent_type(n, size) {}
#       endif // MORT_DONT_USE_VALARRAY

        numvec( const parent_type& rhs ) : parent_type( rhs ) {}

#       ifndef MORT_DONT_USE_VALARRAY
        numvec( const std::slice_array<double>& rhs ) :  parent_type( rhs ) {}
        template< typename T > numvec( const std::_Expr<T,double>& rhs ) : parent_type( rhs ) {}
#       endif // MORT_DONT_USE_VALARRAY

        numvec& operator=( const numvec& rhs ) 
        {
            if( this!=&rhs )
            {
                resize( rhs.size() );
                for( unsigned int i=0; i < size(); ++i )
                {
                    (*this)[i] = rhs[i];
                }
            }

            return *this;
        }

//        numvec subvec(int bgn, int end) { return (*this)[std::slice(bgn,end,1)]; }

#       ifdef MORT_DONT_USE_VALARRAY
        double sum() const
        {
            return std::accumulate(begin(), end(), 0);
        }

        double max() const
        {
            parent_type::const_iterator i = std::max_element(begin(), end());
             assert(i != end());
            return *i;
        }

        numvec& operator/=(const double& x)
        {
            assert(x != 0.0);
            for (parent_type::iterator i = begin(); i != end(); ++i)
                *i /= x;

            return *this;
        }

        numvec& operator*=(const double& x)
        {
            for (parent_type::iterator i = begin(); i != end(); ++i)
                *i *= x;

            return *this;
        }

        numvec& operator+=(const numvec& x)
        {
            assert(size() == x.size());

            if (this != &x) {
                for (size_t n(0); n < size(); ++n)
                    this->operator[](n) += x[n];
            } else {
                for (size_t n(0); n < size(); ++n)
                    (this)[n] *= 2;
            }

            return *this;
        }

        numvec& operator-=(const numvec& x)
        {
            assert(size() == x.size());

            if (this != &x) {
                for (size_t n(0); n < size(); ++n)
                    this->operator[](n) -= x[n];
            } else {
                std::fill(begin(), end(), 0.0);
            }

            return *this;
        }

        numvec operator-() const
        {
            numvec result(size());

            for (size_t n(0); n < size(); ++n)
                result[n] = -this->operator[](n);

            return result;
        }

        using parent_type::operator[];

        numvec operator[](const std::slice& s) const
        {
            numvec result(s.size());
            for (size_t n(0); n < s.size(); ++n)
                result[n] = (*this)[s.start() + n*s.stride()];
            return result;
        }
#       endif // MORT_DONT_USE_VALARRAY
    };
  
    /// construct numvec from elements
    numvec makevec(double x, double y);
    
    numvec makevec(double x, double y, double z);

    numvec makevec(double x, double y, double z, double s);
 
    numvec makevec(double x, double y, double z, double s, double u);
        
    numvec makevec(double x, double y, double z, double s, double u, double v);        
 
    // maximum element of a numvec
    inline double max( const numvec& v ) { return v.max(); }
 
    /// norm of a vector 
    double norm( const numvec& vec );
  
    /// dot product of two vectors.
    double dotpd( const numvec& v1, const numvec& v2 );

    /// cross product of two vector
    numvec cross( const numvec& v1, const numvec& v2 );
 
    /// normalize a vector in place
    numvec& normalize( numvec& vec );
 
    /// return a normalized copy of a vector, the orginal one remain unchanged.
    numvec  normalcpy( const numvec& v );    

#   ifdef MORT_DONT_USE_VALARRAY

    // the following does *not* implement all valarray's operators;
    // just those needed by the rest of the gleap are here

    inline numvec operator-(const numvec& a, const numvec& b)
    {
        assert(a.size() == b.size());

        numvec result(a.size());
        for (size_t n(0); n < a.size(); ++n)
            result[n] = a[n] - b[n];
        return result;
    }

    inline numvec operator+(const numvec& a, const numvec& b)
    {
        assert(a.size() == b.size());

        numvec result(a.size());
        for (size_t n(0); n < a.size(); ++n)
            result[n] = a[n] + b[n];
        return result;
    }

    inline numvec operator*(const numvec& a, const double& b)
    {
        numvec result(a.size());
        for (size_t n(0); n < a.size(); ++n)
            result[n] = a[n]*b;
        return result;
    }

    inline numvec operator*(const double& a, const numvec& b)
    {
        return b*a;
    }

    inline numvec operator/(const numvec& a, const double& b)
    {
        assert(b != 0.0);
        numvec result(a.size());
        for (size_t n(0); n < a.size(); ++n)
            result[n] = a[n]/b;
        return result;
    }

#   endif // MORT_DONT_USE_VALARRAY

} // namespace mort

#endif