Vector.h

// ===========================================================================
// Imagine++ Libraries
// Copyright (C) Imagine
// For detailed information: http://imagine.enpc.fr/software
// ===========================================================================

namespace Imagine {

    template <typename T> class Matrix;
    template <typename T> class SymMatrix;


    template <typename T>
    class Vector : public Array<T> {
        typedef Array<T> Base;
    public:
        typedef typename Base::iterator iterator;
        typedef typename Base::const_iterator const_iterator;
        Vector() : Base() {}
        explicit Vector(size_t n) : Base(n) {}
        Vector(const Base& v) : Base(v) {}
        Vector(T* ptr, size_t n,bool handleDelete=false):Base(ptr,n,handleDelete) {}
        ~Vector() {}
        Vector& operator=(const Vector& v) { Base::operator=(v); return *this; }
        Vector clone() const
        {
            Vector v(this->size());
            equalize((int)this->size(),this->data(),1,v.data(),1);
            return v;
        }
        Vector& fill(T x) { Base::fill(x); return *this; }
        static Vector Zero(size_t n) {
            Vector Z(n);
            memset(Z.data(),0,n*sizeof(T));
            return Z;
        }
        Vector getSubVect(size_t offset, size_t size) const {
            return Vector(Base::getSubArray(offset,size));
        }
        Vector getSubVectRef(size_t offset, size_t size) {
            assert(offset+size<=this->size());
            return Vector(this->data()+offset,size);
        }
        Vector operator+(const Vector& v) const {
            assert(this->size()==v.size());
            Vector p=this->clone();
            return (p+=v);
        }
        Vector& operator+=(const Vector& v) {
            assert(this->size()==v.size());
            combine((int)this->size(), 1, v.data(), this->data());
            return *this;
        }
        Vector operator-(const Vector& v) const {
            assert(this->size()==v.size());
            Vector p=this->clone();
            return (p-=v);
        }
        Vector& operator-=(const Vector& v) {
            assert(this->size()==v.size());
            combine((int)this->size(), -1, v.data(), this->data());
            return *this;
        }
        Vector operator+(T x) const {
            Vector p=this->clone();
            return (p+=x);
        }
        Vector& operator+=(T x) {
            for (iterator it = this->begin() ; it != this->end() ; ++it)
                *it+=x;
            return *this;
        }
        Vector operator-(T x) const {
            Vector p=this->clone();
            return (p-=x);
        }
        Vector& operator-=(T x) {
            for (iterator it = this->begin() ; it != this->end() ; ++it)
                *it-=x;
            return *this;
        }
        Vector operator - () const {
            return (*this)*T(-1);
        }
        Vector operator*(T x) const {
            Vector p=this->clone();
            return (p*=x);
        }
        Vector& operator*=(T x) {
            multiply((int)this->size(), x, this->data());
            return *this;
        }
        Vector operator/(T x) const {
            return (*this)*(1/x);
        }       
        Vector& operator/=(T x) {
            (*this)*=(1/x);
            return *this;
        }
        friend inline Vector operator * (T x, const Vector& v) {return v*x;}
        friend inline Vector operator + (T x, const Vector& v) {return v+x;}
        friend inline Vector operator - (T x, const Vector& v) {return (-v)+x;}
        T operator*(const Vector& v) const {
            assert(this->size()==v.size());
            return scalarProduct((int)this->size(), this->data(), v.data());
        }
        friend T norm(const Vector& v) {
            return squaredNorm((int) v.size(), v.data());
        }
        friend T norm2(const Vector& v) {
            return norm(v)*norm(v);
        }
        Vector& normalize() {
            T n=norm(*this);
            if (n!=0) return operator /= (n);
            return *this;
        }
        friend inline Vector normalized(const Vector& v) {
            T n=norm(v);
            assert(n!=0);
            return v/n;
        }
        friend inline T maxNorm(const Vector& v) {
            if (v.empty()) return 0;
            const_iterator it = v.begin();
            T mx=std::abs(*it);
            for ( ++it ; it != v.end() ; ++it) {
                T a = std::abs(*it);
                if (a>mx)
                    mx=a;
            }
            return mx;
        }
        friend inline T sum(const Vector& v) {
            return std::accumulate(v.begin(),v.end(),T(0));
        }
        friend inline T mean(const Vector& v) {
            return sum(v)/v.size();
        }
        friend Vector convolution(const Vector& K,const Vector& v)  {
            if (v.size()<K.size()) return convolution(v,K);
            Vector p(K.size()+v.size()-1);
            p.fill(T(0));
            for (size_t i=0; i<v.size(); i++) {
                combine((int)K.size(), v[i], K.data(), p.data()+i);
            }
            return p;
        }
        friend Vector truncConvolution(const Vector& K,const Vector& v, size_t n=0)  {
            assert(n<K.size());
            Vector p(v.size());
            p.fill(T(0));
            for (size_t i=0; i<v.size(); i++) {
                int start = std::max((int)(0),(int)(n)-(int)(i));
                int end = std::min((int)(K.size()),(int)(v.size()+n)-(int)(i));
                combine((end - start), v[i], K.data()+start, p.data() + start + (i - n));
            }
            return p;
        }
        friend inline T correlation(const Vector& V1,const Vector& V2) {
            Vector v1=V1-mean(V1);
            Vector v2=V2-mean(V2);
            return v1*v2;
        }
        friend inline T normalizedCorrelation(const Vector& V1,const Vector& V2) {// Non optimized version
            T c1=correlation(V1,V1);
            if (c1==0)
                return 0;
            T c2=correlation(V2,V2);
            if (c2==0)
                return 0;
            return correlation(V1,V2)/::sqrt(c1*c2);
        }   
        friend inline double dist(const Vector& V1,const Vector& V2) {
            double d=0;
            for (size_t i=0;i<V1.size();i++) {
                double e=double(V1[i])-double(V2[i]);
                d+=e*e;
            }
            return ::sqrt(d);
        }
    };


}