/*---------------------------------------------------------------------------* * IT++ * *---------------------------------------------------------------------------* * Copyright (c) 1995-2004 by Tony Ottosson, Thomas Eriksson, Pål Frenger, * * Tobias Ringström, and Jonas Samuelsson. * * * * Permission to use, copy, modify, and distribute this software and its * * documentation under the terms of the GNU General Public License is hereby * * granted. No representations are made about the suitability of this * * software for any purpose. It is provided "as is" without expressed or * * implied warranty. See the GNU General Public License for more details. * *---------------------------------------------------------------------------*/ /*! \file \brief Templated Vector Class Definitions \author Tony Ottosson and Tobias Ringström 1.44 2004/10/28 13:34:47 */ #ifndef __vec_h #define __vec_h #include #include #include #include #include #include #include #include "itconfig.h" //#include "base/mat.h" is done below Vec class definition (needed for g++ 3.4) #include "base/binary.h" #include "base/itassert.h" #include "base/scalfunc.h" #include "base/factory.h" #include "../src/base/copy_vector.h" //#ifndef NO_CBLAS //#include "base/cblas.h" //#endif using std::cout; using std::endl; using std::string; using std::ostream; using std::istream; using std::istringstream; using std::getline; using std::complex; namespace itpp { // Declaration of Vec template class Vec; // Declaration of Mat template class Mat; // Declaration of bin class bin; //----------------------------------------------------------------------------------- // Declaration of Vec Friends //----------------------------------------------------------------------------------- // Addition of two vectors template Vec operator+(const Vec &v1, const Vec &v2); // Addition of a vector and a scalar template Vec operator+(const Vec &v, Num_T t); // Addition of a scalar and a vector template Vec operator+(Num_T t, const Vec &v); // Subtraction of a vector from a vector template Vec operator-(const Vec &v1, const Vec &v2); // Subtraction of a scalar from a vector template Vec operator-(const Vec &v, Num_T t); // Subtraction of vector from scalar. Results in a vector template Vec operator-(Num_T t, const Vec &v); // Negation of vector template Vec operator-(const Vec &v); // Inner (dot) product of two vectors v1 and v2 template Num_T dot(const Vec &v1, const Vec &v2); // Inner (dot) product of two vectors v1 and v2 template Num_T operator*(const Vec &v1, const Vec &v2) { return dot(v1, v2); } // Outer product of two vectors v1 and v2 template Mat outer_product(const Vec &v1, const Vec &v2); // Multiplication of a vector and a scalar template Vec operator*(const Vec &v, Num_T t); // Multiplication of a scalar and a vector. Results in a vector template Vec operator*(Num_T t, const Vec &v); // Elementwise multiplication of the two vectors template Vec elem_mult(const Vec &v1, const Vec &v2); // Elementwise multiplication of the three vectors template Vec elem_mult(const Vec &v1, const Vec &v2, const Vec &v3); // Elementwise multiplication of the four vectors template Vec elem_mult(const Vec &v1, const Vec &v2, const Vec &v3, const Vec &v4); // Division of all elements in \c v with \c t template Vec operator/(const Vec &v, Num_T t); // Division of \c t with all elements in \c v template Vec operator/(const Num_T t, const Vec &v); // Elementwise division template Vec elem_div(const Vec &v1, const Vec &v2); // Elementwise division of scalar \c t and vector \c v template Vec elem_div(const Num_T t, const Vec &v); // Append element \c a to the end of the vector \c v template Vec concat(const Vec &v, const Num_T a); // Concat element \c a to the beginning of the vector \c v template Vec concat(const Num_T a, const Vec &v); // Concat vectors \c v1 and \c v2 template Vec concat(const Vec &v1,const Vec &v2); // Concat vectors \c v1, \c v2 and \c v3 template Vec concat(const Vec &v1, const Vec &v2, const Vec &v3); // Concat vectors \c v1, \c v2, \c v3 and \c v4 template Vec concat(const Vec &v1, const Vec &v2, const Vec &v3, const Vec &v4); // Concat vectors \c v1, \c v2 \c v3, \c v4 and \c v5 template Vec concat(const Vec &v1, const Vec &v2, const Vec &v3, const Vec &v4, const Vec &v5); //----------------------------------------------------------------------------------- // Declaration of Vec //----------------------------------------------------------------------------------- /*! \brief Templated vectors \author Tony Ottosson and Tobias Ringstrom Vectors can be of arbitrarily types, but conversions and functions are prepared for \c bin, \c short, \c int, \c double, and \c complex vectors and these are predefined as: \c bvec, \c svec, \c ivec, \c vec, and \c cvec. \c double and \c complex are \c double and \c complex respectively. Examples: Vector Constructors: When constructing a vector without a length (memory) use \code vec temp; \endcode For construction of a vector of a given length use \code vec temp(length); \endcode It is also possible to assign the constructed vector the value and size of another vector by \code vec temp(invector); \endcode If you have explicit values you would like to assign to the vector it is possible to do this using strings as: \code vec a("0 0.7 5 9.3"); // that is a = [0, 0.7, 5, 9.3] vec a="0 0.7 5 9.3"; // the constructor are called implicitly ivec b="0:5"; // that is b = [0, 1, 2, 3, 4, 5] vec c="3:2.5:13"; // that is c = [3, 5.5, 8, 10.5, 13] \endcode It is also possible to change length by \code temp.set_size(new_length, false); \endcode where \c false is used to indicate that the old values in \c temp is not copied. If you like to preserve the values use \c true. There are a number of methods to access parts of a vector. Examples are \code a(5); // Element number 5 a(5,9); // Elements 5, 6, 7, 8, and 9 a.left(10); // The 10 most left elements (the first) a.right(10); // The 10 most right elements (the last) a.mid(5, 7); // 7 elements starting from element 5 \endcode It is also possible to modify parts of a vector as e.g. in \code a.del(5); // deletes element number 5 a.ins(3.4, 9); // inserts the element 3.4 at position 9 a.replace_mid(12, b); // replaces elements from 12 with the vector b \endcode It is of course also possible to perform the common linear algebra methods such as addition, subtraction, and scalar product (*). Observe though, that vectors are assumed to be column-vectors in operations with matrices. Most elementary functions such as sin(), cosh(), log(), abs(), ..., are available as operations on the individual elements of the vectors. Please see the individual functions for more details. By default, the Vec elements are created using the default constructor for the element type. This can be changed by specifying a suitable Factory in the Vec constructor call; see Detailed Description for Factory. */ template class Vec { public: //! Default constructor. An element factory \c f can be specified explicit Vec(const Factory &f = DEFAULT_FACTORY) : factory(f) { init(); } //! Constructor. An element factory \c f can be specified explicit Vec(int size, const Factory &f = DEFAULT_FACTORY) : factory(f) { it_assert1(size>=0,"Negative size in Vec::Vec(int)"); init(); alloc(size); } //! Copy constructor Vec(const Vec &v); //! Constructor, similar to the copy constructor, but also takes an element factory \c f as argument Vec(const Vec &v, const Factory &f); //! Constructor. An element factory \c f can be specified Vec(const char *values, const Factory &f = DEFAULT_FACTORY) : factory(f) { init(); set(values); } //! Constructor. An element factory \c f can be specified Vec(const std::string &values, const Factory &f = DEFAULT_FACTORY) : factory(f) { init(); set(values); } //! Constructor taking a C-array as input. Copies all data. An element factory \c f can be specified Vec(Num_T *c_array, int size, const Factory &f = DEFAULT_FACTORY) : factory(f) { init(); alloc(size); copy_vector(size, c_array, data); } //! Destructor ~Vec() { free(); } //! The size of the vector int length() const { return datasize; } //! The size of the vector int size() const { return datasize; } //! Set length of vector. if copy = true then keeping the old values void set_length(int size, bool copy=false) { set_size(size,copy); } //! Set length of vector. if copy = true then keeping the old values void set_size(int size, bool copy=false); //! Set the vector to the all zero vector void zeros() { for (int i=0; i=0&&i=0&&i=0&&i=0&&i operator()(int i1, int i2) const; //! Sub-vector where the elements are given by the list \c indexlist const Vec operator()(const Vec &indexlist) const; //! Accessor-style method. First element is 0 const Num_T &get(int i) const { it_assert0(i>=0&&i get(int i1, int i2) const; //! Modifier-style method. First element is 0 void set(int i, const Num_T &v) { it_assert0(i>=0&&i transpose() const; //! Matrix transpose. Converts to a matrix with the vector in the first row Mat T() const { return this->transpose(); } //! Hermitian matrix transpose. Converts to a matrix with the conjugate of the vector in the first row Mat hermitian_transpose() const; //! Hermitian matrix transpose. Converts to a matrix with the conjugate of the vector in the first row Mat H() const { return this->hermitian_transpose(); } //! Addition of vector void operator+=(const Vec &v); //! Addition of scalar void operator+=(Num_T t) { for (int i=0;i operator+<>(const Vec &v1, const Vec &v2); //! Addition of a vector and a scalar friend Vec operator+<>(const Vec &v, Num_T t); //! Addition of a scalar and a vector friend Vec operator+<>(Num_T t, const Vec &v); //! Subtraction of vector void operator-=(const Vec &v); //! Subtraction of scalar void operator-=(Num_T t) { for (int i=0;i operator-<>(const Vec &v1, const Vec &v2); //! Subtraction of scalar from vector friend Vec operator-<>(const Vec &v, Num_T t); //! Sutraction of vector from scalar friend Vec operator-<>(Num_T t, const Vec &v); //! Negation of vector friend Vec operator-<>(const Vec &v); //! Multiply with a scalar void operator*=(Num_T t) { for (int i=0;i(const Vec &v1, const Vec &v2); //! Inner (dot) product friend Num_T dot <>(const Vec &v1, const Vec &v2); //! Outer product of two vectors v1 and v2 friend Mat outer_product <>(const Vec &v1, const Vec &v2); //! Elementwise multiplication of vector and scalar friend Vec operator*<>(const Vec &v, Num_T t); //! Elementwise multiplication of vector and scalar friend Vec operator*<>(Num_T t, const Vec &v); //! Elementwise multiplication friend Vec elem_mult <>(const Vec &v1, const Vec &v2); //! Elementwise multiplication of three vectors friend Vec elem_mult <>(const Vec &v1, const Vec &v2, const Vec &v3); //! Elementwise multiplication of four vectors friend Vec elem_mult <>(const Vec &v1, const Vec &v2, const Vec &v3, const Vec &v4); //! Elementwise division void operator/=(Num_T t) { for (int i=0;i operator/<>(const Vec &v, Num_T t); //! Elementwise division friend Vec operator/<>(const Num_T t, const Vec &v); //! Elementwise division void operator/=(const Vec &v); //! Elementwise division friend Vec elem_div <>(const Vec &v1, const Vec &v2); //! Elementwise division friend Vec elem_div <>(const Num_T t, const Vec &v); //! Get the elements in the vector where \c binlist is \c 1 Vec get(const Vec &binlist) const; //! Get the right \c nr elements from the vector Vec right(int nr) const; //! Get the left \c nr elements from the vector Vec left(int nr) const; //! Get the middle part of vector from \c start including \c nr elements Vec mid(int start, int nr) const; //! Split the vector into two parts at element \c pos. Return the first part and keep the second. Vec split(int pos); //! Shift in element \c In at position 0 \c n times void shift_right(Num_T In, int n=1); //! Shift in vector \c In at position 0 void shift_right(const Vec &In); //! Shift out the \c n left elements and a the same time shift in the element \c at last position \c n times void shift_left(Num_T In, int n=1); //! Shift in vector \c In at last position void shift_left(const Vec &In); //! Append element \c a to the end of the vector \c v friend Vec concat <>(const Vec &v, const Num_T a); //! Concat element \c a to the beginning of the vector \c v friend Vec concat <>(const Num_T a, const Vec &v); //! Concat vectors \c v1 and \c v2 friend Vec concat <>(const Vec &v1,const Vec &v2); //! Concat vectors \c v1, \c v2 and \c v3 friend Vec concat <>(const Vec &v1, const Vec &v2, const Vec &v3); //! Set subvector defined by indicies \c i1 to \c i2 to vector \c v void set_subvector(int i1, int i2, const Vec &v); //! Set subvector defined by first index \c i and size of vector \c v to \c v void set_subvector(int i, const Vec &v); //! Set subvector defined by indicies i1 to i2 to constant t void set_subvector(int i1, int i2, const Num_T t); //! Replace the elements from \c pos by the vector \c v void replace_mid(int pos, const Vec &v); //! Delete element number \c index void del(int index); //! Delete elements from \c i1 to \c i2 void del(int i1, int i2); //! Insert element \c in at \c index void ins(int index, Num_T in); //! Insert vector \c in at \c index void ins(int index, const Vec &in); //! Assign all elements in vector to \c t void operator=(Num_T t) { for (int i=0;i &v); //! Assign vector equal to the 1-dimensional matrix \c m void operator=(const Mat &m); //! Assign vector the values in the string \c values void operator=(const char *values) { set(values); } //! Elementwise equal to the scalar Vec operator==(const Num_T value) const; //! Elementwise not-equal to the scalar Vec operator!=(const Num_T value) const; //! Elementwise less than the scalar Vec operator<(const Num_T value) const; //! Elementwise less than and equal to the scalar Vec operator<=(const Num_T value) const; //! Elementwise greater than the scalar Vec operator>(const Num_T value) const; //! Elementwise greater than and equal to the scalar Vec operator>=(const Num_T value) const; //! Compare two vectors. False if wrong sizes or different values bool operator==(const Vec &v) const; //! Compare two vectors. True if different bool operator!=(const Vec &v) const; //! Index operator without boundary check. Not recommended to use. Num_T &_elem(int i) { return data[i]; } //! Index operator without boundary check. Not recommended to use. const Num_T &_elem(int i) const { return data[i]; } //! Get the pointer to the internal structure. Not recommended to use. Num_T *_data() { return data; } //! Get the pointer to the internal structure. Not recommended to use. const Num_T *_data() const { return data; } protected: //! Allocate storage for a vector of length \c size. void alloc(int size) { if ( datasize == size ) return; free(); // Free memory (if any allocated) if (size == 0) return; create_elements(data, size, factory); datasize=size; it_assert1(data, "Out of memory in Vec::alloc()"); } //! Free the storage space allocated by the vector void free() { delete [] data; init(); } //! The current number of elements in the vector int datasize; //! A pointer to the data area Num_T *data; //! Element factory (set to DEFAULT_FACTORY to use Num_T default constructors only) const Factory &factory; private: void init() { data = 0; datasize = 0; } }; } //namespace itpp #include "base/mat.h" namespace itpp { //----------------------------------------------------------------------------------- // Declaration of input and output streams for Vec //----------------------------------------------------------------------------------- /*! \relates Vec \brief Stream output of vector */ template std::ostream &operator<<(std::ostream &os, const Vec &v); /*! \relates Vec \brief Stream input of vector The input can be on the form "1 2 3" or "[1 2 3]", i.e. with or without brackets. The first form is compatible with the set method, while the second form is compatible with the ostream operator. The elements can be separated by blank space or commas. "[]" means an empty vector. "1:4" means "1 2 3 4". "1:3:10" means every third integer from 1 to 10, i.e. "1 4 7 10". */ template std::istream &operator>>(std::istream &is, Vec &v); //----------------------------------------------------------------------------------- // Type definitions of vec, cvec, ivec, svec, and bvec //----------------------------------------------------------------------------------- /*! \relates Vec \brief Definition of double vector type */ typedef Vec vec; /*! \relates Vec \brief Definition of complex vector type */ typedef Vec > cvec; /*! \relates Vec \brief Definition of integer vector type */ typedef Vec ivec; /*! \relates Vec \brief Definition of short vector type */ typedef Vec svec; /*! \relates Vec \brief Definition of binary vector type */ typedef Vec bvec; //----------------------------------------------------------------------------------- // Implementation of templated Vec members and friends //----------------------------------------------------------------------------------- template inline Vec::Vec(const Vec &v) : factory(v.factory) { init(); alloc(v.datasize); copy_vector(datasize, v.data, data); } template inline Vec::Vec(const Vec &v, const Factory &f) : factory(f) { init(); alloc(v.datasize); copy_vector(datasize, v.data, data); } template void Vec::set_size(int size, bool copy) { it_assert1(size >= 0, "New size must not be negative in Vec::set_size()"); if (size!=datasize) { if (copy) { Vec temp(*this); alloc(size); for (int i=0; i bool cvec::set(const char *values); template<> bool bvec::set(const char *values); template bool Vec::set(const char *values) { istringstream buffer(values); Num_T b, c; int pos=0, maxpos=10; alloc(maxpos); while (buffer.peek()!=EOF) { switch (buffer.peek()) { case ':': // reads format a:b:c or a:b buffer.get(); if (!buffer.eof()) { buffer >> b; } if (!buffer.eof() && buffer.peek() == ':') { buffer.get(); if (!buffer.eof()) { buffer >> c; while (sign(b)*(data[pos-1]+b-c)<=0) { pos++; if (pos > maxpos) { maxpos=maxpos*2; set_size(maxpos, true); } data[pos-1]=data[pos-2]+b; } } } else { while (data[pos-1] maxpos) { maxpos=maxpos*2; set_size(maxpos, true); } data[pos-1]=data[pos-2]+1; } } break; case ',': buffer.get(); break; default: pos++; if (pos > maxpos) { maxpos *= 2; set_size(maxpos, true); } buffer >> data[pos-1]; while (buffer.peek()==' ') { buffer.get(); } break; } } set_size(pos, true); return true; } template bool Vec::set(const std::string &str) { return set(str.c_str()); } template inline const Vec Vec::operator()(int i1, int i2) const { if (i1 == -1) i1 = datasize-1; if (i2 == -1) i2 = datasize-1; it_assert1(i1>=0 && i2>=0 && i1::operator()(i1,i2): indicies out of range"); it_assert1(i2>=i1, "Vec::op(i1,i2): i2 >= i1 necessary"); Vec s(i2-i1+1); copy_vector(s.datasize, data+i1, s.data); return s; } template inline const Vec Vec::get(int i1, int i2) const { return (*this)(i1, i2); } template const Vec Vec::operator()(const Vec &indexlist) const { Vec temp(indexlist.length()); for (int i=0;i=0) && (indexlist(i) < datasize), "Vec::operator()(ivec &): index outside range"); temp(i)=data[indexlist(i)]; } return temp; } template Mat Vec::transpose() const { Mat temp(1, datasize); for (int i=0; i Mat > cvec::hermitian_transpose() const; template Mat Vec::hermitian_transpose() const { Mat temp(1, datasize); for (int i=0; i inline void Vec::operator+=(const Vec &v) { int i; if (datasize == 0) { // if not assigned a size. alloc(v.datasize); for (i=0; i::operator+=: wrong sizes"); for (i=0; i inline Vec operator+(const Vec &v1, const Vec &v2) { int i; Vec r(v1.datasize); it_assert1(v1.datasize==v2.datasize, "Vec::operator+: wrong sizes"); for (i=0; i inline Vec operator+(const Vec &v, Num_T t) { int i; Vec r(v.datasize); for (i=0; i inline Vec operator+(Num_T t, const Vec &v) { int i; Vec r(v.datasize); for (i=0; i inline void Vec::operator-=(const Vec &v) { int i; if (datasize == 0) { // if not assigned a size. alloc(v.datasize); for (i=0; i::operator-=: wrong sizes"); for (i=0; i inline Vec operator-(const Vec &v1, const Vec &v2) { int i; Vec r(v1.datasize); it_assert1(v1.datasize==v2.datasize, "Vec::operator-: wrong sizes"); for (i=0; i inline Vec operator-(const Vec &v, Num_T t) { int i; Vec r(v.datasize); for (i=0; i inline Vec operator-(Num_T t, const Vec &v) { int i; Vec r(v.datasize); for (i=0; i inline Vec operator-(const Vec &v) { int i; Vec r(v.datasize); for (i=0; i double dot(const Vec &v1, const Vec &v2); template<> complex dot(const Vec > &v1, const Vec > &v2); #endif //NO_CBLAS template inline Num_T dot(const Vec &v1, const Vec &v2) { int i; Num_T r=Num_T(0); it_assert1(v1.datasize==v2.datasize, "Vec::dot: wrong sizes"); for (i=0; i inline Mat outer_product(const Vec &v1, const Vec &v2) { int i, j; it_assert1(v1.datasize>0 && v2.datasize>0, "outer_product:: Vector of zero size"); Mat r(v1.datasize, v2.datasize); for (i=0; i inline Vec operator*(const Vec &v, Num_T t) { int i; Vec r(v.datasize); for (i=0; i inline Vec operator*(Num_T t, const Vec &v) { int i; Vec r(v.datasize); for (i=0; i inline Vec elem_mult(const Vec &v1, const Vec &v2) { int i; Vec r(v1.datasize); it_assert1(v1.datasize==v2.datasize, "Vec::elem_mult: wrong sizes"); for (i=0; i inline Vec elem_mult(const Vec &v1, const Vec &v2, const Vec &v3) { int i; Vec r(v1.datasize); it_assert1(v1.datasize==v2.datasize, "Vec::elem_mult: wrong sizes"); it_assert1(v2.datasize==v3.datasize, "Vec::elem_mult: wrong sizes"); for (i=0; i inline Vec elem_mult(const Vec &v1, const Vec &v2, const Vec &v3, const Vec &v4) { int i; Vec r(v1.datasize); it_assert1(v1.datasize==v2.datasize, "Vec::elem_mult: wrong sizes"); it_assert1(v2.datasize==v3.datasize, "Vec::elem_mult: wrong sizes"); it_assert1(v3.datasize==v4.datasize, "Vec::elem_mult: wrong sizes"); for (i=0; i inline Vec operator/(const Vec &v, Num_T t) { int i; Vec r(v.datasize); for (i=0; i inline Vec operator/(const Num_T t, const Vec &v) { int i; Vec r(v.datasize); for (i=0; i inline void Vec::operator/=(const Vec &v) { int i; it_assert1(datasize==v.datasize, "Vec::operator/=: wrong sizes"); for (i=0; i inline Vec elem_div(const Vec &v1, const Vec &v2) { int i; Vec r(v1.datasize); it_assert1(v1.datasize==v2.datasize, "elem_div: wrong sizes"); for (i=0; i inline Vec elem_div(const Num_T t, const Vec &v) { int i; Vec r(v.datasize); for (i=0; i Vec Vec::get(const Vec &binlist) const { it_assert1(datasize == binlist.size(), "Vec::get(bvec &): wrong sizes"); Vec temp(binlist.length()); int j=0; for (int i=0;i inline Vec Vec::right(int nr) const { it_assert1(nr<=datasize, "Vec::right: index out of range"); Vec temp(nr); if (nr!=0) { copy_vector(nr, &data[datasize-nr], &temp[0]); } return temp; } template inline Vec Vec::left(int nr) const { it_assert1(nr<=datasize, "Vec::left: index out of range"); Vec temp(nr); if (nr!=0) { copy_vector(nr, &data[0], &temp[0]); } return temp; } template inline Vec Vec::mid(int start, int nr) const { it_assert1((start>=0)&& ((start+nr)<=datasize), "Vec::mid: indexing out of range"); Vec temp(nr); if (nr!=0) { copy_vector(nr, &data[start], &temp[0]); } return temp; } template Vec Vec::split(int Position) { it_assert1((Position>=0) && (Position<=datasize), "Vec::split: index out of range"); Vec Temp1(Position); Vec Temp2(datasize-Position); int i; for (i=0;i void Vec::shift_right(Num_T In, int n) { int i=datasize; it_assert1(n>=0, "Vec::shift_right: index out of range"); while (--i >= n) data[i] = data[i-n]; while (i >= 0) data[i--] = In; } template void Vec::shift_right(const Vec &In) { int i; for (i=datasize-1; i>=In.datasize; i--) data[i]=data[i-In.datasize]; for (i=0; i