IT++ Logo

vec.h

Go to the documentation of this file.
00001 
00030 #ifndef VEC_H
00031 #define VEC_H
00032 
00033 #ifndef _MSC_VER
00034 #  include <itpp/config.h>
00035 #else
00036 #  include <itpp/config_msvc.h>
00037 #endif
00038 
00039 #include <itpp/base/itassert.h>
00040 #include <itpp/base/math/misc.h>
00041 #include <itpp/base/copy_vector.h>
00042 #include <itpp/base/factory.h>
00043 
00044 
00045 namespace itpp {
00046 
00047   // Declaration of Vec
00048   template<class Num_T> class Vec;
00049   // Declaration of Mat
00050   template<class Num_T> class Mat;
00051   // Declaration of bin
00052   class bin;
00053 
00054   //-----------------------------------------------------------------------------------
00055   // Declaration of Vec Friends
00056   //-----------------------------------------------------------------------------------
00057 
00059   template<class Num_T> const Vec<Num_T> operator+(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
00061   template<class Num_T> const Vec<Num_T> operator+(const Vec<Num_T> &v, const Num_T t);
00063   template<class Num_T> const Vec<Num_T> operator+(const Num_T t, const Vec<Num_T> &v);
00064 
00066   template<class Num_T> const Vec<Num_T> operator-(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
00068   template<class Num_T> const Vec<Num_T> operator-(const Vec<Num_T> &v, const Num_T t);
00070   template<class Num_T> const Vec<Num_T> operator-(const Num_T t, const Vec<Num_T> &v);
00072   template<class Num_T> const Vec<Num_T> operator-(const Vec<Num_T> &v);
00073 
00075   template<class Num_T> Num_T dot(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
00077   template<class Num_T> Num_T operator*(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
00086   template<class Num_T>
00087   const Mat<Num_T> outer_product(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
00088          bool hermitian = false);
00090   template<class Num_T> const Vec<Num_T> operator*(const Vec<Num_T> &v, const Num_T t);
00092   template<class Num_T> const Vec<Num_T> operator*(const Num_T t, const Vec<Num_T> &v);
00093 
00095   template<class Num_T> const Vec<Num_T> elem_mult(const Vec<Num_T> &a, const Vec<Num_T> &b);
00097   template<class Num_T> const Vec<Num_T> elem_mult(const Vec<Num_T> &a, const Vec<Num_T> &b, const Vec<Num_T> &c);
00099   template<class Num_T> const Vec<Num_T> elem_mult(const Vec<Num_T> &a, const Vec<Num_T> &b, const Vec<Num_T> &c, const Vec<Num_T> &d);
00100 
00102   template<class Num_T> void elem_mult_out(const Vec<Num_T> &a, const Vec<Num_T> &b, Vec<Num_T> &out);
00104   template<class Num_T> void elem_mult_out(const Vec<Num_T> &a, const Vec<Num_T> &b, const Vec<Num_T> &c, Vec<Num_T> &out);
00106   template<class Num_T> void elem_mult_out(const Vec<Num_T> &a, const Vec<Num_T> &b, const Vec<Num_T> &c, const Vec<Num_T> &d, Vec<Num_T> &out);
00107 
00109   template<class Num_T> void elem_mult_inplace(const Vec<Num_T> &a, Vec<Num_T> &b);
00111   template<class Num_T> Num_T elem_mult_sum(const Vec<Num_T> &a, const Vec<Num_T> &b);
00112 
00114   template<class Num_T> const Vec<Num_T> operator/(const Vec<Num_T> &v, const Num_T t);
00116   template<class Num_T> const Vec<Num_T> operator/(const Num_T t, const Vec<Num_T> &v);
00117 
00119   template<class Num_T> const Vec<Num_T> elem_div(const Vec<Num_T> &a, const Vec<Num_T> &b);
00121   template<class Num_T> const Vec<Num_T> elem_div(const Num_T t, const Vec<Num_T> &v);
00123   template<class Num_T> void elem_div_out(const Vec<Num_T> &a, const Vec<Num_T> &b, Vec<Num_T> &out);
00125   template<class Num_T> Num_T elem_div_sum(const Vec<Num_T> &a, const Vec<Num_T> &b);
00126 
00128   template<class Num_T> const Vec<Num_T> concat(const Vec<Num_T> &v, const Num_T a);
00130   template<class Num_T> const Vec<Num_T> concat(const Num_T a, const Vec<Num_T> &v);
00132   template<class Num_T> const Vec<Num_T> concat(const Vec<Num_T> &v1,const Vec<Num_T> &v2);
00134   template<class Num_T> const Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2, const Vec<Num_T> &v3);
00136   template<class Num_T> const Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2, const Vec<Num_T> &v3, const Vec<Num_T> &v4);
00138   template<class Num_T> const Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2, const Vec<Num_T> &v3, const Vec<Num_T> &v4, const Vec<Num_T> &v5);
00139 
00140   //-----------------------------------------------------------------------------------
00141   // Declaration of Vec
00142   //-----------------------------------------------------------------------------------
00143 
00207   template<class Num_T>
00208   class Vec {
00209   public:
00211     typedef Num_T value_type;
00212 
00214     explicit Vec(const Factory &f = DEFAULT_FACTORY);
00216     explicit Vec(int size, const Factory &f = DEFAULT_FACTORY);
00218     Vec(const Vec<Num_T> &v);
00220     Vec(const Vec<Num_T> &v, const Factory &f);
00222     Vec(const char *values, const Factory &f = DEFAULT_FACTORY);
00224     Vec(const std::string &values, const Factory &f = DEFAULT_FACTORY);
00226     Vec(const Num_T *c_array, int size, const Factory &f = DEFAULT_FACTORY);
00227 
00229     ~Vec();
00230 
00232     int length() const { return datasize; }
00234     int size() const { return datasize; }
00235 
00237     void set_size(int size, const bool copy=false);
00239     void set_length(int size, const bool copy=false) { set_size(size, copy); }
00241     void zeros();
00243     void clear() { zeros(); }
00245     void ones();
00247     void set(const char *str);
00249     void set(const std::string &str);
00250 
00252     const Num_T &operator[](int i) const;
00254     const Num_T &operator()(int i) const;
00256     Num_T &operator[](int i);
00258     Num_T &operator()(int i);
00260     const Vec<Num_T> operator()(int i1, int i2) const;
00262     const Vec<Num_T> operator()(const Vec<int> &indexlist) const;
00263 
00265     const Num_T &get(int i) const;
00267     const Vec<Num_T> get(int i1, int i2) const;
00269     void set(int i, const Num_T &v);
00270 
00272     Mat<Num_T> transpose() const;
00274     Mat<Num_T> T() const { return this->transpose(); }
00276     Mat<Num_T> hermitian_transpose() const;
00278     Mat<Num_T> H() const { return this->hermitian_transpose(); }
00279 
00281     Vec<Num_T>& operator+=(const Vec<Num_T> &v);
00283     Vec<Num_T>& operator+=(const Num_T t);
00285     friend const Vec<Num_T> operator+<>(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
00287     friend const Vec<Num_T> operator+<>(const Vec<Num_T> &v, const Num_T t);
00289     friend const Vec<Num_T> operator+<>(const Num_T t, const Vec<Num_T> &v);
00290 
00292     Vec<Num_T>& operator-=(const Vec<Num_T> &v);
00294     Vec<Num_T>& operator-=(const Num_T t);
00296     friend const Vec<Num_T> operator-<>(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
00298     friend const Vec<Num_T> operator-<>(const Vec<Num_T> &v, const Num_T t);
00300     friend const Vec<Num_T> operator-<>(const Num_T t, const Vec<Num_T> &v);
00302     friend const Vec<Num_T> operator-<>(const Vec<Num_T> &v);
00303 
00305     Vec<Num_T>& operator*=(const Num_T t);
00307     friend Num_T operator*<>(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
00309     friend Num_T dot <>(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
00311     friend const Mat<Num_T> outer_product <>(const Vec<Num_T> &v1,
00312                const Vec<Num_T> &v2,
00313                bool hermitian);
00315     friend const Vec<Num_T> operator*<>(const Vec<Num_T> &v, const Num_T t);
00317     friend const Vec<Num_T> operator*<>(const Num_T t, const Vec<Num_T> &v);
00318 
00320     friend const Vec<Num_T> elem_mult <>(const Vec<Num_T> &a, const Vec<Num_T> &b);
00322     friend const Vec<Num_T> elem_mult <>(const Vec<Num_T> &a, const Vec<Num_T> &b, const Vec<Num_T> &c);
00324     friend const Vec<Num_T> elem_mult <>(const Vec<Num_T> &a, const Vec<Num_T> &b, const Vec<Num_T> &c, const Vec<Num_T> &d);
00325 
00327     friend void elem_mult_out <>(const Vec<Num_T> &a, const Vec<Num_T> &b, Vec<Num_T> &out);
00329     friend void elem_mult_out <>(const Vec<Num_T> &a, const Vec<Num_T> &b, const Vec<Num_T> &c, Vec<Num_T> &out);
00331     friend void elem_mult_out <>(const Vec<Num_T> &a, const Vec<Num_T> &b, const Vec<Num_T> &c, const Vec<Num_T> &d, Vec<Num_T> &out);
00332 
00334     friend void elem_mult_inplace <>(const Vec<Num_T> &a, Vec<Num_T> &b);
00336     friend Num_T elem_mult_sum <>(const Vec<Num_T> &a, const Vec<Num_T> &b);
00337 
00339     Vec<Num_T>& operator/=(const Num_T t);
00341     Vec<Num_T>& operator/=(const Vec<Num_T> &v);
00342 
00344     friend const Vec<Num_T> operator/<>(const Vec<Num_T> &v, const Num_T t);
00346     friend const Vec<Num_T> operator/<>(const Num_T t, const Vec<Num_T> &v);
00347 
00349     friend const Vec<Num_T> elem_div <>(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
00351     friend const Vec<Num_T> elem_div <>(const Num_T t, const Vec<Num_T> &v);
00353     friend void elem_div_out <>(const Vec<Num_T> &v1, const Vec<Num_T> &v2, Vec<Num_T> &out);
00355     friend Num_T elem_div_sum <>(const Vec<Num_T> &a, const Vec<Num_T> &b);
00356 
00358     Vec<Num_T> get(const Vec<bin> &binlist) const;
00360     Vec<Num_T> right(int nr) const;
00362     Vec<Num_T> left(int nr) const;
00364     Vec<Num_T> mid(int start, int nr) const;
00366     Vec<Num_T> split(int pos);
00368     void shift_right(const Num_T In, int n=1);
00370     void shift_right(const Vec<Num_T> &In);
00372     void shift_left(const Num_T In, int n=1);
00374     void shift_left(const Vec<Num_T> &In);
00375 
00377     friend const Vec<Num_T> concat<>(const Vec<Num_T> &v, const Num_T a);
00379     friend const Vec<Num_T> concat<>(const Num_T a, const Vec<Num_T> &v);
00381     friend const Vec<Num_T> concat<>(const Vec<Num_T> &v1,const Vec<Num_T> &v2);
00383     friend const Vec<Num_T> concat<>(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
00384              const Vec<Num_T> &v3);
00386     friend const Vec<Num_T> concat<>(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
00387              const Vec<Num_T> &v3, const Vec<Num_T> &v4);
00389     friend const Vec<Num_T> concat<>(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
00390              const Vec<Num_T> &v3, const Vec<Num_T> &v4,
00391              const Vec<Num_T> &v5);
00392 
00394     void set_subvector(int i1, int i2, const Vec<Num_T> &v);
00396     void set_subvector(int i, const Vec<Num_T> &v);
00398     void set_subvector(int i1, int i2, const Num_T t);
00400     void replace_mid(int pos, const Vec<Num_T> &v);
00402     void del(int index);
00404     void del(int i1, int i2);
00406     void ins(int index, const Num_T in);
00408     void ins(int index, const Vec<Num_T> &in);
00409 
00411     Vec<Num_T>& operator=(const Num_T t);
00413     Vec<Num_T>& operator=(const Vec<Num_T> &v);
00415     Vec<Num_T>& operator=(const Mat<Num_T> &m);
00417     Vec<Num_T>& operator=(const char *values);
00418 
00420     Vec<bin> operator==(const Num_T value) const;
00422     Vec<bin> operator!=(const Num_T value) const;
00424     Vec<bin> operator<(const Num_T value) const;
00426     Vec<bin> operator<=(const Num_T value) const;
00428     Vec<bin> operator>(const Num_T value) const;
00430     Vec<bin> operator>=(const Num_T value) const;
00431 
00433     bool operator==(const Vec<Num_T> &v) const;
00435     bool operator!=(const Vec<Num_T> &v) const;
00436 
00438     Num_T &_elem(int i) { return data[i]; }
00440     const Num_T &_elem(int i) const { return data[i]; }
00441 
00443     Num_T *_data() { return data; }
00445     const Num_T *_data() const { return data; }
00446 
00447   protected:
00449     void alloc(int size);
00451     void free();
00452 
00454     int datasize;
00456     Num_T *data;
00458     const Factory &factory;
00459   private:
00460     // This function is used in set() methods to replace commas with spaces
00461     std::string replace_commas(const std::string &str);
00462   };
00463 
00464   //-----------------------------------------------------------------------------------
00465   // Type definitions of vec, cvec, ivec, svec, and bvec
00466   //-----------------------------------------------------------------------------------
00467 
00472   typedef Vec<double> vec;
00473 
00478   typedef Vec<std::complex<double> > cvec;
00479 
00484   typedef Vec<int> ivec;
00485 
00490   typedef Vec<short int> svec;
00491 
00496   typedef Vec<bin> bvec;
00497 
00498 } //namespace itpp
00499 
00500 
00501 #include <itpp/base/mat.h>
00502 
00503 namespace itpp {
00504 
00505   //-----------------------------------------------------------------------------------
00506   // Declaration of input and output streams for Vec
00507   //-----------------------------------------------------------------------------------
00508 
00513   template<class Num_T>
00514   std::ostream &operator<<(std::ostream &os, const Vec<Num_T> &v);
00515 
00527   template<class Num_T>
00528   std::istream &operator>>(std::istream &is, Vec<Num_T> &v);
00529 
00530   //-----------------------------------------------------------------------------------
00531   // Implementation of templated Vec members and friends
00532   //-----------------------------------------------------------------------------------
00533 
00534   template<class Num_T> inline
00535   void Vec<Num_T>::alloc(int size)
00536   {
00537     if (size > 0) {
00538       create_elements(data, size, factory);
00539       datasize = size;
00540     }
00541     else {
00542       data = 0;
00543       datasize = 0;
00544     }
00545   }
00546 
00547   template<class Num_T> inline
00548   void Vec<Num_T>::free()
00549   {
00550     destroy_elements(data, datasize);
00551     datasize = 0;
00552   }
00553 
00554 
00555   template<class Num_T> inline
00556   Vec<Num_T>::Vec(const Factory &f) : datasize(0), data(0), factory(f) {}
00557 
00558   template<class Num_T> inline
00559   Vec<Num_T>::Vec(int size, const Factory &f) : datasize(0), data(0), factory(f)
00560   {
00561     it_assert_debug(size>=0, "Negative size in Vec::Vec(int)");
00562     alloc(size);
00563   }
00564 
00565   template<class Num_T> inline
00566   Vec<Num_T>::Vec(const Vec<Num_T> &v) : datasize(0), data(0), factory(v.factory)
00567   {
00568     alloc(v.datasize);
00569     copy_vector(datasize, v.data, data);
00570   }
00571 
00572   template<class Num_T> inline
00573   Vec<Num_T>::Vec(const Vec<Num_T> &v, const Factory &f) : datasize(0), data(0), factory(f)
00574   {
00575     alloc(v.datasize);
00576     copy_vector(datasize, v.data, data);
00577   }
00578 
00579   template<class Num_T> inline
00580   Vec<Num_T>::Vec(const char *values, const Factory &f) : datasize(0), data(0), factory(f)
00581   {
00582     set(values);
00583   }
00584 
00585   template<class Num_T> inline
00586   Vec<Num_T>::Vec(const std::string &values, const Factory &f) : datasize(0), data(0), factory(f)
00587   {
00588     set(values);
00589   }
00590 
00591   template<class Num_T> inline
00592   Vec<Num_T>::Vec(const Num_T *c_array, int size, const Factory &f) : datasize(0), data(0), factory(f)
00593   {
00594     alloc(size);
00595     copy_vector(size, c_array, data);
00596   }
00597 
00598   template<class Num_T> inline
00599   Vec<Num_T>::~Vec()
00600   {
00601     free();
00602   }
00603 
00604   template<class Num_T>
00605   void Vec<Num_T>::set_size(int size, bool copy)
00606   {
00607     it_assert_debug(size >= 0, "Vec::set_size(): New size must not be negative");
00608     if (datasize == size)
00609       return;
00610     if (copy) {
00611       // create a temporary pointer to the allocated data
00612       Num_T* tmp = data;
00613       // store the current number of elements
00614       int old_datasize = datasize;
00615       // check how many elements we need to copy
00616       int min = datasize < size ? datasize : size;
00617       // allocate new memory
00618       alloc(size);
00619       // copy old elements into a new memory region
00620       copy_vector(min, tmp, data);
00621       // initialize the rest of resized vector
00622       for (int i = min; i < size; ++i)
00623   data[i] = Num_T(0);
00624       // delete old elements
00625       destroy_elements(tmp, old_datasize);
00626     }
00627     else {
00628       free();
00629       alloc(size);
00630     }
00631   }
00632 
00633   template<class Num_T> inline
00634   const Num_T& Vec<Num_T>::operator[](int i) const
00635   {
00636     it_assert_debug((i >= 0) && (i < datasize), "Vec::operator[]: Index out of range");
00637     return data[i];
00638   }
00639 
00640   template<class Num_T> inline
00641   const Num_T& Vec<Num_T>::operator()(int i) const
00642   {
00643     it_assert_debug((i >= 0) && (i < datasize), "Vec::operator(): Index out of range");
00644     return data[i];
00645   }
00646 
00647   template<class Num_T> inline
00648   Num_T& Vec<Num_T>::operator[](int i)
00649   {
00650     it_assert_debug((i >= 0) && (i < datasize), "Vec::operator[]: Index out of range");
00651     return data[i];
00652   }
00653 
00654   template<class Num_T> inline
00655   Num_T& Vec<Num_T>::operator()(int i)
00656   {
00657     it_assert_debug((i >= 0) && (i < datasize), "Vec::operator(): Index out of range");
00658     return data[i];
00659   }
00660 
00661   template<class Num_T> inline
00662   const Vec<Num_T> Vec<Num_T>::operator()(int i1, int i2) const
00663   {
00664     int ii1=i1, ii2=i2;
00665 
00666     if (ii1 == -1) ii1 = datasize-1;
00667     if (ii2 == -1) ii2 = datasize-1;
00668 
00669     it_assert_debug(ii1>=0 && ii2>=0 && ii1<datasize && ii2<datasize, "Vec::operator()(i1,i2): indicies out of range");
00670     it_assert_debug(ii2>=ii1, "Vec::op(i1,i2): i2 >= i1 necessary");
00671 
00672     Vec<Num_T> s(ii2-ii1+1);
00673     copy_vector(s.datasize, data+ii1, s.data);
00674 
00675     return s;
00676   }
00677 
00678   template<class Num_T>
00679   const Vec<Num_T> Vec<Num_T>::operator()(const Vec<int> &indexlist) const
00680   {
00681     Vec<Num_T> temp(indexlist.length());
00682     for (int i=0;i<indexlist.length();i++) {
00683       it_assert((indexlist(i)>=0) && (indexlist(i) < datasize), "Vec::operator()(ivec &): index outside range");
00684       temp(i)=data[indexlist(i)];
00685     }
00686     return temp;
00687   }
00688 
00689 
00690   template<class Num_T> inline
00691   const Num_T& Vec<Num_T>::get(int i) const
00692   {
00693     it_assert_debug((i >= 0) && (i < datasize), "method get()");
00694     return data[i];
00695   }
00696 
00697   template<class Num_T> inline
00698   const Vec<Num_T> Vec<Num_T>::get(int i1, int i2) const
00699   {
00700     return (*this)(i1, i2);
00701   }
00702 
00703 
00704   template<class Num_T> inline
00705   void Vec<Num_T>::zeros()
00706   {
00707     for (int i = 0; i < datasize; i++)
00708       data[i] = Num_T(0);
00709   }
00710 
00711   template<class Num_T> inline
00712   void Vec<Num_T>::ones()
00713   {
00714     for (int i = 0; i < datasize; i++)
00715       data[i] = Num_T(1);
00716   }
00717 
00718   template<class Num_T> inline
00719   void Vec<Num_T>::set(int i, const Num_T &v)
00720   {
00721     it_assert_debug((i >= 0) && (i < datasize), "method set()");
00722     data[i] = v;
00723   }
00724 
00726   template<>
00727   void Vec<double>::set(const std::string &str);
00728   template<>
00729   void Vec<std::complex<double> >::set(const std::string &str);
00730   template<>
00731   void Vec<int>::set(const std::string &str);
00732   template<>
00733   void Vec<short int>::set(const std::string &str);
00734   template<>
00735   void Vec<bin>::set(const std::string &str);
00737 
00738   template<class Num_T>
00739   void Vec<Num_T>::set(const std::string &str)
00740   {
00741     it_error("Vec::set(): Only `double', `complex<double>', `int', "
00742        "`short int' and `bin' types supported");
00743   }
00744 
00745   template<class Num_T>
00746   void Vec<Num_T>::set(const char *str)
00747   {
00748     set(std::string(str));
00749   }
00750 
00751 
00752   template<class Num_T>
00753   Mat<Num_T> Vec<Num_T>::transpose() const
00754   {
00755     Mat<Num_T> temp(1, datasize);
00756     copy_vector(datasize, data, temp._data());
00757     return temp;
00758   }
00759 
00760   template<>
00761   Mat<std::complex<double> > Vec<std::complex<double> >::hermitian_transpose() const;
00762 
00763   template<class Num_T>
00764   Mat<Num_T> Vec<Num_T>::hermitian_transpose() const
00765   {
00766     Mat<Num_T> temp(1, datasize);
00767     copy_vector(datasize, data, temp._data());
00768     return temp;
00769   }
00770 
00771   template<class Num_T> inline
00772   Vec<Num_T>& Vec<Num_T>::operator+=(const Vec<Num_T> &v)
00773   {
00774     if (datasize == 0) { // if not assigned a size.
00775       if (this != &v) { // check for self addition
00776   alloc(v.datasize);
00777   copy_vector(datasize, v.data, data);
00778       }
00779     } else {
00780       it_assert_debug(datasize == v.datasize, "Vec::operator+=: Wrong sizes");
00781       for (int i = 0; i < datasize; i++)
00782   data[i] += v.data[i];
00783     }
00784     return *this;
00785   }
00786 
00787   template<class Num_T> inline
00788   Vec<Num_T>& Vec<Num_T>::operator+=(const Num_T t)
00789   {
00790     for (int i=0;i<datasize;i++)
00791       data[i]+=t;
00792     return *this;
00793   }
00794 
00795   template<class Num_T> inline
00796   const Vec<Num_T> operator+(const Vec<Num_T> &v1, const Vec<Num_T> &v2)
00797   {
00798     int i;
00799     Vec<Num_T> r(v1.datasize);
00800 
00801     it_assert_debug(v1.datasize==v2.datasize, "Vec::operator+: wrong sizes");
00802     for (i=0; i<v1.datasize; i++)
00803       r.data[i] = v1.data[i] + v2.data[i];
00804 
00805     return r;
00806   }
00807 
00808   template<class Num_T> inline
00809   const Vec<Num_T> operator+(const Vec<Num_T> &v, const Num_T t)
00810   {
00811     int i;
00812     Vec<Num_T> r(v.datasize);
00813 
00814     for (i=0; i<v.datasize; i++)
00815       r.data[i] = v.data[i] + t;
00816 
00817     return r;
00818   }
00819 
00820   template<class Num_T> inline
00821   const Vec<Num_T> operator+(const Num_T t, const Vec<Num_T> &v)
00822   {
00823     int i;
00824     Vec<Num_T> r(v.datasize);
00825 
00826     for (i=0; i<v.datasize; i++)
00827       r.data[i] = t + v.data[i];
00828 
00829     return r;
00830   }
00831 
00832   template<class Num_T> inline
00833   Vec<Num_T>& Vec<Num_T>::operator-=(const Vec<Num_T> &v)
00834   {
00835     if (datasize == 0) { // if not assigned a size.
00836       if (this != &v) { // check for self decrementation
00837   alloc(v.datasize);
00838   for (int i = 0; i < v.datasize; i++)
00839     data[i] = -v.data[i];
00840       }
00841     } else {
00842       it_assert_debug(datasize == v.datasize, "Vec::operator-=: Wrong sizes");
00843       for (int i = 0; i < datasize; i++)
00844   data[i] -= v.data[i];
00845     }
00846     return *this;
00847   }
00848 
00849   template<class Num_T> inline
00850   Vec<Num_T>& Vec<Num_T>::operator-=(const Num_T t)
00851   {
00852     for (int i=0;i<datasize;i++)
00853       data[i]-=t;
00854     return *this;
00855   }
00856 
00857   template<class Num_T> inline
00858   const Vec<Num_T> operator-(const Vec<Num_T> &v1, const Vec<Num_T> &v2)
00859   {
00860     int i;
00861     Vec<Num_T> r(v1.datasize);
00862 
00863     it_assert_debug(v1.datasize==v2.datasize, "Vec::operator-: wrong sizes");
00864     for (i=0; i<v1.datasize; i++)
00865       r.data[i] = v1.data[i] - v2.data[i];
00866 
00867     return r;
00868   }
00869 
00870   template<class Num_T> inline
00871   const Vec<Num_T> operator-(const Vec<Num_T> &v, const Num_T t)
00872   {
00873     int i;
00874     Vec<Num_T> r(v.datasize);
00875 
00876     for (i=0; i<v.datasize; i++)
00877       r.data[i] = v.data[i] - t;
00878 
00879     return r;
00880   }
00881 
00882   template<class Num_T> inline
00883   const Vec<Num_T> operator-(const Num_T t, const Vec<Num_T> &v)
00884   {
00885     int i;
00886     Vec<Num_T> r(v.datasize);
00887 
00888     for (i=0; i<v.datasize; i++)
00889       r.data[i] = t - v.data[i];
00890 
00891     return r;
00892   }
00893 
00894   template<class Num_T> inline
00895   const Vec<Num_T> operator-(const Vec<Num_T> &v)
00896   {
00897     int i;
00898     Vec<Num_T> r(v.datasize);
00899 
00900     for (i=0; i<v.datasize; i++)
00901       r.data[i] = -v.data[i];
00902 
00903     return r;
00904   }
00905 
00906   template<class Num_T> inline
00907   Vec<Num_T>& Vec<Num_T>::operator*=(const Num_T t)
00908   {
00909     scal_vector(datasize, t, data);
00910     return *this;
00911   }
00912 
00913 #if defined(HAVE_BLAS)
00914   template<> inline
00915   double dot(const vec &v1, const vec &v2)
00916   {
00917     it_assert_debug(v1.datasize == v2.datasize, "vec::dot: wrong sizes");
00918     int incr = 1;
00919     return blas::ddot_(&v1.datasize, v1.data, &incr, v2.data, &incr);
00920   }
00921 
00922   template<> inline
00923   std::complex<double> dot(const cvec &v1, const cvec &v2)
00924   {
00925     it_assert_debug(v1.datasize == v2.datasize, "cvec::dot: wrong sizes");
00926     int incr = 1;
00927     std::complex<double> output;
00928     blas::zdotusub_(&output, &v1.datasize, v1.data, &incr, v2.data, &incr);
00929     return output;
00930   }
00931 #endif // HAVE_BLAS
00932 
00933   template<class Num_T> inline
00934   Num_T dot(const Vec<Num_T> &v1, const Vec<Num_T> &v2)
00935   {
00936     int i;
00937     Num_T r=Num_T(0);
00938 
00939     it_assert_debug(v1.datasize==v2.datasize, "Vec::dot: wrong sizes");
00940     for (i=0; i<v1.datasize; i++)
00941       r += v1.data[i] * v2.data[i];
00942 
00943     return r;
00944   }
00945 
00946   template<class Num_T> inline
00947   Num_T operator*(const Vec<Num_T> &v1, const Vec<Num_T> &v2)
00948   {
00949     return dot(v1, v2);
00950   }
00951 
00952 
00953 #if defined(HAVE_BLAS)
00954   template<> inline
00955   const mat outer_product(const vec &v1, const vec &v2, bool hermitian)
00956   {
00957     it_assert_debug((v1.datasize > 0) && (v2.datasize > 0),
00958         "Vec::outer_product():: Input vector of zero size");
00959 
00960     mat out(v1.datasize, v2.datasize);
00961     out.zeros();
00962     double alpha = 1.0;
00963     int incr = 1;
00964     blas::dger_(&v1.datasize, &v2.datasize, &alpha, v1.data, &incr,
00965     v2.data, &incr, out._data(), &v1.datasize);
00966     return out;
00967   }
00968 
00969   template<> inline
00970   const cmat outer_product(const cvec &v1, const cvec &v2, bool hermitian)
00971   {
00972     it_assert_debug((v1.datasize > 0) && (v2.datasize > 0),
00973         "Vec::outer_product():: Input vector of zero size");
00974 
00975     cmat out(v1.datasize, v2.datasize);
00976     out.zeros();
00977     std::complex<double> alpha(1.0);
00978     int incr = 1;
00979     if (hermitian) {
00980       blas::zgerc_(&v1.datasize, &v2.datasize, &alpha, v1.data, &incr,
00981        v2.data, &incr, out._data(), &v1.datasize);
00982     }
00983     else {
00984       blas::zgeru_(&v1.datasize, &v2.datasize, &alpha, v1.data, &incr,
00985        v2.data, &incr, out._data(), &v1.datasize);
00986     }
00987     return out;
00988   }
00989 #else
00991   template<> inline
00992   const cmat outer_product(const cvec &v1, const cvec &v2, bool hermitian)
00993   {
00994     it_assert_debug((v1.datasize > 0) && (v2.datasize > 0),
00995         "Vec::outer_product():: Input vector of zero size");
00996 
00997     cmat out(v1.datasize, v2.datasize);
00998     if (hermitian) {
00999       for (int i = 0; i < v1.datasize; ++i) {
01000   for (int j = 0; j < v2.datasize; ++j) {
01001     out(i, j) = v1.data[i] * conj(v2.data[j]);
01002   }
01003       }
01004     }
01005     else {
01006       for (int i = 0; i < v1.datasize; ++i) {
01007   for (int j = 0; j < v2.datasize; ++j) {
01008     out(i, j) = v1.data[i] * v2.data[j];
01009   }
01010       }
01011     }
01012     return out;
01013   }
01014 #endif // HAVE_BLAS
01015 
01016   template<class Num_T> inline
01017   const Mat<Num_T> outer_product(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
01018          bool hermitian)
01019   {
01020     int i, j;
01021 
01022     it_assert_debug((v1.datasize > 0) && (v2.datasize > 0),
01023         "Vec::outer_product:: Input vector of zero size");
01024 
01025     Mat<Num_T> r(v1.datasize, v2.datasize);
01026     for (i=0; i<v1.datasize; i++) {
01027       for (j=0; j<v2.datasize; j++) {
01028   r(i,j) = v1.data[i] * v2.data[j];
01029       }
01030     }
01031 
01032     return r;
01033   }
01034 
01035   template<class Num_T> inline
01036   const Vec<Num_T> operator*(const Vec<Num_T> &v, const Num_T t)
01037   {
01038     int i;
01039     Vec<Num_T> r(v.datasize);
01040     for (i=0; i<v.datasize; i++)
01041       r.data[i] = v.data[i] * t;
01042 
01043     return r;
01044   }
01045 
01046   template<class Num_T> inline
01047   const Vec<Num_T> operator*(const Num_T t, const Vec<Num_T> &v)
01048   {
01049     return operator*(v, t);
01050   }
01051 
01052   template<class Num_T> inline
01053   const Vec<Num_T> elem_mult(const Vec<Num_T> &a, const Vec<Num_T> &b)
01054   {
01055     Vec<Num_T> out;
01056     elem_mult_out(a,b,out);
01057     return out;
01058   }
01059 
01060   template<class Num_T> inline
01061   const Vec<Num_T> elem_mult(const Vec<Num_T> &a, const Vec<Num_T> &b, const Vec<Num_T> &c)
01062   {
01063     Vec<Num_T> out;
01064     elem_mult_out(a,b,c,out);
01065     return out;
01066   }
01067 
01068   template<class Num_T> inline
01069   const Vec<Num_T> elem_mult(const Vec<Num_T> &a, const Vec<Num_T> &b, const Vec<Num_T> &c, const Vec<Num_T> &d)
01070   {
01071     Vec<Num_T> out;
01072     elem_mult_out(a,b,c,d,out);
01073     return out;
01074   }
01075 
01076   template<class Num_T> inline
01077   void elem_mult_out(const Vec<Num_T> &a, const Vec<Num_T> &b, Vec<Num_T> &out)
01078   {
01079     it_assert_debug(a.datasize==b.datasize, "Vec::elem_mult_out: wrong sizes");
01080 
01081     if(out.datasize != a.datasize)
01082       out.set_size(a.size());
01083 
01084     for(int i=0; i<a.datasize; i++)
01085       out.data[i] = a.data[i] * b.data[i];
01086   }
01087 
01088   template<class Num_T> inline
01089   void elem_mult_out(const Vec<Num_T> &a, const Vec<Num_T> &b, const Vec<Num_T> &c, Vec<Num_T> &out)
01090   {
01091     it_assert_debug(a.datasize==b.datasize==c.datasize, "Vec::elem_mult_out: wrong sizes");
01092 
01093     if(out.datasize != a.datasize)
01094       out.set_size(a.size());
01095 
01096     for(int i=0; i<a.datasize; i++)
01097       out.data[i] = a.data[i] * b.data[i] * c.data[i];
01098   }
01099 
01100   template<class Num_T> inline
01101   void elem_mult_out(const Vec<Num_T> &a, const Vec<Num_T> &b, const Vec<Num_T> &c, const Vec<Num_T> &d, Vec<Num_T> &out)
01102   {
01103     it_assert_debug(a.datasize==b.datasize==c.datasize==d.datasize, "Vec::elem_mult_out: wrong sizes");
01104 
01105     if(out.datasize != a.datasize)
01106       out.set_size(a.size());
01107 
01108     for(int i=0; i<a.datasize; i++)
01109       out.data[i] = a.data[i] * b.data[i] * c.data[i] * d.data[i];
01110   }
01111 
01112   template<class Num_T> inline
01113   void elem_mult_inplace(const Vec<Num_T> &a, Vec<Num_T> &b)
01114   {
01115     it_assert_debug(a.datasize==b.datasize, "Vec::elem_mult_inplace: wrong sizes");
01116 
01117     for(int i=0; i<a.datasize; i++)
01118       b.data[i] *= a.data[i];
01119   }
01120 
01121   template<class Num_T> inline
01122   Num_T elem_mult_sum(const Vec<Num_T> &a, const Vec<Num_T> &b)
01123   {
01124     it_assert_debug(a.datasize==b.datasize, "Vec::elem_mult_sum: wrong sizes");
01125 
01126     Num_T acc = 0;
01127 
01128     for(int i=0; i<a.datasize; i++)
01129       acc += a.data[i] * b.data[i];
01130 
01131     return acc;
01132   }
01133 
01134   template<class Num_T> inline
01135   const Vec<Num_T> operator/(const Vec<Num_T> &v, const Num_T t)
01136   {
01137     int i;
01138     Vec<Num_T> r(v.datasize);
01139 
01140     for (i=0; i<v.datasize; i++)
01141       r.data[i] = v.data[i] / t;
01142 
01143     return r;
01144   }
01145 
01146   template<class Num_T> inline
01147   const Vec<Num_T> operator/(const Num_T t, const Vec<Num_T> &v)
01148   {
01149     int i;
01150     Vec<Num_T> r(v.datasize);
01151 
01152     for (i=0; i<v.datasize; i++)
01153       r.data[i] = t / v.data[i];
01154 
01155     return r;
01156   }
01157 
01158   template<class Num_T> inline
01159   Vec<Num_T>& Vec<Num_T>::operator/=(const Num_T t)
01160   {
01161     for (int i = 0; i < datasize; ++i) {
01162       data[i] /= t;
01163     }
01164     return *this;
01165   }
01166 
01167   template<class Num_T> inline
01168   Vec<Num_T>& Vec<Num_T>::operator/=(const Vec<Num_T> &v)
01169   {
01170     it_assert_debug(datasize == v.datasize, "Vec::operator/=(): wrong sizes");
01171     for (int i = 0; i < datasize; ++i) {
01172       data[i] /= v.data[i];
01173     }
01174     return *this;
01175   }
01176 
01177   template<class Num_T> inline
01178   const Vec<Num_T> elem_div(const Vec<Num_T> &a, const Vec<Num_T> &b)
01179   {
01180     Vec<Num_T> out;
01181     elem_div_out(a,b,out);
01182     return out;
01183   }
01184 
01185   template<class Num_T> inline
01186   const Vec<Num_T> elem_div(const Num_T t, const Vec<Num_T> &v)
01187   {
01188     int i;
01189     Vec<Num_T> r(v.datasize);
01190 
01191     for (i=0; i<v.datasize; i++)
01192       r.data[i] = t / v.data[i];
01193 
01194     return r;
01195   }
01196 
01197   template<class Num_T> inline
01198   void elem_div_out(const Vec<Num_T> &a, const Vec<Num_T> &b, Vec<Num_T> &out)
01199   {
01200     it_assert_debug(a.datasize==b.datasize, "Vecelem_div_out: wrong sizes");
01201 
01202     if(out.datasize != a.datasize)
01203       out.set_size(a.size());
01204 
01205     for(int i=0; i<a.datasize; i++)
01206       out.data[i] = a.data[i] / b.data[i];
01207   }
01208 
01209   template<class Num_T> inline
01210   Num_T elem_div_sum(const Vec<Num_T> &a, const Vec<Num_T> &b)
01211   {
01212     it_assert_debug(a.datasize==b.datasize, "Vec::elem_div_sum: wrong sizes");
01213 
01214     Num_T acc = 0;
01215 
01216     for(int i=0; i<a.datasize; i++)
01217       acc += a.data[i] / b.data[i];
01218 
01219     return acc;
01220   }
01221 
01222   template<class Num_T>
01223   Vec<Num_T> Vec<Num_T>::get(const Vec<bin> &binlist) const
01224   {
01225     int size = binlist.size();
01226     it_assert_debug(datasize == size, "Vec::get(bvec &): wrong sizes");
01227     Vec<Num_T> temp(size);
01228     int j = 0;
01229     for (int i = 0; i < size; ++i) {
01230       if (binlist(i) == bin(1)) {
01231   temp(j) = data[i];
01232   j++;
01233       }
01234     }
01235     temp.set_size(j, true);
01236     return temp;
01237   }
01238 
01239   template<class Num_T> inline
01240   Vec<Num_T> Vec<Num_T>::right(int nr) const
01241   {
01242     it_assert_debug(nr <= datasize, "Vec::right(): index out of range");
01243     Vec<Num_T> temp(nr);
01244     if (nr > 0) {
01245       copy_vector(nr, &data[datasize-nr], temp.data);
01246     }
01247     return temp;
01248   }
01249 
01250   template<class Num_T> inline
01251   Vec<Num_T> Vec<Num_T>::left(int nr) const
01252   {
01253     it_assert_debug(nr <= datasize, "Vec::left(): index out of range");
01254     Vec<Num_T> temp(nr);
01255     if (nr > 0) {
01256       copy_vector(nr, data, temp.data);
01257     }
01258     return temp;
01259   }
01260 
01261   template<class Num_T> inline
01262   Vec<Num_T> Vec<Num_T>::mid(int start, int nr) const
01263   {
01264     it_assert_debug((start >= 0) && ((start+nr) <= datasize),
01265         "Vec::mid(): indexing out of range");
01266     Vec<Num_T> temp(nr);
01267     if (nr > 0) {
01268       copy_vector(nr, &data[start], temp.data);
01269     }
01270     return temp;
01271   }
01272 
01273   template<class Num_T>
01274   Vec<Num_T> Vec<Num_T>::split(int pos)
01275   {
01276     it_assert_debug((pos >= 0) && (pos <= datasize), "Vec::split(): index out of range");
01277     Vec<Num_T> temp1(pos);
01278     Vec<Num_T> temp2(datasize-pos);
01279     copy_vector(pos, data, temp1.data);
01280     copy_vector(datasize-pos, &data[pos], temp2.data);
01281     (*this) = temp2;
01282     return temp1;
01283   }
01284 
01285   template<class Num_T>
01286   void Vec<Num_T>::shift_right(const Num_T In, int n)
01287   {
01288     int i=datasize;
01289 
01290     it_assert_debug(n>=0, "Vec::shift_right: index out of range");
01291     while (--i >= n)
01292       data[i] = data[i-n];
01293     while (i >= 0)
01294       data[i--] = In;
01295   }
01296 
01297   template<class Num_T>
01298   void Vec<Num_T>::shift_right(const Vec<Num_T> &In)
01299   {
01300     int i;
01301 
01302     for (i=datasize-1; i>=In.datasize; i--)
01303       data[i]=data[i-In.datasize];
01304     for (i=0; i<In.datasize; i++)
01305       data[i]=In[i];
01306   }
01307 
01308   template<class Num_T>
01309   void Vec<Num_T>::shift_left(const Num_T In, int n)
01310   {
01311     int i;
01312 
01313     it_assert_debug(n>=0, "Vec::shift_left: index out of range");
01314     for (i=0; i<datasize-n; i++)
01315       data[i] = data[i+n];
01316     while (i < datasize)
01317       data[i++] = In;
01318   }
01319 
01320   template<class Num_T>
01321   void Vec<Num_T>::shift_left(const Vec<Num_T> &In)
01322   {
01323     int i;
01324 
01325     for (i=0; i<datasize-In.datasize; i++)
01326       data[i]=data[i+In.datasize];
01327     for (i=datasize-In.datasize; i<datasize; i++)
01328       data[i]=In[i-datasize+In.datasize];
01329   }
01330 
01331   template<class Num_T>
01332   const Vec<Num_T> concat(const Vec<Num_T> &v, const Num_T a)
01333   {
01334     int size = v.size();
01335     Vec<Num_T> temp(size+1);
01336     copy_vector(size, v.data, temp.data);
01337     temp(size) = a;
01338     return temp;
01339   }
01340 
01341   template<class Num_T>
01342   const Vec<Num_T> concat(const Num_T a, const Vec<Num_T> &v)
01343   {
01344     int size = v.size();
01345     Vec<Num_T> temp(size+1);
01346     temp(0) = a;
01347     copy_vector(size, v.data, &temp.data[1]);
01348     return temp;
01349   }
01350 
01351   template<class Num_T>
01352   const Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2)
01353   {
01354     int size1 = v1.size();
01355     int size2 = v2.size();
01356     Vec<Num_T> temp(size1+size2);
01357     copy_vector(size1, v1.data, temp.data);
01358     copy_vector(size2, v2.data, &temp.data[size1]);
01359     return temp;
01360   }
01361 
01362   template<class Num_T>
01363   const Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
01364         const Vec<Num_T> &v3)
01365   {
01366     int size1 = v1.size();
01367     int size2 = v2.size();
01368     int size3 = v3.size();
01369     Vec<Num_T> temp(size1+size2+size3);
01370     copy_vector(size1, v1.data, temp.data);
01371     copy_vector(size2, v2.data, &temp.data[size1]);
01372     copy_vector(size3, v3.data, &temp.data[size1+size2]);
01373     return temp;
01374   }
01375 
01376   template<class Num_T>
01377   const Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
01378         const Vec<Num_T> &v3, const Vec<Num_T> &v4)
01379   {
01380     int size1 = v1.size();
01381     int size2 = v2.size();
01382     int size3 = v3.size();
01383     int size4 = v4.size();
01384     Vec<Num_T> temp(size1+size2+size3+size4);
01385     copy_vector(size1, v1.data, temp.data);
01386     copy_vector(size2, v2.data, &temp.data[size1]);
01387     copy_vector(size3, v3.data, &temp.data[size1+size2]);
01388     copy_vector(size4, v4.data, &temp.data[size1+size2+size3]);
01389     return temp;
01390   }
01391 
01392   template<class Num_T>
01393   const Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
01394         const Vec<Num_T> &v3, const Vec<Num_T> &v4,
01395         const Vec<Num_T> &v5)
01396   {
01397     int size1 = v1.size();
01398     int size2 = v2.size();
01399     int size3 = v3.size();
01400     int size4 = v4.size();
01401     int size5 = v5.size();
01402     Vec<Num_T> temp(size1+size2+size3+size4+size5);
01403     copy_vector(size1, v1.data, temp.data);
01404     copy_vector(size2, v2.data, &temp.data[size1]);
01405     copy_vector(size3, v3.data, &temp.data[size1+size2]);
01406     copy_vector(size4, v4.data, &temp.data[size1+size2+size3]);
01407     copy_vector(size5, v5.data, &temp.data[size1+size2+size3+size4]);
01408     return temp;
01409   }
01410 
01411   template<class Num_T> inline
01412   void Vec<Num_T>::set_subvector(int i1, int i2, const Vec<Num_T> &v)
01413   {
01414     if (i1 == -1) i1 = datasize-1;
01415     if (i2 == -1) i2 = datasize-1;
01416 
01417     it_assert_debug(i1>=0 && i2>=0 && i1<datasize && i2<datasize, "Vec::set_subvector(): indicies out of range");
01418     it_assert_debug(i2>=i1, "Vec::set_subvector(): i2 >= i1 necessary");
01419     it_assert_debug(i2-i1+1 == v.datasize, "Vec::set_subvector(): wrong sizes");
01420 
01421     copy_vector(v.datasize, v.data, data+i1);
01422   }
01423 
01424   template<class Num_T> inline
01425   void Vec<Num_T>:: set_subvector(int i, const Vec<Num_T> &v)
01426   {
01427     it_assert_debug(i>=0, "Vec::set_subvector(): index out of range");
01428     it_assert_debug(i+v.datasize <= datasize, "Vec::set_subvector(): too long input vector");
01429     copy_vector(v.datasize, v.data, data+i);
01430   }
01431 
01432   template<class Num_T>
01433   void Vec<Num_T>::set_subvector(int i1, int i2, const Num_T t)
01434   {
01435     if (i1 == -1) i1 = datasize-1;
01436     if (i2 == -1) i2 = datasize-1;
01437 
01438     it_assert_debug(i1>=0 && i2>=0 && i1<datasize && i2<datasize, "Vec::set_subvector(): indicies out of range");
01439     it_assert_debug(i2>=i1, "Vec::set_subvector(): i2 >= i1 necessary");
01440 
01441     for (int i=i1;i<=i2;i++)
01442       data[i] = t;
01443   }
01444 
01445   template<class Num_T>
01446   void Vec<Num_T>::replace_mid(int pos, const Vec<Num_T> &v)
01447   {
01448     it_assert_debug((pos>=0) && ((pos+v.length())<=datasize), "Vec::replace_mid: indexing out of range");
01449     copy_vector(v.datasize, v.data, &data[pos]);
01450   }
01451 
01452   template<class Num_T>
01453   void Vec<Num_T>::del(int index)
01454   {
01455     it_assert_debug((index >= 0) && (index < datasize),
01456         "Vec::del(): index out of range");
01457     Vec<Num_T> temp(*this);
01458     set_size(datasize-1, false);
01459     copy_vector(index, temp.data, data);
01460     copy_vector(datasize-index, &temp.data[index+1], &data[index]);
01461   }
01462 
01463   template<class Num_T>
01464   void  Vec<Num_T>::del(int i1, int i2)
01465   {
01466     it_assert_debug((i1 >= 0) && (i2 < datasize) && (i1 < i2),
01467         "Vec::del(): index out of range");
01468     Vec<Num_T> temp(*this);
01469     int new_size = datasize-(i2-i1+1);
01470     set_size(new_size, false);
01471     copy_vector(i1, temp.data, data);
01472     copy_vector(datasize-i1, &temp.data[i2+1], &data[i1]);
01473   }
01474 
01475   template<class Num_T>
01476   void Vec<Num_T>::ins(int index, const Num_T in)
01477   {
01478     it_assert_debug((index>=0) && (index<=datasize), "Vec::ins: index out of range");
01479     Vec<Num_T> Temp(*this);
01480 
01481     set_size(datasize+1, false);
01482     copy_vector(index, Temp.data, data);
01483     data[index]=in;
01484     copy_vector(Temp.datasize-index, Temp.data+index, data+index+1);
01485   }
01486 
01487   template<class Num_T>
01488   void Vec<Num_T>::ins(int index, const Vec<Num_T> &in)
01489   {
01490     it_assert_debug((index>=0) && (index<=datasize), "Vec::ins: index out of range");
01491     Vec<Num_T> Temp(*this);
01492 
01493     set_size(datasize+in.length(), false);
01494     copy_vector(index, Temp.data, data);
01495     copy_vector(in.size(), in.data, &data[index]);
01496     copy_vector(Temp.datasize-index, Temp.data+index, data+index+in.size());
01497   }
01498 
01499   template<class Num_T> inline
01500   Vec<Num_T>& Vec<Num_T>::operator=(const Num_T t)
01501   {
01502     for (int i=0;i<datasize;i++)
01503       data[i] = t;
01504     return *this;
01505   }
01506 
01507   template<class Num_T> inline
01508   Vec<Num_T>& Vec<Num_T>::operator=(const Vec<Num_T> &v)
01509   {
01510     if (this != &v) {
01511       set_size(v.datasize, false);
01512       copy_vector(datasize, v.data, data);
01513     }
01514     return *this;
01515   }
01516 
01517   template<class Num_T> inline
01518   Vec<Num_T>& Vec<Num_T>::operator=(const Mat<Num_T> &m)
01519   {
01520     it_assert_debug( (m.cols() == 1 && datasize == m.rows()) ||
01521     (m.rows() == 1 && datasize == m.cols()), "Vec::operator=(Mat<Num_T>): wrong size");
01522 
01523     if (m.cols() == 1) {
01524       set_size(m.rows(), false);
01525       copy_vector(m.rows(), m._data(), data);
01526     } else if (m.rows() == 1) {
01527       set_size(m.cols(), false);
01528       copy_vector(m.cols(), m._data(), m.rows(), data, 1);
01529     } else
01530       it_error("Vec::operator=(Mat<Num_T>): wrong size");
01531     return *this;
01532   }
01533 
01534   template<class Num_T> inline
01535   Vec<Num_T>& Vec<Num_T>::operator=(const char *values)
01536   {
01537     set(values);
01538     return *this;
01539   }
01540 
01541   template<>
01542   bvec Vec<std::complex<double> >::operator==(std::complex<double>) const;
01543 
01544   template<class Num_T>
01545   bvec Vec<Num_T>::operator==(const Num_T value) const
01546   {
01547     it_assert(datasize > 0, "Vec::operator==: vector must have size > 0");
01548     Vec<Num_T> invector(*this);
01549     bvec temp(invector.length());
01550 
01551     for (int i=0;i<invector.length();i++)
01552       temp(i)=(invector(i)==value);
01553 
01554     return temp;
01555   }
01556 
01557   template<>
01558   bvec Vec<std::complex<double> >::operator!=(std::complex<double>) const;
01559 
01560   template<class Num_T>
01561   bvec Vec<Num_T>::operator!=(const Num_T value) const
01562   {
01563     it_assert(datasize > 0, "Vec::operator!=: vector must have size > 0");
01564     Vec<Num_T> invector(*this);
01565     bvec temp(invector.length());
01566 
01567     for (int i=0;i<invector.length();i++)
01568       temp(i)=(invector(i)!=value);
01569 
01570     return temp;
01571   }
01572 
01573   template<>
01574   bvec Vec<std::complex<double> >::operator<(std::complex<double>) const;
01575 
01576   template<class Num_T>
01577   bvec Vec<Num_T>::operator<(const Num_T value) const
01578   {
01579     it_assert(datasize > 0, "Vec::operator<: vector must have size > 0");
01580     Vec<Num_T> invector(*this);
01581     bvec temp(invector.length());
01582 
01583     for (int i=0;i<invector.length();i++)
01584       temp(i)=(invector(i)<value);
01585 
01586     return temp;
01587   }
01588 
01589   template<>
01590   bvec Vec<std::complex<double> >::operator<=(std::complex<double>) const;
01591 
01592   template<class Num_T>
01593   bvec Vec<Num_T>::operator<=(const Num_T value) const
01594   {
01595     it_assert(datasize > 0, "Vec::operator<=: vector must have size > 0");
01596     Vec<Num_T> invector(*this);
01597     bvec temp(invector.length());
01598 
01599     for (int i=0;i<invector.length();i++)
01600       temp(i)=(invector(i)<=value);
01601 
01602     return temp;
01603   }
01604 
01605   template<>
01606   bvec Vec<std::complex<double> >::operator>(std::complex<double>) const;
01607 
01608   template<class Num_T>
01609   bvec Vec<Num_T>::operator>(const Num_T value) const
01610   {
01611     it_assert(datasize > 0, "Vec::operator>: vector must have size > 0");
01612     Vec<Num_T> invector(*this);
01613     bvec temp(invector.length());
01614 
01615     for (int i=0;i<invector.length();i++)
01616       temp(i)=(invector(i)>value);
01617 
01618     return temp;
01619   }
01620 
01621   template<>
01622   bvec Vec<std::complex<double> >::operator>=(std::complex<double>) const;
01623 
01624   template<class Num_T>
01625   bvec Vec<Num_T>::operator>=(const Num_T value) const
01626   {
01627     it_assert(datasize > 0, "Vec::operator>=: vector must have size > 0");
01628     Vec<Num_T> invector(*this);
01629     bvec temp(invector.length());
01630 
01631     for (int i=0;i<invector.length();i++)
01632       temp(i)=(invector(i)>=value);
01633 
01634     return temp;
01635   }
01636 
01637   template<class Num_T>
01638   bool Vec<Num_T>::operator==(const Vec<Num_T> &invector) const
01639   {
01640     // OBS ! if wrong size, return false
01641     if (datasize!=invector.datasize) return false;
01642     for (int i=0;i<datasize;i++) {
01643       if (data[i]!=invector.data[i]) return false;
01644     }
01645     return true;
01646   }
01647 
01648   template<class Num_T>
01649   bool Vec<Num_T>::operator!=(const Vec<Num_T> &invector) const
01650   {
01651     if (datasize!=invector.datasize) return true;
01652     for (int i=0;i<datasize;i++) {
01653       if (data[i]!=invector.data[i]) return true;
01654     }
01655     return false;
01656   }
01657 
01659   template<class Num_T>
01660   std::ostream &operator<<(std::ostream &os, const Vec<Num_T> &v)
01661   {
01662     int i, sz=v.length();
01663 
01664     os << "[" ;
01665     for (i=0; i<sz; i++) {
01666       os << v(i) ;
01667       if (i < sz-1)
01668   os << " ";
01669     }
01670     os << "]" ;
01671 
01672     return os;
01673   }
01674 
01676   template<class Num_T>
01677   std::istream &operator>>(std::istream &is, Vec<Num_T> &v)
01678   {
01679     std::ostringstream buffer;
01680     bool started = false;
01681     bool finished = false;
01682     bool brackets = false;
01683     char c;
01684 
01685     while (!finished) {
01686       if (is.eof()) {
01687         finished = true;
01688       } else {
01689         c = is.get();
01690 
01691         if (is.eof() || (c == '\n')) {
01692           if (brackets) {
01693             // Right bracket missing
01694             is.setstate(std::ios_base::failbit);
01695             finished = true;
01696           } else if (!((c == '\n') && !started)) {
01697             finished = true;
01698           }
01699         } else if ((c == ' ') || (c == '\t')) {
01700           if (started) {
01701             buffer << ' ';
01702           }
01703         } else if (c == '[') {
01704           if (started) {
01705             // Unexpected left bracket
01706             is.setstate(std::ios_base::failbit);
01707             finished = true;
01708           } else {
01709             started = true;
01710             brackets = true;
01711           }
01712         } else if (c == ']') {
01713           if (!started || !brackets) {
01714             // Unexpected right bracket
01715             is.setstate(std::ios_base::failbit);
01716             finished = true;
01717           } else {
01718             finished = true;
01719           }
01720           while (!is.eof() && (((c = is.peek()) == ' ') || (c == '\t'))) {
01721             is.get();
01722           }
01723           if (!is.eof() && (c == '\n')) {
01724             is.get();
01725           }
01726         } else {
01727           started = true;
01728           buffer << c;
01729         }
01730       }
01731     }
01732 
01733     if (!started) {
01734       v.set_size(0, false);
01735     } else {
01736       v.set(buffer.str());
01737     }
01738 
01739     return is;
01740   }
01741 
01743 
01744   // ----------------------------------------------------------------------
01745   // Instantiations
01746   // ----------------------------------------------------------------------
01747 
01748 #ifdef HAVE_EXTERN_TEMPLATE
01749 
01750   extern template class Vec<double>;
01751   extern template class Vec<int>;
01752   extern template class Vec<short int>;
01753   extern template class Vec<std::complex<double> >;
01754   extern template class Vec<bin>;
01755 
01756   // addition operator
01757 
01758   extern template const vec operator+(const vec &v1, const vec &v2);
01759   extern template const cvec operator+(const cvec &v1, const cvec &v2);
01760   extern template const ivec operator+(const ivec &v1, const ivec &v2);
01761   extern template const svec operator+(const svec &v1, const svec &v2);
01762   extern template const bvec operator+(const bvec &v1, const bvec &v2);
01763 
01764   extern template const vec operator+(const vec &v1, double t);
01765   extern template const cvec operator+(const cvec &v1, std::complex<double> t);
01766   extern template const ivec operator+(const ivec &v1, int t);
01767   extern template const svec operator+(const svec &v1, short t);
01768   extern template const bvec operator+(const bvec &v1, bin t);
01769 
01770   extern template const vec operator+(double t, const vec &v1);
01771   extern template const cvec operator+(std::complex<double> t, const cvec &v1);
01772   extern template const ivec operator+(int t, const ivec &v1);
01773   extern template const svec operator+(short t, const svec &v1);
01774   extern template const bvec operator+(bin t, const bvec &v1);
01775 
01776   // subraction operator
01777 
01778   extern template const vec operator-(const vec &v1, const vec &v2);
01779   extern template const cvec operator-(const cvec &v1, const cvec &v2);
01780   extern template const ivec operator-(const ivec &v1, const ivec &v2);
01781   extern template const svec operator-(const svec &v1, const svec &v2);
01782   extern template const bvec operator-(const bvec &v1, const bvec &v2);
01783 
01784   extern template const vec operator-(const vec &v, double t);
01785   extern template const cvec operator-(const cvec &v, std::complex<double> t);
01786   extern template const ivec operator-(const ivec &v, int t);
01787   extern template const svec operator-(const svec &v, short t);
01788   extern template const bvec operator-(const bvec &v, bin t);
01789 
01790   extern template const vec operator-(double t, const vec &v);
01791   extern template const cvec operator-(std::complex<double> t, const cvec &v);
01792   extern template const ivec operator-(int t, const ivec &v);
01793   extern template const svec operator-(short t, const svec &v);
01794   extern template const bvec operator-(bin t, const bvec &v);
01795 
01796   // unary minus
01797 
01798   extern template const vec operator-(const vec &v);
01799   extern template const cvec operator-(const cvec &v);
01800   extern template const ivec operator-(const ivec &v);
01801   extern template const svec operator-(const svec &v);
01802   extern template const bvec operator-(const bvec &v);
01803 
01804   // multiplication operator
01805 
01806 #if !defined(HAVE_BLAS)
01807   extern template double dot(const vec &v1, const vec &v2);
01808   extern template std::complex<double> dot(const cvec &v1, const cvec &v2);
01809 #endif
01810   extern template int dot(const ivec &v1, const ivec &v2);
01811   extern template short dot(const svec &v1, const svec &v2);
01812   extern template bin dot(const bvec &v1, const bvec &v2);
01813 
01814 #if !defined(HAVE_BLAS)
01815   extern template double operator*(const vec &v1, const vec &v2);
01816   extern template std::complex<double> operator*(const cvec &v1,
01817                                                  const cvec &v2);
01818 #endif
01819   extern template int operator*(const ivec &v1, const ivec &v2);
01820   extern template short operator*(const svec &v1, const svec &v2);
01821   extern template bin operator*(const bvec &v1, const bvec &v2);
01822 
01823 #if !defined(HAVE_BLAS)
01824   extern template const mat outer_product(const vec &v1, const vec &v2,
01825             bool hermitian);
01826 #endif
01827   extern template const imat outer_product(const ivec &v1, const ivec &v2,
01828              bool hermitian);
01829   extern template const smat outer_product(const svec &v1, const svec &v2,
01830              bool hermitian);
01831   extern template const bmat outer_product(const bvec &v1, const bvec &v2,
01832              bool hermitian);
01833 
01834   extern template const vec operator*(const vec &v, double t);
01835   extern template const cvec operator*(const cvec &v, std::complex<double> t);
01836   extern template const ivec operator*(const ivec &v, int t);
01837   extern template const svec operator*(const svec &v, short t);
01838   extern template const bvec operator*(const bvec &v, bin t);
01839 
01840   extern template const vec operator*(double t, const vec &v);
01841   extern template const cvec operator*(std::complex<double> t, const cvec &v);
01842   extern template const ivec operator*(int t, const ivec &v);
01843   extern template const svec operator*(short t, const svec &v);
01844   extern template const bvec operator*(bin t, const bvec &v);
01845 
01846   // elementwise multiplication
01847 
01848   extern template const vec elem_mult(const vec &a, const vec &b);
01849   extern template const cvec elem_mult(const cvec &a, const cvec &b);
01850   extern template const ivec elem_mult(const ivec &a, const ivec &b);
01851   extern template const svec elem_mult(const svec &a, const svec &b);
01852   extern template const bvec elem_mult(const bvec &a, const bvec &b);
01853 
01854   extern template void elem_mult_out(const vec &a, const vec &b, vec &out);
01855   extern template void elem_mult_out(const cvec &a, const cvec &b, cvec &out);
01856   extern template void elem_mult_out(const ivec &a, const ivec &b, ivec &out);
01857   extern template void elem_mult_out(const svec &a, const svec &b, svec &out);
01858   extern template void elem_mult_out(const bvec &a, const bvec &b, bvec &out);
01859 
01860   extern template const vec elem_mult(const vec &a, const vec &b,
01861                                       const vec &c);
01862   extern template const cvec elem_mult(const cvec &a, const cvec &b,
01863                                        const cvec &c);
01864   extern template const ivec elem_mult(const ivec &a, const ivec &b,
01865                                        const ivec &c);
01866   extern template const svec elem_mult(const svec &a, const svec &b,
01867                                        const svec &c);
01868   extern template const bvec elem_mult(const bvec &a, const bvec &b,
01869                                        const bvec &c);
01870 
01871   extern template void elem_mult_out(const vec &a, const vec &b, const vec &c,
01872                                      vec &out);
01873   extern template void elem_mult_out(const cvec &a, const cvec &b,
01874                                      const cvec &c, cvec &out);
01875   extern template void elem_mult_out(const ivec &a, const ivec &b,
01876                                      const ivec &c, ivec &out);
01877   extern template void elem_mult_out(const svec &a, const svec &b,
01878                                      const svec &c, svec &out);
01879   extern template void elem_mult_out(const bvec &a, const bvec &b,
01880                                      const bvec &c, bvec &out);
01881 
01882   extern template const vec elem_mult(const vec &a, const vec &b, const vec &c,
01883                                       const vec &d);
01884   extern template const cvec elem_mult(const cvec &a, const cvec &b,
01885                                        const cvec &c, const cvec &d);
01886   extern template const ivec elem_mult(const ivec &a, const ivec &b,
01887                                        const ivec &c, const ivec &d);
01888   extern template const svec elem_mult(const svec &a, const svec &b,
01889                                        const svec &c, const svec &d);
01890   extern template const bvec elem_mult(const bvec &a, const bvec &b,
01891                                        const bvec &c, const bvec &d);
01892 
01893   extern template void elem_mult_out(const vec &a, const vec &b, const vec &c,
01894                                      const vec &d, vec &out);
01895   extern template void elem_mult_out(const cvec &a, const cvec &b,
01896                                      const cvec &c, const cvec &d, cvec &out);
01897   extern template void elem_mult_out(const ivec &a, const ivec &b,
01898                                      const ivec &c, const ivec &d, ivec &out);
01899   extern template void elem_mult_out(const svec &a, const svec &b,
01900                                      const svec &c, const svec &d, svec &out);
01901   extern template void elem_mult_out(const bvec &a, const bvec &b,
01902                                      const bvec &c, const bvec &d, bvec &out);
01903 
01904   // in-place elementwise multiplication
01905 
01906   extern template void elem_mult_inplace(const vec &a, vec &b);
01907   extern template void elem_mult_inplace(const cvec &a, cvec &b);
01908   extern template void elem_mult_inplace(const ivec &a, ivec &b);
01909   extern template void elem_mult_inplace(const svec &a, svec &b);
01910   extern template void elem_mult_inplace(const bvec &a, bvec &b);
01911 
01912   // elementwise multiplication followed by summation
01913 
01914   extern template double elem_mult_sum(const vec &a, const vec &b);
01915   extern template std::complex<double> elem_mult_sum(const cvec &a,
01916                                                      const cvec &b);
01917   extern template int elem_mult_sum(const ivec &a, const ivec &b);
01918   extern template short elem_mult_sum(const svec &a, const svec &b);
01919   extern template bin elem_mult_sum(const bvec &a, const bvec &b);
01920 
01921   // division operator
01922 
01923   extern template const vec operator/(const vec &v, double t);
01924   extern template const cvec operator/(const cvec &v, std::complex<double> t);
01925   extern template const ivec operator/(const ivec &v, int t);
01926   extern template const svec operator/(const svec &v, short t);
01927   extern template const bvec operator/(const bvec &v, bin t);
01928 
01929   extern template const vec operator/(double t, const vec &v);
01930   extern template const cvec operator/(std::complex<double> t, const cvec &v);
01931   extern template const ivec operator/(int t, const ivec &v);
01932   extern template const svec operator/(short t, const svec &v);
01933   extern template const bvec operator/(bin t, const bvec &v);
01934 
01935   // elementwise division operator
01936 
01937   extern template const vec elem_div(const vec &a, const vec &b);
01938   extern template const cvec elem_div(const cvec &a, const cvec &b);
01939   extern template const ivec elem_div(const ivec &a, const ivec &b);
01940   extern template const svec elem_div(const svec &a, const svec &b);
01941   extern template const bvec elem_div(const bvec &a, const bvec &b);
01942 
01943   extern template const vec elem_div(double t, const vec &v);
01944   extern template const cvec elem_div(std::complex<double> t, const cvec &v);
01945   extern template const ivec elem_div(int t, const ivec &v);
01946   extern template const svec elem_div(short t, const svec &v);
01947   extern template const bvec elem_div(bin t, const bvec &v);
01948 
01949   extern template void elem_div_out(const vec &a, const vec &b, vec &out);
01950   extern template void elem_div_out(const cvec &a, const cvec &b, cvec &out);
01951   extern template void elem_div_out(const ivec &a, const ivec &b, ivec &out);
01952   extern template void elem_div_out(const svec &a, const svec &b, svec &out);
01953   extern template void elem_div_out(const bvec &a, const bvec &b, bvec &out);
01954 
01955   // elementwise division followed by summation
01956 
01957   extern template double elem_div_sum(const vec &a, const vec &b);
01958   extern template std::complex<double> elem_div_sum(const cvec &a,
01959                                                     const cvec &b);
01960   extern template int elem_div_sum(const ivec &a, const ivec &b);
01961   extern template short elem_div_sum(const svec &a, const svec &b);
01962   extern template bin elem_div_sum(const bvec &a, const bvec &b);
01963 
01964   // concat operator
01965 
01966   extern template const vec concat(const vec &v, double a);
01967   extern template const cvec concat(const cvec &v, std::complex<double> a);
01968   extern template const ivec concat(const ivec &v, int a);
01969   extern template const svec concat(const svec &v, short a);
01970   extern template const bvec concat(const bvec &v, bin a);
01971 
01972   extern template const vec concat(double a, const vec &v);
01973   extern template const cvec concat(std::complex<double> a, const cvec &v);
01974   extern template const ivec concat(int a, const ivec &v);
01975   extern template const svec concat(short a, const svec &v);
01976   extern template const bvec concat(bin a, const bvec &v);
01977 
01978   extern template const vec concat(const vec &v1, const vec &v2);
01979   extern template const cvec concat(const cvec &v1, const cvec &v2);
01980   extern template const ivec concat(const ivec &v1, const ivec &v2);
01981   extern template const svec concat(const svec &v1, const svec &v2);
01982   extern template const bvec concat(const bvec &v1, const bvec &v2);
01983 
01984   extern template const vec concat(const vec &v1, const vec &v2,
01985                                    const vec &v3);
01986   extern template const cvec concat(const cvec &v1, const cvec &v2,
01987                                     const cvec &v3);
01988   extern template const ivec concat(const ivec &v1, const ivec &v2,
01989                                     const ivec &v3);
01990   extern template const svec concat(const svec &v1, const svec &v2,
01991                                     const svec &v3);
01992   extern template const bvec concat(const bvec &v1, const bvec &v2,
01993                                     const bvec &v3);
01994 
01995   extern template const vec concat(const vec &v1, const vec &v2, const vec &v3,
01996                                    const vec &v4);
01997   extern template const cvec concat(const cvec &v1, const cvec &v2,
01998                                     const cvec &v3, const cvec &v4);
01999   extern template const ivec concat(const ivec &v1, const ivec &v2,
02000                                     const ivec &v3, const ivec &v4);
02001   extern template const svec concat(const svec &v1, const svec &v2,
02002                                     const svec &v3, const svec &v4);
02003   extern template const bvec concat(const bvec &v1, const bvec &v2,
02004                                     const bvec &v3, const bvec &v4);
02005 
02006   extern template const vec concat(const vec &v1, const vec &v2, const vec &v3,
02007                                    const vec &v4, const vec &v5);
02008   extern template const cvec concat(const cvec &v1, const cvec &v2,
02009                                     const cvec &v3, const cvec &v4,
02010                                     const cvec &v5);
02011   extern template const ivec concat(const ivec &v1, const ivec &v2,
02012                                     const ivec &v3, const ivec &v4,
02013                                     const ivec &v5);
02014   extern template const svec concat(const svec &v1, const svec &v2,
02015                                     const svec &v3, const svec &v4,
02016                                     const svec &v5);
02017   extern template const bvec concat(const bvec &v1, const bvec &v2,
02018                                     const bvec &v3, const bvec &v4,
02019                                     const bvec &v5);
02020 
02021   // I/O streams
02022 
02023   extern template std::ostream &operator<<(std::ostream& os, const vec &vect);
02024   extern template std::ostream &operator<<(std::ostream& os, const cvec &vect);
02025   extern template std::ostream &operator<<(std::ostream& os, const svec &vect);
02026   extern template std::ostream &operator<<(std::ostream& os, const ivec &vect);
02027   extern template std::ostream &operator<<(std::ostream& os, const bvec &vect);
02028   extern template std::istream &operator>>(std::istream& is, vec &vect);
02029   extern template std::istream &operator>>(std::istream& is, cvec &vect);
02030   extern template std::istream &operator>>(std::istream& is, svec &vect);
02031   extern template std::istream &operator>>(std::istream& is, ivec &vect);
02032   extern template std::istream &operator>>(std::istream& is, bvec &vect);
02033 
02034 #endif // HAVE_EXTERN_TEMPLATE
02035 
02037 
02038 } // namespace itpp
02039 
02040 #endif // #ifndef VEC_H
SourceForge Logo

Generated on Sun Dec 9 17:31:01 2007 for IT++ by Doxygen 1.5.4