Changeset 995 for branches/peter-dev

branches/peter-dev/test/vector_test.cc

r988	r995
34	34
35	35	using namespace theplu::yat;
	36
36	37
37	38	void check_file_access(std::string& str)

branches/peter-dev/yat/classifier/NBC.cc

-                      r963
+                      r995
     // normalize each row (label) to sum up to unity (probability)
     for (size_t i=0; i<prediction.rows(); ++i)
       utility::vector(prediction,i) *=
 .0/utility::sum(utility::vector(prediction,i));
+    for (size_t i=0; i<prediction.rows(); ++i){
+      prediction.row_vec(i) *=
+.0/sum(prediction.row_vec(i));
+    }
+  }

branches/peter-dev/yat/classifier/SVM.cc

r964	r995
198	198	{
199	199	trained_=false;
200		alpha_~~.clone(utility::vector(target_.size(), 0)~~);
	200	alpha_ = utility::vector(target_.size(), 0);
201	201	}
202	202

branches/peter-dev/yat/classifier/utility.cc

-                      r865
+                      r995
   void convert(const DataLookup1D& lookup, utility::vector& vector)
+  {
     vector.clone(utility::vector(lookup.size()));
+    vector = utility::vector(lookup.size());
     for(u_int i=0; i<lookup.size(); i++)
       vector(i)=lookup(i);
 …
+  {
     value.clone(utility::vector(lookup.size()));
     weight.clone(utility::vector(lookup.size()));
+    value = utility::vector(lookup.size());
+    weight = utility::vector(lookup.size());
     for(u_int i=0; i<lookup.size(); i++){
       value(i)=lookup.data(i);

branches/peter-dev/yat/regression/Local.cc

-                      r865
+                      r995
 #include "OneDimensionalWeighted.h"
 #include "yat/utility/vector.h"
+#include "yat/utility/vectorView.h"
 #include <algorithm>
 …
     size_t nof_fits=data_.size()/step_size;
     x_.clone(utility::vector(nof_fits));
     y_predicted_.clone(utility::vector(x_.size()));
     y_err_.clone(utility::vector(x_.size()));
+    x_ = utility::vector(nof_fits);
+    y_predicted_ = utility::vector(x_.size());
+    y_err_ = utility::vector(x_.size());
     sort(data_.begin(), data_.end());
 …
       assert(max_index<data_.size());
       utility::vector x_local(x, min_index, max_index-min_index+1);
       utility::vector y_local(y, min_index, max_index-min_index+1);
+      utility::vectorView x_local(x, min_index, max_index-min_index+1);
+      utility::vectorView y_local(y, min_index, max_index-min_index+1);
       // calculating weights

branches/peter-dev/yat/regression/MultiDimensional.cc

r865	r995
58	58	assert(x.rows()==y.size());
59	59	covariance_.clone(utility::matrix(x.columns(),x.columns()));
60		fit_parameters_~~.clone(utility::vector(x.columns()~~));
	60	fit_parameters_ = utility::vector(x.columns());
61	61	if (work_)
62	62	gsl_multifit_linear_free(work_);

branches/peter-dev/yat/regression/MultiDimensionalWeighted.cc

r916	r995
59	59
60	60	covariance_.clone(utility::matrix(x.columns(),x.columns()));
61		fit_parameters_~~.clone(utility::vector(x.columns()~~));
	61	fit_parameters_ = utility::vector(x.columns());
62	62	if (work_)
63	63	gsl_multifit_linear_free(work_);

branches/peter-dev/yat/statistics/AveragerPairWeighted.h

r937	r995
201	201	{
202	202	for (size_t i=0; i<x.size(); ++i){
203		utility::yat_assert<std::runtime_error>(!std::isnan(x~~[i]~~));
	203	utility::yat_assert<std::runtime_error>(!std::isnan(x(i)));
204	204	add(x[i],y[i],wx[i],wy[i]);
205	205	}

branches/peter-dev/yat/utility/Makefile.am

-                      r981
+                      r995
   Alignment.cc ColumnStream.cc CommandLine.cc FileUtil.cc kNNI.cc \
   matrix.cc NNI.cc Option.cc OptionFile.cc OptionHelp.cc OptionSwitch.cc \
+  PCA.cc stl_utility.cc SVD.cc TypeInfo.cc utility.cc vector.cc WeNNI.cc
+  PCA.cc stl_utility.cc SVD.cc TypeInfo.cc utility.cc vector.cc \
+  vectorBase.cc vectorConstView.cc vectorMutable.cc vectorView.cc WeNNI.cc Peter.cc
 include_utilitydir = $(includedir)/yat/utility
 …
   IteratorWeighted.h kNNI.h matrix.h NNI.h \
   Option.h OptionArg.h OptionFile.h OptionHelp.h OptionSwitch.h \
+  PCA.h stl_utility.h SVD.h TypeInfo.h utility.h vector.h WeNNI.h \
+  yat_assert.h
+  PCA.h stl_utility.h SVD.h TypeInfo.h utility.h vector.h
+  vectorBase.h vectorMutable.h vectorView.h vectorConstView.h
+  WeNNI.h yat_assert.h

branches/peter-dev/yat/utility/PCA.cc

-                      r865
+                      r995
 #include "SVD.h"
 #include "utility.h"
+#include "vectorConstView.h"
 #include <iostream>
 …
     eigenvectors_.clone(U);
     eigenvectors_ .transpose();
+    eigenvalues_ = pSVD->s();
+    // T
+    for( size_t i = 0; i < eigenvalues_.size(); ++i )
+      eigenvalues_(i) = eigenvalues_(i)*eigenvalues_(i);
+    eigenvalues_ *= 1.0/A_center.rows();
+    // Sort the eigenvectors in order of eigenvalues
+    // Simple (not efficient) algorithm that always
+    // make sure that the i:th element is in its correct
+    // position (N element --> Ordo( N*N ))
+    for ( size_t i = 0; i < eigenvalues_.size(); ++i )
+      for ( size_t j = i + 1; j < eigenvalues_.size(); ++j )
+        if ( eigenvalues_(j) > eigenvalues_(i) ) {
+          std::swap( eigenvalues_(i), eigenvalues_(j) );
+          eigenvectors_.swap_rows(i,j);
+        }
+  }
+  /*
+  void PCA::process_transposed_problem(void)
+  {
+    // Row-center the data matrix
+    utility::matrix A_center( A_.rows(), A_.columns() );
+    this->row_center( A_center );
+    // Transform into SVD friendly dimension
+    A_.transpose();
+    A_center.transpose();
+    // Single value decompose the data matrix
+    std::auto_ptr<SVD> pSVD( new SVD( A_center ) );
+    pSVD->decompose();
+    utility::matrix U(pSVD->U());
+    utility::matrix V(pSVD->V());
+    // Read the eigenvectors and eigenvalues
+    eigenvectors_.clone(V);
+    eigenvectors_.transpose();
     eigenvalues_.clone(pSVD->s());
+    // Transform back when done with SVD!
+    // (used V insted of U now for eigenvectors)
+    A_.transpose();
+    A_center.transpose();
     // T
 …
+        }
+  }
-  /*
-  void PCA::process_transposed_problem(void)
+  {
-    // Row-center the data matrix
-    utility::matrix A_center( A_.rows(), A_.columns() );
-    this->row_center( A_center );
-    // Transform into SVD friendly dimension
-    A_.transpose();
-    A_center.transpose();
-    // Single value decompose the data matrix
-    std::auto_ptr<SVD> pSVD( new SVD( A_center ) );
-    pSVD->decompose();
-    utility::matrix U(pSVD->U());
-    utility::matrix V(pSVD->V());
-    // Read the eigenvectors and eigenvalues
-    eigenvectors_.clone(V);
-    eigenvectors_.transpose();
-    eigenvalues_.clone(pSVD->s());
-    // Transform back when done with SVD!
-    // (used V insted of U now for eigenvectors)
-    A_.transpose();
-    A_center.transpose();
-    // T
-    for( size_t i = 0; i < eigenvalues_.size(); ++i )
-      eigenvalues_[ i ] = eigenvalues_[ i ]*eigenvalues_[ i ];
-    eigenvalues_ *= 1.0/A_center.rows();
-    // Sort the eigenvectors in order of eigenvalues
-    // Simple (not efficient) algorithm that always
-    // make sure that the i:th element is in its correct
-    // position (N element --> Ordo( N*N ))
-    for ( size_t i = 0; i < eigenvalues_.size(); ++i )
-      for ( size_t j = i + 1; j < eigenvalues_.size(); ++j )
-        if ( eigenvalues_[ j ] > eigenvalues_[ i ] ) {
-          std::swap( eigenvalues_[ i ], eigenvalues_[ j ] );
-          eigenvectors_.swap_rows(i,j);
+        }
+  }
   */
 …
   void PCA::row_center(utility::matrix& A_center)
+  {
     meanvalues_.clone(utility::vector(A_.rows()));
+    meanvalues_ = vector(A_.rows());
     utility::vector A_row_sum(A_.rows());
+    for (size_t i=0; i<A_row_sum.size(); ++i)
+      A_row_sum(i)=utility::sum(utility::vector(A_,i));
+    for (size_t i=0; i<A_row_sum.size(); ++i){
+      vectorConstView tmp;
+      tmp = A_.const_row(i);
+      A_row_sum(i) = sum(tmp);
+    }
     for( size_t i = 0; i < A_center.rows(); ++i ) {
       meanvalues_[i] = A_row_sum(i) / static_cast<double>(A_.columns());
+      meanvalues_(i) = A_row_sum(i) / static_cast<double>(A_.columns());
       for( size_t j = 0; j < A_center.columns(); ++j )
         A_center(i,j) = A_(i,j) - meanvalues_(i);

branches/peter-dev/yat/utility/matrix.cc

-                      r916
+                      r995
 #include "matrix.h"
 #include "vector.h"
+#include "vectorView.h"
+#include "vectorConstView.h"
 #include "utility.h"
 …
+  void matrix::column(const size_t column, const vector& vec)
+  {
+    assert(m_);
+    int status=gsl_matrix_set_col(m_, column, vec.gsl_vector_p());
+    if (status)
+      throw utility::GSL_error(std::string("matrix::set_column",status));
+  }
+  void matrix::row(const size_t row, const vector& vec)
+  {
+    assert(m_);
+    int status=gsl_matrix_set_row(m_, row, vec.gsl_vector_p());
+    if (status)
+      throw utility::GSL_error(std::string("matrix::set_row",status));
+  vectorView matrix::column_vec(size_t col)
+  {
+    vectorView res(*this, col, false);
+    return res;
+  }
+  vectorConstView matrix::const_column(size_t col) const
+  {
+    return vectorConstView(*this, col, false);
+  }
+  vectorConstView matrix::const_row(size_t col) const
+  {
+    return vectorConstView(*this, col, true);
+  }
+  vectorView matrix::row_vec(size_t row)
+  {
+    vectorView res(*this, row, true);
+    return res;
+  }
 …
   vector operator*(const matrix& m, const vector& v)
+  vector operator*(const matrix& m, const vectorBase& v)
+  {
     utility::vector res(m.rows());
     for (size_t i=0; i<res.size(); ++i)
       res(i) = vector(m,i) * v;
+      res(i) = vectorConstView(m,i) * v;
     return res;
+  }
 …
     utility::vector res(m.columns());
     for (size_t i=0; i<res.size(); ++i)
       res(i) = v * vector(m,i,false);
+      res(i) = v * vectorConstView(m,i,false);
     return res;
+  }

branches/peter-dev/yat/utility/matrix.h

-                      r865
+                      r995
 #include "vector.h"
+#include "vectorView.h"
+#include "vectorConstView.h"
 #include "Exception.h"
 …
     ~matrix(void);
+    ///
+    /// Set all elements to \a value.
+    ///
+    void all(const double value);
     /**
        \brief Make a copy of \a other.
 …
     const matrix& clone(const matrix& other);
+    /**
+     */
+    vectorView column_vec(size_t);
     ///
     /// @return The number of columns in the matrix.
     ///
     size_t columns(void) const;
+    /**
+     */
+    vectorConstView const_column(size_t) const;
+    /**
+     */
+    vectorConstView const_row(size_t) const;
     /**
 …
     size_t rows(void) const;
+    ///
+    /// Set all elements to \a value.
+    ///
+    void all(const double value);
+    /**
+       \brief Set \a column values to values in \a vec.
+       \note No check on size is done.
+       \throw GSL_error if index is out of range or mis-match in
+       sizes.
+    */
+    void column(const size_t column, const vector& vec);
+    /**
+       \brief Set \a row values to values in \a vec.
+       \note No check on size is done.
+       \throw GSL_error if index is out of range or mis-match in
+       sizes.
+    */
+    void row(const size_t row, const vector& vec);
+    /**
+     */
+    vectorView row_vec(size_t);
     /**

branches/peter-dev/yat/utility/vector.cc

-                      r978
+                      r995
   vector::vector(void)
     : v_(NULL), view_(NULL), view_const_(NULL), proxy_v_(NULL)
+    : vectorMutable(NULL)
+  {
+  }
 …
   vector::vector(size_t n, double init_value)
+    : v_(gsl_vector_alloc(n)), view_(NULL), view_const_(NULL),
+      proxy_v_(v_)
+  {
+    if (!v_)
+    : vectorMutable(gsl_vector_alloc(n))
+  {
+    if (!vec_)
       throw utility::GSL_error("vector::vector failed to allocate memory");
 …
   vector::vector(const vector& other)
+    : v_(other.create_gsl_vector_copy()), view_(NULL),
+      view_const_(NULL), proxy_v_(v_)
+  {
+  }
+  vector::vector(vector& v, size_t offset, size_t n, size_t stride)
+    : view_const_(NULL)
+  {
+    view_ = new gsl_vector_view(gsl_vector_subvector_with_stride(v.v_,offset,
+                                                                 stride,n));
+    if (!view_)
+      throw utility::GSL_error("vector::vector failed to setup view");
+    proxy_v_ = v_ = &(view_->vector);
+  }
+  vector::vector(const vector& v, size_t offset, size_t n, size_t stride)
+    : v_(NULL), view_(NULL)
+  {
+    view_const_ = new gsl_vector_const_view(
+                   gsl_vector_const_subvector_with_stride(v.v_,offset,stride,n));
+    if (!view_const_)
+      throw utility::GSL_error("vector::vector failed to setup view");
+    proxy_v_ = &(view_const_->vector);
+  }
+  vector::vector(matrix& m, size_t i, bool row)
+    : view_const_(NULL)
+  {
+    view_=new gsl_vector_view(row ?
+                              gsl_matrix_row   (m.gsl_matrix_p(),i) :
+                              gsl_matrix_column(m.gsl_matrix_p(),i)  );
+    if (!view_)
+      throw utility::GSL_error("vector::vector failed to setup view");
+    proxy_v_ = v_ = &(view_->vector);
+  }
+  vector::vector(const matrix& m, size_t i, bool row)
+    : v_(NULL), view_(NULL)
+  {
+    view_const_ = new gsl_vector_const_view(row ?
+                                   gsl_matrix_const_row   (m.gsl_matrix_p(),i) :
+                                   gsl_matrix_const_column(m.gsl_matrix_p(),i) );
+    if (!view_const_)
+      throw utility::GSL_error("vector::vector failed to setup view");
+    proxy_v_ = &(view_const_->vector);
+    : vectorMutable(create_gsl_vector_copy(other))
+  {
+  }
+  vector::vector(const vectorBase& other)
+    : vectorMutable(create_gsl_vector_copy(other))
+  {
+  }
 …
   vector::vector(std::istream& is, char sep)
     throw (utility::IO_error, std::exception)
     : view_(NULL), view_const_(NULL)
+    : vectorMutable(NULL)
+  {
     // read the data file and store in stl vectors (dynamically
 …
     is.clear(std::ios::goodbit);
     // convert the data to a gsl vector
     proxy_v_ = v_ = gsl_vector_alloc(nof_rows*nof_columns);
     if (!v_)
+    vec_ = gsl_vector_alloc(nof_rows*nof_columns);
+    if (!vec_)
       throw utility::GSL_error("vector::vector failed to allocate memory");
     size_t n=0;
 …
     for (size_t i=0; i<nof_rows; i++)
       for (size_t j=0; j<nof_columns; j++)
         gsl_vector_set( v_, n++, data_matrix[i][j] );
+        gsl_vector_set( vec_, n++, data_matrix[i][j] );
+  }
 …
   vector::~vector(void)
+  {
     delete_allocated_memory();
+  }
   vector::iterator vector::begin(void)
+  {
     return vector::iterator(*this, 0);
+  }
   vector::const_iterator vector::begin(void) const
+  {
     return vector::const_iterator(*this, 0);
+  }
   const vector& vector::clone(const vector& other)
+  {
     if (this!=&other) {
       delete_allocated_memory();
       if (other.view_) {
+        view_ = new gsl_vector_view(*other.view_);
+        proxy_v_ = v_ = &(view_->vector);
+      }
       else if (other.view_const_) {
         view_const_ = new gsl_vector_const_view(*other.view_const_);
         proxy_v_ = &(view_const_->vector);
         v_=NULL;
+      }
       else if (other.v_)
         proxy_v_ = v_ = other.create_gsl_vector_copy();
+    gsl_vector_free(vec_);
+  }
+  gsl_vector* vector::create_gsl_vector_copy(const vectorBase& other) const
+  {
+    gsl_vector* vec = gsl_vector_alloc(other.size());
+    if (!vec)
+      throw utility::GSL_error("vector::create_gsl_vector_copy failed to allocate memory");
+    if (gsl_vector_memcpy(vec, other.gsl_vector_p()))
+      throw utility::GSL_error("vector::create_gsl_matrix_copy dimension mis-match");
+    return vec;
+  }
+  void vector::resize(size_t n, double init_value)
+  {
+    gsl_vector_free(vec_);
+    vec_ = gsl_vector_alloc(n);
+    if (!vec_)
+      throw utility::GSL_error("vector::vector failed to allocate memory");
+    all(init_value);
+  }
+  //Peter use swap idiom
+  const vector& vector::operator=( const vectorBase& other )
+  {
+    if (!other.size())
+      vec_=NULL;
+    else  {
+      if (size()!=other.size())
+        resize(other.size(),0.0);
+      gsl_vector_memcpy(vec_,other.gsl_vector_p());
+    }
     return *this;
 …
+  gsl_vector* vector::create_gsl_vector_copy(void) const
+  {
+    gsl_vector* vec = gsl_vector_alloc(size());
+    if (!vec)
+      throw utility::GSL_error("vector::create_gsl_vector_copy failed to allocate memory");
+    if (gsl_vector_memcpy(vec, proxy_v_))
+      throw utility::GSL_error("vector::create_gsl_matrix_copy dimension mis-match");
+    return vec;
+  }
+  void vector::delete_allocated_memory(void)
+  {
+    if (view_)
+      delete view_;
+    else if (view_const_)
+      delete view_const_;
+    else if (v_)
+      gsl_vector_free(v_);
+    proxy_v_=v_=NULL;
+  }
+  void vector::div(const vector& other)
+  {
+    assert(v_);
+    int status=gsl_vector_div(v_,other.gsl_vector_p());
+    if (status)
+      throw utility::GSL_error(std::string("vector::div",status));
+  }
+  vector::iterator vector::end(void)
+  {
+    return vector::iterator(*this, size());
+  }
+  vector::const_iterator vector::end(void) const
+  {
+    return vector::const_iterator(*this, size());
+  }
+  bool vector::equal(const vector& other, const double d) const
+  {
+    if (this==&other)
+      return true;
+    if (size()!=other.size())
+      return false;
+    // if gsl error handler disabled, out of bounds index will not
+    // abort the program.
+    for (size_t i=0; i<size(); ++i)
+      // The two last condition checks are needed for NaN detection
+      if (fabs( (*this)(i)-other(i) ) > d ||
+          (*this)(i)!=(*this)(i) || other(i)!=other(i))
+        return false;
+    return true;
+  }
+  const gsl_vector* vector::gsl_vector_p(void) const
+  {
+    return proxy_v_;
+  }
+  gsl_vector* vector::gsl_vector_p(void)
+  {
+    return v_;
+  }
+  bool vector::isview(void) const
+  {
+    return view_ || view_const_;
+  }
+  void vector::mul(const vector& other)
+  {
+    assert(v_);
+    int status=gsl_vector_mul(v_,other.gsl_vector_p());
+    if (status)
+      throw utility::GSL_error(std::string("vector::mul",status));
+  }
+  void vector::resize(size_t n, double init_value)
+  {
+    delete_allocated_memory();
+    proxy_v_ = v_ = gsl_vector_alloc(n);
+    if (!v_)
+      throw utility::GSL_error("vector::vector failed to allocate memory");
+    all(init_value);
+  }
+  void vector::reverse(void)
+  {
+    assert(v_);
+    gsl_vector_reverse(v_);
+  }
+  void vector::all(const double& value)
+  {
+    assert(v_);
+    gsl_vector_set_all(v_,value);
+  }
+  size_t vector::size(void) const
+  {
+    if (!proxy_v_)
+      return 0;
+    return proxy_v_->size;
+  }
+  void vector::swap(size_t i, size_t j)
+  {
+    assert(v_);
+    int status=gsl_vector_swap_elements(v_, i, j);
+    if (status)
+      throw utility::GSL_error(std::string("vector::swap_elements",status));
+  }
+  double& vector::operator()(size_t i)
+  {
+    assert(v_);
+    double* d=gsl_vector_ptr(v_, i);
+    if (!d)
+      throw utility::GSL_error("vector::operator()",GSL_EINVAL);
+    return *d;
+  }
+  const double& vector::operator()(size_t i) const
+  {
+    const double* d=gsl_vector_const_ptr(proxy_v_, i);
+    if (!d)
+      throw utility::GSL_error("vector::operator()",GSL_EINVAL);
+    return *d;
+  }
+  double& vector::operator[](size_t i)
+  {
+    return this->operator()(i);
+  }
+  const double& vector::operator[](size_t i) const
+  {
+    return this->operator()(i);
+  }
+  bool vector::operator==(const vector& other) const
+  {
+    return equal(other);
+  }
+  bool vector::operator!=(const vector& other) const
+  {
+    return !equal(other);
+  }
+  double vector::operator*( const vector &other ) const
+  {
+    double res = 0.0;;
+    for ( size_t i = 0; i < size(); ++i )
+      res += other(i) * (*this)(i);
+    return res;
+  }
+  const vector& vector::operator=( const vector& other )
+  {
+    assert(v_);
+    if (gsl_vector_memcpy(v_,other.gsl_vector_p()))
+      throw utility::GSL_error("vector::set dimension mis-match");
+  //Peter use swap idiom
+  const vectorMutable& vector::operator=( vectorMutable& other )
+  {
+    if (!other.size())
+      vec_=NULL;
+    else  {
+      if (size()!=other.size())
+        resize(other.size(),0.0);
+      gsl_vector_memcpy(vec_,other.gsl_vector_p());
+    }
     return *this;
+  }
-  const vector& vector::operator+=(const vector& other)
+  {
-    assert(v_);
-    int status=gsl_vector_add(v_, other.gsl_vector_p());
-    if (status)
-      throw utility::GSL_error(std::string("vector::add", status));
-    return *this;
+  }
-  const vector& vector::operator+=(double d)
+  {
-    assert(v_);
-    gsl_vector_add_constant(v_, d);
-    return *this;
+  }
-  const vector& vector::operator-=(const vector& other)
+  {
-    assert(v_);
-    int status=gsl_vector_sub(v_, other.gsl_vector_p());
-    if (status)
-      throw utility::GSL_error(std::string("vector::sub", status));
-    return *this;
+  }
-  const vector& vector::operator-=(const double d)
+  {
-    assert(v_);
-    gsl_vector_add_constant(v_, -d);
-    return *this;
+  }
-  const vector& vector::operator*=(const double d)
+  {
-    assert(v_);
-    gsl_vector_scale(v_, d);
-    return *this;
+  }
-  bool isnull(const vector& v)
+  {
-    return gsl_vector_isnull(v.gsl_vector_p());
+  }
-  double max(const vector& v)
+  {
-    return gsl_vector_max(v.gsl_vector_p());
+  }
-  size_t max_index(const vector& v)
+  {
-    return gsl_vector_max_index(v.gsl_vector_p());
+  }
-  double min(const vector& v)
+  {
-    return gsl_vector_min(v.gsl_vector_p());
+  }
-  size_t min_index(const vector& v)
+  {
-    return gsl_vector_min_index(v.gsl_vector_p());
+  }
-  bool nan(const vector& templat, vector& flag)
+  {
-    size_t vsize=templat.size();
-    if (vsize!=flag.size())
-      flag.clone(vector(vsize,1.0));
-    else
-      flag.all(1.0);
-    bool nan=false;
-    for (size_t i=0; i<vsize; i++)
-      if (std::isnan(templat(i))) {
-        flag(i)=0;
-        nan=true;
+      }
-    return nan;
+  }
-  void set_basis(vector& v, size_t i)
+  {
-    assert(v.gsl_vector_p());
-    gsl_vector_set_basis(v.gsl_vector_p(),i);
+  }
-  void shuffle(vector& invec)
+  {
-    random::DiscreteUniform rnd;
-    std::random_shuffle(invec.begin(), invec.end(), rnd);
+  }
-  void sort(vector& v)
+  {
-    assert(v.gsl_vector_p());
-    std::sort(v.begin(), v.end());
+  }
-  void sort_index(std::vector<size_t>& sort_index, const vector& invec)
+  {
-    assert(invec.gsl_vector_p());
-    gsl_permutation* p = gsl_permutation_alloc(invec.size());
-    int status=gsl_sort_vector_index(p,invec.gsl_vector_p());
-    if (status) {
-      gsl_permutation_free(p);
-      throw utility::GSL_error(std::string("sort_index(vector&,const vector&)",status));
+    }
-    sort_index=std::vector<size_t>(p->data,p->data+p->size);
-    gsl_permutation_free(p);
+  }
-  void sort_smallest_index(std::vector<size_t>& sort_index, size_t k,
-                            const vector& invec)
+  {
-    assert(invec.gsl_vector_p());
-    assert(k<=invec.size());
-    sort_index.resize(k);
-    gsl_sort_vector_smallest_index(&sort_index[0],k,invec.gsl_vector_p());
+  }
-  void sort_largest_index(std::vector<size_t>& sort_index, size_t k,
-                            const vector& invec)
+  {
-    assert(invec.gsl_vector_p());
-    assert(k<=invec.size());
-    sort_index.resize(k);
-    gsl_sort_vector_largest_index(&sort_index[0],k,invec.gsl_vector_p());
+  }
-  double sum(const vector& v)
+  {
-    double sum = 0;
-    size_t vsize=v.size();
-    for (size_t i=0; i<vsize; ++i)
-      sum += v[i];
-    return sum;
+  }
 …
-  std::ostream& operator<<(std::ostream& s, const vector& a)
+  {
-    s.setf(std::ios::dec);
-    s.precision(12);
-    for (size_t j = 0; j < a.size(); ++j) {
-      s << a[j];
-      if ( (j+1)<a.size() )
-        s << s.fill();
+    }
-    return s;
+  }
 }}} // of namespace utility, yat, and thep

branches/peter-dev/yat/utility/vector.h

-                      r966
+                      r995
 */
+#include "vectorMutable.h"
 #include "Exception.h"
 #include "Iterator.h"
 …
   */
   class vector
+  class vector : public vectorMutable
+  {
   public:
-    /// \brief vector::iterator
-    typedef Iterator<double&, vector> iterator;
-    /// \brief vector::const_iterator
-    typedef Iterator<const double, const vector> const_iterator;
     /**
 …
     /**
        \brief Vector view constructor.
+       \brief The copy constructor.
        Create a view of vector \a v, with starting index \a offset,
        size \a n, and an optional \a stride.
+       \note If the object to be copied is a vector view, the values
+       of the view will be copied, i.e. the view is not copied.
+       A vector view can be used as any vector with the difference
+       that changes made to the view will also change the object that
+       is viewed. Also, using the copy constructor will create a new
+       vector object that is a copy of whatever is viewed. If a copy
+       of the view is needed then you should use this constructor to
+       obtain a copy.
+       \note If the object viewed by the view goes out of scope or is
+       deleted, the view becomes invalid and the result of further use
+       is undefined.
+       \throw GSL_error if a view cannot be set up.
+       \throw A GSL_error is indirectly thrown if memory allocation
+       fails.
     */
+    vector(vector& v, size_t offset, size_t n, size_t stride=1);
+    /**
+       \brief Vector const view constructor.
+       Create a view of vector \a v, with starting index \a offset,
+       size \a n, and an optional \a stride.
+       A vector view can be used as any const vector. Using the copy
+       constructor will create a new vector object that is a copy of
+       whatever is viewed. If a copy of the view is needed then you
+       should use this constructor to obtain a copy.
+       \note If the object viewed by the view goes out of scope or is
+       deleted, the view becomes invalid and the result of further use
+       is undefined.
+       \throw GSL_error if a view cannot be set up.
+    */
+    vector(const vector& v, size_t offset, size_t n, size_t stride=1);
+    ///
+    /// Matrix row/column view constructor.
+    ///
+    /// Create a row/column vector view of matrix \a m, pointing at
+    /// row/column \a i. The parameter \a row is used to set whether
+    /// the view should be a row or column view. If \a row is set to
+    /// true, the view will be a row view (default behaviour), and,
+    /// naturally, a column view otherwise.
+    ///
+    /// A vector view can be used as any vector with the difference
+    /// that changes made to the view will also change the object that
+    /// is viewed. Also, using the copy constructor will create a new
+    /// vector object that is a copy of whatever is viewed. If a copy
+    /// of the view is needed then you should use the vector view
+    /// constructor to obtain a copy.
+    ///
+    /// @note If the object viewed by the view goes out of scope or is
+    /// deleted, the view becomes invalid and the result of further
+    /// use is undefined.
+    ///
+    vector(matrix& m, size_t i, bool row=true);
+    ///
+    /// Matrix row/column const view constructor.
+    ///
+    /// Create a row/column vector view of matrix \a m, pointing at
+    /// row/column \a i. The parameter \a row is used to set whether
+    /// the view should be a row or column view. If \a row is set to
+    /// true, the view will be a row view (default behaviour), and,
+    /// naturally, a column view otherwise.
+    ///
+    /// A const vector view can be used as any const vector. Using the
+    /// copy constructor will create a new vector object that is a
+    /// copy of whatever is viewed. If a copy of the view is needed
+    /// then you should use the vector view constructor to obtain a
+    /// copy.
+    ///
+    /// @note If the object viewed by the view goes out of scope or is
+    /// deleted, the view becomes invalid and the result of further
+    /// use is undefined.
+    ///
+    vector(const matrix& m, size_t i, bool row=true);
+    vector(const vectorBase& other);
     /**
 …
     /**
-       \return mutable iterator to start of vector
-     */
-    iterator begin(void);
-    /**
-       \return read-only iterator to start of vector
-     */
-    const_iterator begin(void) const;
-    /**
-       \brief Make a copy of \a other.
-       This function will make a deep copy of \a other. Memory is
-       resized and view state is changed if needed.
-    */
-    const vector& clone(const vector& other);
-    /**
-       \brief This function performs element-wise division, \f$ this_i =
-       this_i/other_i \; \forall i \f$.
-       \throw GSL_error if dimensions mis-match.
-    */
-    void div(const vector& other);
-    /**
-       \return mutable iterator to end of vector
-     */
-    iterator end(void);
-    /**
-       \return read-only iterator to end of vector
-     */
-    const_iterator end(void) const;
-    /**
-       \brief Check whether vectors are equal within a user defined
-       precision, set by \a precision.
-       \return True if each element deviates less or equal than \a
-       d. If any vector contain a NaN, false is always returned.
-       \see operator== and operator!=
-    */
-    bool equal(const vector&, const double precision=0) const;
-    ///
-    /// @return A const pointer to the internal GSL vector,
-    ///
-    const gsl_vector* gsl_vector_p(void) const;
-    ///
-    /// @return A pointer to the internal GSL vector,
-    ///
-    gsl_vector* gsl_vector_p(void);
-    ///
-    /// Check if the vector object is a view (sub-vector) to another
-    /// vector.
-    ///
-    /// @return True if the object is a view, false othwerwise.
-    ///
-    bool isview(void) const;
-    /**
-       \brief This function performs element-wise multiplication, \f$
-       this_i = this_i * other_i \; \forall i \f$.
-       \throw GSL_error if dimensions mis-match.
-    */
-    void mul(const vector& other);
-    /**
        @brief Resize vector
 …
        vector becomes invalid.
     */
+    void resize(size_t, double init_value=0);
+    /**
+       \brief Reverse the order of elements in the vector.
+    */
+    void reverse(void);
+    ///
+    /// Set all elements to \a value.
+    ///
+    void all(const double& value);
+    ///
+    /// @return the number of elements in the vector.
+    ///
+    size_t size(void) const;
+    /**
+       \brief Exchange elements \a i and \a j.
+       \throw GSL_error if vector lengths differs.
+    */
+    void swap(size_t i, size_t j);
+    /**
+       \brief Element access operator.
+       \return Reference to element \a i.
+       \throw If GSL range checks are enabled in the underlying GSL
+       library a GSL_error exception is thrown if either index is out
+       of range.
+    */
+    double& operator()(size_t i);
+    /**
+       \brief Element access operator.
+       \return Const reference to element \a i.
+       \throw If GSL range checks are enabled in the underlying GSL
+       library a GSL_error exception is thrown if either index is out
+       of range.
+    */
+    const double& operator()(size_t i) const;
+    ///
+    /// Element access operator.
+    ///
+    /// @return Reference to element \a i.
+    ///
+    double& operator[](size_t i);
+    ///
+    /// Const element access operator.
+    ///
+    /// @return The value of element \a i.
+    ///
+    const double& operator[](size_t i) const;
+    /**
+       \brief Comparison operator. Takes linear time.
+       Checks are performed with exact matching, i.e., rounding off
+       effects may destroy comparison. Use the equal function for
+       comparing elements within a user defined precision.
+       \return True if all elements are equal otherwise false.
+       \see equal(const vector&, const double precision=0)
+    */
+    bool operator==(const vector&) const;
+    /**
+       \brief Comparison operator. Takes linear time.
+       Checks are performed with exact matching, i.e., rounding off
+       effects may destroy comparison. Use the equal function for
+       comparing elements within a user defined precision.
+       \return False if all elements are equal otherwise true.
+       \see equal(const vector&, const double precision=0)
+    */
+    bool operator!=(const vector&) const;
+    ///
+    /// @return The dot product.
+    ///
+    double operator*(const vector&) const;
+    void resize(size_t n, double init_value);
     /**
        \brief The assignment operator.
-       Dimensions of the vectors must match. If the LHS vector is a
-       view, the underlying data will be changed.
-       \return A const reference to the resulting vector.
-       \see void set(const vector&).
-       \throw GSL_error if dimensions mis-match.
-    */
-    const vector& operator=(const vector&);
-    /**
-       \brief Addition and assign operator. Vector addition, \f$
-       this_i = this_i + other_i \; \forall i \f$.
-       \return A const reference to the resulting vector.
-       \throw GSL_error if dimensions mis-match.
-    */
-    const vector& operator+=(const vector&);
-    /**
-       \brief Add a constant to a vector, \f$ this_i = this_i + d \;
-       \forall i \f$.
        \return A const reference to the resulting vector.
     */
+    const vector& operator+=(double d);
+    /**
+       \brief Subtract and assign operator. Vector subtraction, \f$
+       this_i = this_i - other_i \; \forall i \f$.
+       \return A const reference to the resulting vector.
+       \throw GSL_error if dimensions mis-match.
+    */
+    const vector& operator-=(const vector&);
+    /**
+       \brief Subtract a constant to a vector, \f$ this_i = this_i - d
+       \; \forall i \f$.
+       \return A const reference to the resulting vector.
+    */
+    const vector& operator-=(double d);
+    /**
+       \brief Multiply with scalar and assign operator, \f$ this_i =
+       this_i * d \; \forall i \f$.
+       \return A const reference to the resulting vector.
+    */
+    const vector& operator*=(double d);
+    const vectorMutable& operator=(vectorMutable&);
+    const vector& operator=(const vectorBase&);
   private:
 …
        copy, or if dimensions mis-match.
     */
     gsl_vector* create_gsl_vector_copy(void) const;
+    gsl_vector* create_gsl_vector_copy(const vectorBase&) const;
-    /**
-       \brief Clear all dynamically allocated memory.
-       Internal utility function.
-    */
-    void delete_allocated_memory(void);
-    gsl_vector* v_;
-    gsl_vector_view* view_;
-    const gsl_vector_const_view* view_const_;
-    // proxy_v_ is used to access the proper underlying gsl_vector
-    // in all const member functions. It is not used by design for
-    // non-const vector functions and operators. This is to make sure
-    // that runtime errors occur if a const vector is used in an
-    // inappropriate manner such as on left hand side in assignment
-    // (remember, v_ is null for const vector views).
-    const gsl_vector* proxy_v_;
   };
-  /**
-     \brief Check if all elements of the vector are zero.
-     \return True if all elements in the vector is zero, false
-     othwerwise.
-  */
-  bool isnull(const vector&);
-  /**
-     \brief Get the maximum value of the vector.
-     \return The maximum value of the vector.
-  */
-  double max(const vector&);
-  /**
-     \brief Locate the maximum value in the vector.
-     \return The index to the maximum value of the vector.
-     \note Lower index has precedence.
-  */
-  size_t max_index(const vector&);
-  /**
-     \brief Get the minimum value of the vector.
-     \return The minimum value of the vector.
-  */
-  double min(const vector&);
-  /**
-     \brief Locate the minimum value in the vector.
-     \return The index to the minimum value of the vector.
-     \note Lower index has precedence.
-  */
-  size_t min_index(const vector&);
-  /**
-     \brief Create a vector \a flag indicating NaN's in another vector
-     \a templat.
-     The \a flag vector is changed to contain 1's and 0's only. A 1
-     means that the corresponding element in the \a templat vector is
-     valid and a zero means that the corresponding element is a NaN.
-     \note Space for vector \a flag is reallocated to fit the size of
-     vector \a templat if sizes mismatch.
-     \return True if the \a templat vector contains at least one NaN.
-  */
-  bool nan(const vector& templat, vector& flag);
-  /**
-     \brief Transforms a vector to a basis vector.
-     All elements are set to zero except the \a i-th element which is
-     set to one.
-  */
-  void set_basis(vector&, size_t i);
-  /**
-     Randomly shuffles the elements in vector \a invec
-  */
-  void shuffle(vector& invec);
-  /**
-     Sort the elements in the vector.
-  */
-  void sort(vector&);
-  /**
-     Create a vector \a sort_index containing the indeces of
-     elements in a another vector \a invec.  The elements of \a
-     sort_index give the index of the vector element which would
-     have been stored in that position if the vector had been sorted
-     in place. The first element of \a sort_index gives the index of the least
-     element in \a invec, and the last element of \a sort_index gives the index of the
-     greatest element in \a invec . The vector \a invec is not changed.
-  */
-  void sort_index(std::vector<size_t>& sort_index, const vector& invec);
-  /** Similar to sort_index but creates a vector with indices to the \a k
-  smallest elements in \a invec.
-  */
-  void sort_smallest_index(std::vector<size_t>& sort_index, size_t k, const
-  vector& invec);
-  /** Similar to sort_index but creates a vector with indices to the \a k
-  largest elements in \a invec.
-  */
-  void sort_largest_index(std::vector<size_t>& sort_index, size_t k, const
-  vector& invec);
-  /**
-     \brief Calculate the sum of all vector elements.
-     \return The sum.
-  */
-  double sum(const vector&);
   /**
 …
   void swap(vector&, vector&);
-  /**
-     \brief The output operator for the vector class.
-  */
-  std::ostream& operator<<(std::ostream&, const vector&);
 }}} // of namespace utility, yat, and theplu

branches/peter-dev/yat/utility/vectorBase.cc

-                      r994
+                      r995
 */
 #include "vector.h"
+#include "vectorBase.h"
 #include "matrix.h"
 #include "utility.h"
 …
+  vector::vector(void)
+    : v_(NULL), view_(NULL), view_const_(NULL), proxy_v_(NULL)
+  {
+  }
+  vector::vector(size_t n, double init_value)
+    : v_(gsl_vector_alloc(n)), view_(NULL), view_const_(NULL),
+      proxy_v_(v_)
+  {
+    if (!v_)
+      throw utility::GSL_error("vector::vector failed to allocate memory");
+    all(init_value);
+  }
+  vector::vector(const vector& other)
+    : v_(other.create_gsl_vector_copy()), view_(NULL),
+      view_const_(NULL), proxy_v_(v_)
+  {
+  }
+  vector::vector(vector& v, size_t offset, size_t n, size_t stride)
+    : view_const_(NULL)
+  {
+    view_ = new gsl_vector_view(gsl_vector_subvector_with_stride(v.v_,offset,
+                                                                 stride,n));
+    if (!view_)
+      throw utility::GSL_error("vector::vector failed to setup view");
+    proxy_v_ = v_ = &(view_->vector);
+  }
+  vector::vector(const vector& v, size_t offset, size_t n, size_t stride)
+    : v_(NULL), view_(NULL)
+  {
+    view_const_ = new gsl_vector_const_view(
+                   gsl_vector_const_subvector_with_stride(v.v_,offset,stride,n));
+    if (!view_const_)
+      throw utility::GSL_error("vector::vector failed to setup view");
+    proxy_v_ = &(view_const_->vector);
+  }
+  vector::vector(matrix& m, size_t i, bool row)
+    : view_const_(NULL)
+  {
+    view_=new gsl_vector_view(row ?
+                              gsl_matrix_row   (m.gsl_matrix_p(),i) :
+                              gsl_matrix_column(m.gsl_matrix_p(),i)  );
+    if (!view_)
+      throw utility::GSL_error("vector::vector failed to setup view");
+    proxy_v_ = v_ = &(view_->vector);
+  }
+  vector::vector(const matrix& m, size_t i, bool row)
+    : v_(NULL), view_(NULL)
+  {
+    view_const_ = new gsl_vector_const_view(row ?
+                                   gsl_matrix_const_row   (m.gsl_matrix_p(),i) :
+                                   gsl_matrix_const_column(m.gsl_matrix_p(),i) );
+    if (!view_const_)
+      throw utility::GSL_error("vector::vector failed to setup view");
+    proxy_v_ = &(view_const_->vector);
+  }
+  vector::vector(std::istream& is, char sep)
+    throw (utility::IO_error, std::exception)
+    : view_(NULL), view_const_(NULL)
+  {
+    // read the data file and store in stl vectors (dynamically
+    // expandable)
+    std::vector<std::vector<double> > data_matrix;
+    u_int nof_columns=0;
+    u_int nof_rows=0;
+    std::string line;
+    while(getline(is, line, '\n')) {
+      // Empty lines
+      if (!line.size())
+          continue;
+      nof_rows++;
+      std::vector<double> v;
+      std::string element;
+      std::stringstream ss(line);
+      bool ok=true;
+      while(ok) {
+        if(sep=='\0')
+          ok=(ss>>element);
+        else
+          ok=getline(ss, element, sep);
+        if(!ok)
+          break;
+        if(utility::is_double(element)) {
+          v.push_back(atof(element.c_str()));
+        }
+        else if (!element.size() || utility::is_nan(element)) {
+          v.push_back(std::numeric_limits<double>::quiet_NaN());
+        }
+        else {
+          std::stringstream ss("Warning: '");
+          ss << element << "' is not a double.";
+          throw IO_error(ss.str());
+        }
+      }
+      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
+          v.push_back(std::numeric_limits<double>::quiet_NaN());
+      if (!nof_columns)
+        nof_columns=v.size();
+      else if (nof_rows && (nof_columns>1)) {
+        std::ostringstream s;
+        s << "vector::vector(std::istream&) data file error:\n"
+          << "    File has inconsistent number of rows (" << nof_rows
+          << ") and columns (" << nof_columns
+          << ").\n    Expected a row or a column vector.";
+        throw utility::IO_error(s.str());
+      }
+      else if (v.size()!=nof_columns) {
+        std::ostringstream s;
+        s << "vector::vector(std::istream&) data file error:\n"
+          << "    Line " << nof_rows << " has " << v.size()
+          << " columns; expected " << nof_columns << " column.";
+        throw utility::IO_error(s.str());
+      }
+      data_matrix.push_back(v);
+    }
+    // manipulate the state of the stream to be good
+    is.clear(std::ios::goodbit);
+    // convert the data to a gsl vector
+    proxy_v_ = v_ = gsl_vector_alloc(nof_rows*nof_columns);
+    if (!v_)
+      throw utility::GSL_error("vector::vector failed to allocate memory");
+    size_t n=0;
+    // if gsl error handler disabled, out of bounds index will not
+    // abort the program.
+    for (size_t i=0; i<nof_rows; i++)
+      for (size_t j=0; j<nof_columns; j++)
+        gsl_vector_set( v_, n++, data_matrix[i][j] );
+  }
+  vector::~vector(void)
+  {
+    delete_allocated_memory();
+  }
+  vector::iterator vector::begin(void)
+  {
+    return vector::iterator(*this, 0);
+  }
+  vector::const_iterator vector::begin(void) const
+  {
+    return vector::const_iterator(*this, 0);
+  }
+  const vector& vector::clone(const vector& other)
+  {
+    if (this!=&other) {
+      delete_allocated_memory();
+      if (other.view_) {
+        view_ = new gsl_vector_view(*other.view_);
+        proxy_v_ = v_ = &(view_->vector);
+      }
+      else if (other.view_const_) {
+        view_const_ = new gsl_vector_const_view(*other.view_const_);
+        proxy_v_ = &(view_const_->vector);
+        v_=NULL;
+      }
+      else if (other.v_)
+        proxy_v_ = v_ = other.create_gsl_vector_copy();
+    }
+    return *this;
+  }
+  gsl_vector* vector::create_gsl_vector_copy(void) const
+  {
+    gsl_vector* vec = gsl_vector_alloc(size());
+    if (!vec)
+      throw utility::GSL_error("vector::create_gsl_vector_copy failed to allocate memory");
+    if (gsl_vector_memcpy(vec, proxy_v_))
+      throw utility::GSL_error("vector::create_gsl_matrix_copy dimension mis-match");
+    return vec;
+  }
+  void vector::delete_allocated_memory(void)
+  {
+    if (view_)
+      delete view_;
+    else if (view_const_)
+      delete view_const_;
+    else if (v_)
+      gsl_vector_free(v_);
+    proxy_v_=v_=NULL;
+  }
+  void vector::div(const vector& other)
+  {
+    assert(v_);
+    int status=gsl_vector_div(v_,other.gsl_vector_p());
+    if (status)
+      throw utility::GSL_error(std::string("vector::div",status));
+  }
+  vector::iterator vector::end(void)
+  {
+    return vector::iterator(*this, size());
+  }
+  vector::const_iterator vector::end(void) const
+  {
+    return vector::const_iterator(*this, size());
+  }
+  bool vector::equal(const vector& other, const double d) const
+  vectorBase::vectorBase(const gsl_vector* v)
+    : const_vec_(v)
+  {
+  }
+  vectorBase::~vectorBase(void)
+  {
+  }
+  vectorBase::const_iterator vectorBase::begin(void) const
+  {
+    return vectorBase::const_iterator(*this, 0);
+  }
+  vectorBase::const_iterator vectorBase::end(void) const
+  {
+    return vectorBase::const_iterator(*this, size());
+  }
+  bool vectorBase::equal(const vectorBase& other, const double d) const
+  {
     if (this==&other)
 …
     // abort the program.
     for (size_t i=0; i<size(); ++i)
-      // The two last condition checks are needed for NaN detection
       if (fabs( (*this)(i)-other(i) ) > d ||
           (*this)(i)!=(*this)(i) || other(i)!=other(i))
+          std::isnan((*this)(i)) || std::isnan(other(i)) )
         return false;
     return true;
 …
+  const gsl_vector* vector::gsl_vector_p(void) const
+  {
+    return proxy_v_;
+  }
+  gsl_vector* vector::gsl_vector_p(void)
+  {
+    return v_;
+  }
+  bool vector::isview(void) const
+  {
+    return view_ || view_const_;
+  }
+  void vector::mul(const vector& other)
+  {
+    assert(v_);
+    int status=gsl_vector_mul(v_,other.gsl_vector_p());
+    if (status)
+      throw utility::GSL_error(std::string("vector::mul",status));
+  }
+  void vector::resize(size_t n, double init_value)
+  {
+    delete_allocated_memory();
+    proxy_v_ = v_ = gsl_vector_alloc(n);
+    if (!v_)
+      throw utility::GSL_error("vector::vector failed to allocate memory");
+    all(init_value);
+  }
+  void vector::reverse(void)
+  {
+    assert(v_);
+    gsl_vector_reverse(v_);
+  }
+  void vector::all(const double& value)
+  {
+    assert(v_);
+    gsl_vector_set_all(v_,value);
+  }
+  size_t vector::size(void) const
+  {
+    if (!proxy_v_)
+  const gsl_vector* vectorBase::gsl_vector_p(void) const
+  {
+    return const_vec_;
+  }
+  size_t vectorBase::size(void) const
+  {
+    if (const_vec_)
       return 0;
+    return proxy_v_->size;
+  }
+  void vector::swap(size_t i, size_t j)
+  {
+    assert(v_);
+    int status=gsl_vector_swap_elements(v_, i, j);
+    if (status)
+      throw utility::GSL_error(std::string("vector::swap_elements",status));
+  }
+  double& vector::operator()(size_t i)
+  {
+    assert(v_);
+    double* d=gsl_vector_ptr(v_, i);
+    return const_vec_->size;
+  }
+  const double& vectorBase::operator()(size_t i) const
+  {
+    const double* d=gsl_vector_const_ptr(const_vec_, i);
     if (!d)
       throw utility::GSL_error("vector::operator()",GSL_EINVAL);
+      throw utility::GSL_error("vectorBase::operator()",GSL_EINVAL);
     return *d;
+  }
   const double& vector::operator()(size_t i) const
+  {
     const double* d=gsl_vector_const_ptr(proxy_v_, i);
+  const double& vectorBase::operator[](size_t i) const
+  {
+    const double* d=gsl_vector_const_ptr(const_vec_, i);
     if (!d)
       throw utility::GSL_error("vector::operator()",GSL_EINVAL);
+      throw utility::GSL_error("vectorBase::operator[]",GSL_EINVAL);
     return *d;
+  }
+  double& vector::operator[](size_t i)
+  {
+    return this->operator()(i);
+  }
+  const double& vector::operator[](size_t i) const
+  {
+    return this->operator()(i);
+  }
+  bool vector::operator==(const vector& other) const
+  bool vectorBase::operator==(const vectorBase& other) const
+  {
     return equal(other);
 …
   bool vector::operator!=(const vector& other) const
+  bool vectorBase::operator!=(const vectorBase& other) const
+  {
     return !equal(other);
 …
   double vector::operator*( const vector &other ) const
+  double vectorBase::operator*( const vectorBase &other ) const
+  {
     double res = 0.0;;
 …
+  const vector& vector::operator=( const vector& other )
+  {
+    assert(v_);
+    if (gsl_vector_memcpy(v_,other.gsl_vector_p()))
+      throw utility::GSL_error("vector::set dimension mis-match");
+    return *this;
+  }
+  const vector& vector::operator+=(const vector& other)
+  {
+    assert(v_);
+    int status=gsl_vector_add(v_, other.gsl_vector_p());
+    if (status)
+      throw utility::GSL_error(std::string("vector::add", status));
+    return *this;
+  }
+  const vector& vector::operator+=(double d)
+  {
+    assert(v_);
+    gsl_vector_add_constant(v_, d);
+    return *this;
+  }
+  const vector& vector::operator-=(const vector& other)
+  {
+    assert(v_);
+    int status=gsl_vector_sub(v_, other.gsl_vector_p());
+    if (status)
+      throw utility::GSL_error(std::string("vector::sub", status));
+    return *this;
+  }
+  const vector& vector::operator-=(const double d)
+  {
+    assert(v_);
+    gsl_vector_add_constant(v_, -d);
+    return *this;
+  }
+  const vector& vector::operator*=(const double d)
+  {
+    assert(v_);
+    gsl_vector_scale(v_, d);
+    return *this;
+  }
+  bool isnull(const vector& v)
+  bool isnull(const vectorBase& v)
+  {
     return gsl_vector_isnull(v.gsl_vector_p());
 …
   double max(const vector& v)
+  double max(const vectorBase& v)
+  {
     return gsl_vector_max(v.gsl_vector_p());
 …
   size_t max_index(const vector& v)
+  size_t max_index(const vectorBase& v)
+  {
     return gsl_vector_max_index(v.gsl_vector_p());
 …
   double min(const vector& v)
+  double min(const vectorBase& v)
+  {
     return gsl_vector_min(v.gsl_vector_p());
 …
   size_t min_index(const vector& v)
+  size_t min_index(const vectorBase& v)
+  {
     return gsl_vector_min_index(v.gsl_vector_p());
 …
+  bool nan(const vector& templat, vector& flag)
+  {
+    size_t vsize=templat.size();
+    if (vsize!=flag.size())
+      flag.clone(vector(vsize,1.0));
+    else
+      flag.all(1.0);
+  bool nan(const vectorBase& templat, vector& flag)
+  {
+    size_t vsize(templat.size());
+    flag = vector(vsize, 1.0);
     bool nan=false;
     for (size_t i=0; i<vsize; i++)
 …
+  void set_basis(vector& v, size_t i)
+  {
+    assert(v.gsl_vector_p());
+    gsl_vector_set_basis(v.gsl_vector_p(),i);
+  }
+  void shuffle(vector& invec)
+  {
+    random::DiscreteUniform rnd;
+    std::random_shuffle(invec.begin(), invec.end(), rnd);
+  }
+  void sort(vector& v)
+  {
+    assert(v.gsl_vector_p());
+    std::sort(v.begin(), v.end());
+  }
+  void sort_index(std::vector<size_t>& sort_index, const vector& invec)
+  void sort_index(std::vector<size_t>& sort_index, const vectorBase& invec)
+  {
     assert(invec.gsl_vector_p());
 …
     if (status) {
       gsl_permutation_free(p);
       throw utility::GSL_error(std::string("sort_index(vector&,const vector&)",status));
+      throw utility::GSL_error(std::string("sort_index(vector&,const vectorBase&)",status));
+    }
     sort_index=std::vector<size_t>(p->data,p->data+p->size);
 …
   void sort_smallest_index(std::vector<size_t>& sort_index, size_t k,
                             const vector& invec)
+                            const vectorBase& invec)
+  {
     assert(invec.gsl_vector_p());
 …
   void sort_largest_index(std::vector<size_t>& sort_index, size_t k,
                             const vector& invec)
+                            const vectorBase& invec)
+  {
     assert(invec.gsl_vector_p());
 …
   double sum(const vector& v)
+  double sum(const vectorBase& v)
+  {
     double sum = 0;
     size_t vsize=v.size();
     for (size_t i=0; i<vsize; ++i)
       sum += v[i];
+      sum += v(i);
     return sum;
+  }
+  void swap(vector& v, vector& w)
+  {
+    assert(v.gsl_vector_p()); assert(w.gsl_vector_p());
+    int status=gsl_vector_swap(v.gsl_vector_p(),w.gsl_vector_p());
+    if (status)
+      throw utility::GSL_error(std::string("swap(vector&,vector&)",status));
+  }
+  std::ostream& operator<<(std::ostream& s, const vector& a)
+  std::ostream& operator<<(std::ostream& s, const vectorBase& a)
+  {
     s.setf(std::ios::dec);
     s.precision(12);
     for (size_t j = 0; j < a.size(); ++j) {
       s << a[j];
+      s << a(j);
       if ( (j+1)<a.size() )
         s << s.fill();

branches/peter-dev/yat/utility/vectorBase.h

-                      r994
+                      r995
 #ifndef _theplu_yat_utility_vector_
 #define _theplu_yat_utility_vector_
+#ifndef _theplu_yat_utility_vector_base_
+#define _theplu_yat_utility_vector_base_
 // $Id$
 …
   class matrix;
+  class vector;
   /**
 …
   */
   class vector
+  class vectorBase
+  {
   public:
+    /// \brief vector::iterator
+    typedef Iterator<double&, vector> iterator;
+    /// \brief vector::const_iterator
+    typedef Iterator<const double, const vector> const_iterator;
+    /**
+       \brief The default constructor.
+    */
+    vector(void);
+    /**
+       \brief Allocates memory space for \a n elements, and sets all
+       elements to \a init_value.
+       \throw GSL_error if memory allocation fails.
+    */
+    vector(size_t n, double init_value=0);
+    /**
+       \brief The copy constructor.
+       \note If the object to be copied is a vector view, the values
+       of the view will be copied, i.e. the view is not copied.
+       \throw A GSL_error is indirectly thrown if memory allocation
+       fails.
+    */
+    vector(const vector& other);
+    /**
+       \brief Vector view constructor.
+       Create a view of vector \a v, with starting index \a offset,
+       size \a n, and an optional \a stride.
+       A vector view can be used as any vector with the difference
+       that changes made to the view will also change the object that
+       is viewed. Also, using the copy constructor will create a new
+       vector object that is a copy of whatever is viewed. If a copy
+       of the view is needed then you should use this constructor to
+       obtain a copy.
+       \note If the object viewed by the view goes out of scope or is
+       deleted, the view becomes invalid and the result of further use
+       is undefined.
+       \throw GSL_error if a view cannot be set up.
+    */
+    vector(vector& v, size_t offset, size_t n, size_t stride=1);
+    /**
+       \brief Vector const view constructor.
+       Create a view of vector \a v, with starting index \a offset,
+       size \a n, and an optional \a stride.
+       A vector view can be used as any const vector. Using the copy
+       constructor will create a new vector object that is a copy of
+       whatever is viewed. If a copy of the view is needed then you
+       should use this constructor to obtain a copy.
+       \note If the object viewed by the view goes out of scope or is
+       deleted, the view becomes invalid and the result of further use
+       is undefined.
+       \throw GSL_error if a view cannot be set up.
+    */
+    vector(const vector& v, size_t offset, size_t n, size_t stride=1);
+    ///
+    /// Matrix row/column view constructor.
+    ///
+    /// Create a row/column vector view of matrix \a m, pointing at
+    /// row/column \a i. The parameter \a row is used to set whether
+    /// the view should be a row or column view. If \a row is set to
+    /// true, the view will be a row view (default behaviour), and,
+    /// naturally, a column view otherwise.
+    ///
+    /// A vector view can be used as any vector with the difference
+    /// that changes made to the view will also change the object that
+    /// is viewed. Also, using the copy constructor will create a new
+    /// vector object that is a copy of whatever is viewed. If a copy
+    /// of the view is needed then you should use the vector view
+    /// constructor to obtain a copy.
+    ///
+    /// @note If the object viewed by the view goes out of scope or is
+    /// deleted, the view becomes invalid and the result of further
+    /// use is undefined.
+    ///
+    vector(matrix& m, size_t i, bool row=true);
+    ///
+    /// Matrix row/column const view constructor.
+    ///
+    /// Create a row/column vector view of matrix \a m, pointing at
+    /// row/column \a i. The parameter \a row is used to set whether
+    /// the view should be a row or column view. If \a row is set to
+    /// true, the view will be a row view (default behaviour), and,
+    /// naturally, a column view otherwise.
+    ///
+    /// A const vector view can be used as any const vector. Using the
+    /// copy constructor will create a new vector object that is a
+    /// copy of whatever is viewed. If a copy of the view is needed
+    /// then you should use the vector view constructor to obtain a
+    /// copy.
+    ///
+    /// @note If the object viewed by the view goes out of scope or is
+    /// deleted, the view becomes invalid and the result of further
+    /// use is undefined.
+    ///
+    vector(const matrix& m, size_t i, bool row=true);
+    /**
+       \brief The istream constructor.
+       Either elements should be separated with white space characters
+       (default), or elements should be separated by the delimiter \a
+       sep. When delimiter \a sep is used empty elements are stored as
+       NaN's (except that empty lines are ignored). The end of input
+       to the vector is at end of file marker.
+       \throw GSL_error if memory allocation fails, IO_error if
+       unexpected input is found in the input stream.
+    */
+    explicit vector(std::istream &, char sep='\0')
+      throw(utility::IO_error, std::exception);
+    /// \brief vectorBase::const_iterator
+    typedef Iterator<const double, const vectorBase> const_iterator;
+    /**
+       \brief Constructor.
+    */
+    vectorBase(const gsl_vector*);
     ///
     /// The destructor.
     ///
+    ~vector(void);
+    /**
+       \return mutable iterator to start of vector
+     */
+    iterator begin(void);
+    /**
+       \return read-only iterator to start of vector
+    virtual ~vectorBase(void);
+    /**
+       \return read-only iterator to start of vectorBase
      */
     const_iterator begin(void) const;
     /**
+       \brief Make a copy of \a other.
+       This function will make a deep copy of \a other. Memory is
+       resized and view state is changed if needed.
+    */
+    const vector& clone(const vector& other);
+    /**
+       \brief This function performs element-wise division, \f$ this_i =
+       this_i/other_i \; \forall i \f$.
+       \throw GSL_error if dimensions mis-match.
+    */
+    void div(const vector& other);
+    /**
+       \return mutable iterator to end of vector
+     */
+    iterator end(void);
+    /**
+       \return read-only iterator to end of vector
+       \return read-only iterator to end of vectorBase
      */
     const_iterator end(void) const;
     /**
        \brief Check whether vectors are equal within a user defined
+       \brief Check whether vectorBases are equal within a user defined
        precision, set by \a precision.
        \return True if each element deviates less or equal than \a
        d. If any vector contain a NaN, false is always returned.
+       d. If any vectorBase contain a NaN, false is always returned.
        \see operator== and operator!=
     */
     bool equal(const vector&, const double precision=0) const;
+    bool equal(const vectorBase&, const double precision=0) const;
     ///
 …
     ///
+    /// @return A pointer to the internal GSL vector,
+    ///
+    gsl_vector* gsl_vector_p(void);
+    ///
+    /// Check if the vector object is a view (sub-vector) to another
+    /// vector.
+    ///
+    /// @return True if the object is a view, false othwerwise.
+    ///
+    bool isview(void) const;
+    /**
+       \brief This function performs element-wise multiplication, \f$
+       this_i = this_i * other_i \; \forall i \f$.
+       \throw GSL_error if dimensions mis-match.
+    */
+    void mul(const vector& other);
+    /**
+       @brief Resize vector
+       All elements are set to @a init_value.
+       \note Underlying GSL vector is destroyed and a view into this
+       vector becomes invalid.
+    */
+    void resize(size_t, double init_value=0);
+    /**
+       \brief Reverse the order of elements in the vector.
+    */
+    void reverse(void);
+    ///
+    /// Set all elements to \a value.
+    ///
+    void all(const double& value);
+    ///
+    /// @return the number of elements in the vector.
+    /// @return the number of elements in the vectorBase.
     ///
     size_t size(void) const;
     /**
-       \brief Exchange elements \a i and \a j.
-       \throw GSL_error if vector lengths differs.
-    */
-    void swap(size_t i, size_t j);
-    /**
        \brief Element access operator.
        \return Reference to element \a i.
+       \return Const reference to element \a i.
        \throw If GSL range checks are enabled in the underlying GSL
 …
        of range.
     */
-    double& operator()(size_t i);
-    /**
-       \brief Element access operator.
-       \return Const reference to element \a i.
-       \throw If GSL range checks are enabled in the underlying GSL
-       library a GSL_error exception is thrown if either index is out
-       of range.
-    */
     const double& operator()(size_t i) const;
-    ///
-    /// Element access operator.
-    ///
-    /// @return Reference to element \a i.
-    ///
-    double& operator[](size_t i);
-    ///
-    /// Const element access operator.
-    ///
-    /// @return The value of element \a i.
-    ///
     const double& operator[](size_t i) const;
 …
        \return True if all elements are equal otherwise false.
        \see equal(const vector&, const double precision=0)
     */
     bool operator==(const vector&) const;
+       \see equal(const vectorBase&, const double precision=0)
+    */
+    bool operator==(const vectorBase&) const;
     /**
 …
        \return False if all elements are equal otherwise true.
        \see equal(const vector&, const double precision=0)
     */
     bool operator!=(const vector&) const;
+       \see equal(const vectorBase&, const double precision=0)
+    */
+    bool operator!=(const vectorBase&) const;
     ///
     /// @return The dot product.
     ///
+    double operator*(const vector&) const;
+    /**
+       \brief The assignment operator.
+       Dimensions of the vectors must match. If the LHS vector is a
+       view, the underlying data will be changed.
+       \return A const reference to the resulting vector.
+       \see void set(const vector&).
+       \throw GSL_error if dimensions mis-match.
+    */
+    const vector& operator=(const vector&);
+    /**
+       \brief Addition and assign operator. Vector addition, \f$
+       this_i = this_i + other_i \; \forall i \f$.
+       \return A const reference to the resulting vector.
+       \throw GSL_error if dimensions mis-match.
+    */
+    const vector& operator+=(const vector&);
+    /**
+       \brief Add a constant to a vector, \f$ this_i = this_i + d \;
+       \forall i \f$.
+       \return A const reference to the resulting vector.
+    */
+    const vector& operator+=(double d);
+    /**
+       \brief Subtract and assign operator. Vector subtraction, \f$
+       this_i = this_i - other_i \; \forall i \f$.
+       \return A const reference to the resulting vector.
+       \throw GSL_error if dimensions mis-match.
+    */
+    const vector& operator-=(const vector&);
+    /**
+       \brief Subtract a constant to a vector, \f$ this_i = this_i - d
+       \; \forall i \f$.
+       \return A const reference to the resulting vector.
+    */
+    const vector& operator-=(double d);
+    /**
+       \brief Multiply with scalar and assign operator, \f$ this_i =
+       this_i * d \; \forall i \f$.
+       \return A const reference to the resulting vector.
+    */
+    const vector& operator*=(double d);
+  private:
+    /**
+       \brief Create a new copy of the internal GSL vector.
+       Necessary memory for the new GSL vector is allocated and the
+       caller is responsible for freeing the allocated memory.
+       \return A pointer to a copy of the internal GSL vector.
+       \throw GSL_error if memory cannot be allocated for the new
+       copy, or if dimensions mis-match.
+    */
+    gsl_vector* create_gsl_vector_copy(void) const;
+    /**
+       \brief Clear all dynamically allocated memory.
+       Internal utility function.
+    */
+    void delete_allocated_memory(void);
+    gsl_vector* v_;
+    gsl_vector_view* view_;
+    const gsl_vector_const_view* view_const_;
+    // proxy_v_ is used to access the proper underlying gsl_vector
+    // in all const member functions. It is not used by design for
+    // non-const vector functions and operators. This is to make sure
+    // that runtime errors occur if a const vector is used in an
+    // inappropriate manner such as on left hand side in assignment
+    // (remember, v_ is null for const vector views).
+    const gsl_vector* proxy_v_;
+    double operator*(const vectorBase&) const;
+  protected:
+    const gsl_vector* const_vec_;
   };
   /**
      \brief Check if all elements of the vector are zero.
      \return True if all elements in the vector is zero, false
+     \brief Check if all elements of the vectorBase are zero.
+     \return True if all elements in the vectorBase is zero, false
      othwerwise.
   */
   bool isnull(const vector&);
   /**
      \brief Get the maximum value of the vector.
      \return The maximum value of the vector.
   */
   double max(const vector&);
   /**
      \brief Locate the maximum value in the vector.
      \return The index to the maximum value of the vector.
+  bool isnull(const vectorBase&);
+  /**
+     \brief Get the maximum value of the vectorBase.
+     \return The maximum value of the vectorBase.
+  */
+  double max(const vectorBase&);
+  /**
+     \brief Locate the maximum value in the vectorBase.
+     \return The index to the maximum value of the vectorBase.
      \note Lower index has precedence.
   */
   size_t max_index(const vector&);
   /**
      \brief Get the minimum value of the vector.
      \return The minimum value of the vector.
   */
   double min(const vector&);
   /**
      \brief Locate the minimum value in the vector.
      \return The index to the minimum value of the vector.
+  size_t max_index(const vectorBase&);
+  /**
+     \brief Get the minimum value of the vectorBase.
+     \return The minimum value of the vectorBase.
+  */
+  double min(const vectorBase&);
+  /**
+     \brief Locate the minimum value in the vectorBase.
+     \return The index to the minimum value of the vectorBase.
      \note Lower index has precedence.
   */
   size_t min_index(const vector&);
   /**
      \brief Create a vector \a flag indicating NaN's in another vector
+  size_t min_index(const vectorBase&);
+  /**
+     \brief Create a vectorBase \a flag indicating NaN's in another vectorBase
      \a templat.
      The \a flag vector is changed to contain 1's and 0's only. A 1
      means that the corresponding element in the \a templat vector is
+     The \a flag vectorBase is changed to contain 1's and 0's only. A 1
+     means that the corresponding element in the \a templat vectorBase is
      valid and a zero means that the corresponding element is a NaN.
+     \note Space for vector \a flag is reallocated to fit the size of
+     vector \a templat if sizes mismatch.
+     \return True if the \a templat vector contains at least one NaN.
+  */
+  bool nan(const vector& templat, vector& flag);
+  /**
+     \brief Transforms a vector to a basis vector.
+     All elements are set to zero except the \a i-th element which is
+     set to one.
+  */
+  void set_basis(vector&, size_t i);
+  /**
+     Randomly shuffles the elements in vector \a invec
+  */
+  void shuffle(vector& invec);
+  /**
+     Sort the elements in the vector.
+  */
+  void sort(vector&);
+  /**
+     Create a vector \a sort_index containing the indeces of
+     elements in a another vector \a invec.  The elements of \a
+     sort_index give the index of the vector element which would
+     have been stored in that position if the vector had been sorted
+     \note Space for vectorBase \a flag is reallocated to fit the size of
+     vectorBase \a templat if sizes mismatch.
+     \return True if the \a templat vectorBase contains at least one NaN.
+  */
+  bool nan(const vectorBase& templat, vector& flag);
+  /**
+     Create a vectorBase \a sort_index containing the indeces of
+     elements in a another vectorBase \a invec.  The elements of \a
+     sort_index give the index of the vectorBase element which would
+     have been stored in that position if the vectorBase had been sorted
      in place. The first element of \a sort_index gives the index of the least
+     element in \a invec, and the last element of \a sort_index gives the index of the
+     greatest element in \a invec . The vector \a invec is not changed.
+  */
+  void sort_index(std::vector<size_t>& sort_index, const vector& invec);
+  /** Similar to sort_index but creates a vector with indices to the \a k
+     element in \a invec, and the last element of \a sort_index gives the
+     index of the greatest element in \a invec . The vectorBase \a invec
+     is not changed.
+  */
+  void sort_index(std::vector<size_t>& sort_index, const vectorBase& invec);
+  /**
+      Similar to sort_index but creates a vectorBase with indices to the \a k
   smallest elements in \a invec.
   */
   void sort_smallest_index(std::vector<size_t>& sort_index, size_t k, const
+  vector& invec);
+  /** Similar to sort_index but creates a vector with indices to the \a k
+  vectorBase& invec);
+  /**
+      Similar to sort_index but creates a vectorBase with indices to the \a k
   largest elements in \a invec.
   */
   void sort_largest_index(std::vector<size_t>& sort_index, size_t k, const
   vector& invec);
+  vectorBase& invec);
   /**
      \brief Calculate the sum of all vector elements.
+     \brief Calculate the sum of all vectorBase elements.
      \return The sum.
   */
+  double sum(const vector&);
+  /**
+     \brief Swap vector elements by copying.
+     The two vectors must have the same length.
+     \throw GSL_error if vector lengths differs.
+  */
+  void swap(vector&, vector&);
+  /**
+     \brief The output operator for the vector class.
+  */
+  std::ostream& operator<<(std::ostream&, const vector&);
+  double sum(const vectorBase&);
+  /**
+     \brief The output operator for the vectorBase class.
+  */
+  std::ostream& operator<<(std::ostream&, const vectorBase&);
 }}} // of namespace utility, yat, and theplu

branches/peter-dev/yat/utility/vectorConstView.cc

-                      r994
+                      r995
 */
 #include "vector.h"
+#include "vectorConstView.h"
 #include "matrix.h"
 #include "utility.h"
 …
   vector::vector(void)
     : v_(NULL), view_(NULL), view_const_(NULL), proxy_v_(NULL)
+  vectorConstView::vectorConstView(void)
+    : vectorBase(NULL)
+  {
+  }
+  vector::vector(size_t n, double init_value)
+    : v_(gsl_vector_alloc(n)), view_(NULL), view_const_(NULL),
+      proxy_v_(v_)
+  {
+    if (!v_)
+      throw utility::GSL_error("vector::vector failed to allocate memory");
+    all(init_value);
+  }
+  vector::vector(const vector& other)
+    : v_(other.create_gsl_vector_copy()), view_(NULL),
+      view_const_(NULL), proxy_v_(v_)
+  vectorConstView::vectorConstView(const vectorBase& other)
+    : vectorBase(other.gsl_vector_p())
+  {
+  }
+  vector::vector(vector& v, size_t offset, size_t n, size_t stride)
+    : view_const_(NULL)
+  vectorConstView::vectorConstView(const vectorBase& v, size_t offset,
+                                   size_t n, size_t stride)
+    : vectorBase(NULL)
+  {
+    view_ = new gsl_vector_view(gsl_vector_subvector_with_stride(v.v_,offset,
+                                                                 stride,n));
+    if (!view_)
+      throw utility::GSL_error("vector::vector failed to setup view");
+    proxy_v_ = v_ = &(view_->vector);
+    const_view_ = new gsl_vector_const_view(
+                   gsl_vector_const_subvector_with_stride(v.gsl_vector_p(),
+                                                          offset, stride, n));
+    const_vec_ = &(const_view_->vector);
+  }
   vector::vector(const vector& v, size_t offset, size_t n, size_t stride)
     : v_(NULL), view_(NULL)
+  vectorConstView::vectorConstView(const matrix& m, size_t i, bool row)
+    : vectorBase(NULL)
+  {
     view_const_ = new gsl_vector_const_view(
                    gsl_vector_const_subvector_with_stride(v.v_,offset,stride,n));
     if (!view_const_)
       throw utility::GSL_error("vector::vector failed to setup view");
     proxy_v_ = &(view_const_->vector);
+    const_view_ =
+      new gsl_vector_const_view(row ?
+                                gsl_matrix_const_row(m.gsl_matrix_p(),i) :
+                                gsl_matrix_const_column(m.gsl_matrix_p(),i));
+    const_vec_ = &(const_view_->vector);
+  }
+  vector::vector(matrix& m, size_t i, bool row)
+    : view_const_(NULL)
+  vectorConstView::~vectorConstView(void)
+  {
-    view_=new gsl_vector_view(row ?
-                              gsl_matrix_row   (m.gsl_matrix_p(),i) :
-                              gsl_matrix_column(m.gsl_matrix_p(),i)  );
-    if (!view_)
-      throw utility::GSL_error("vector::vector failed to setup view");
-    proxy_v_ = v_ = &(view_->vector);
+  }
-  vector::vector(const matrix& m, size_t i, bool row)
-    : v_(NULL), view_(NULL)
+  {
-    view_const_ = new gsl_vector_const_view(row ?
-                                   gsl_matrix_const_row   (m.gsl_matrix_p(),i) :
-                                   gsl_matrix_const_column(m.gsl_matrix_p(),i) );
-    if (!view_const_)
-      throw utility::GSL_error("vector::vector failed to setup view");
-    proxy_v_ = &(view_const_->vector);
+  }
-  vector::vector(std::istream& is, char sep)
-    throw (utility::IO_error, std::exception)
-    : view_(NULL), view_const_(NULL)
+  {
-    // read the data file and store in stl vectors (dynamically
-    // expandable)
-    std::vector<std::vector<double> > data_matrix;
-    u_int nof_columns=0;
-    u_int nof_rows=0;
-    std::string line;
-    while(getline(is, line, '\n')) {
-      // Empty lines
-      if (!line.size())
-          continue;
-      nof_rows++;
-      std::vector<double> v;
-      std::string element;
-      std::stringstream ss(line);
-      bool ok=true;
-      while(ok) {
-        if(sep=='\0')
-          ok=(ss>>element);
-        else
-          ok=getline(ss, element, sep);
-        if(!ok)
-          break;
-        if(utility::is_double(element)) {
-          v.push_back(atof(element.c_str()));
+        }
-        else if (!element.size() || utility::is_nan(element)) {
-          v.push_back(std::numeric_limits<double>::quiet_NaN());
+        }
-        else {
-          std::stringstream ss("Warning: '");
-          ss << element << "' is not a double.";
-          throw IO_error(ss.str());
+        }
+      }
-      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
-          v.push_back(std::numeric_limits<double>::quiet_NaN());
-      if (!nof_columns)
-        nof_columns=v.size();
-      else if (nof_rows && (nof_columns>1)) {
-        std::ostringstream s;
-        s << "vector::vector(std::istream&) data file error:\n"
-          << "    File has inconsistent number of rows (" << nof_rows
-          << ") and columns (" << nof_columns
-          << ").\n    Expected a row or a column vector.";
-        throw utility::IO_error(s.str());
+      }
-      else if (v.size()!=nof_columns) {
-        std::ostringstream s;
-        s << "vector::vector(std::istream&) data file error:\n"
-          << "    Line " << nof_rows << " has " << v.size()
-          << " columns; expected " << nof_columns << " column.";
-        throw utility::IO_error(s.str());
+      }
-      data_matrix.push_back(v);
+    }
-    // manipulate the state of the stream to be good
-    is.clear(std::ios::goodbit);
-    // convert the data to a gsl vector
-    proxy_v_ = v_ = gsl_vector_alloc(nof_rows*nof_columns);
-    if (!v_)
-      throw utility::GSL_error("vector::vector failed to allocate memory");
-    size_t n=0;
-    // if gsl error handler disabled, out of bounds index will not
-    // abort the program.
-    for (size_t i=0; i<nof_rows; i++)
-      for (size_t j=0; j<nof_columns; j++)
-        gsl_vector_set( v_, n++, data_matrix[i][j] );
+  }
-  vector::~vector(void)
+  {
-    delete_allocated_memory();
+  }
-  vector::iterator vector::begin(void)
+  {
-    return vector::iterator(*this, 0);
+  }
-  vector::const_iterator vector::begin(void) const
+  {
-    return vector::const_iterator(*this, 0);
+  }
-  const vector& vector::clone(const vector& other)
+  {
-    if (this!=&other) {
-      delete_allocated_memory();
-      if (other.view_) {
-        view_ = new gsl_vector_view(*other.view_);
-        proxy_v_ = v_ = &(view_->vector);
+      }
-      else if (other.view_const_) {
-        view_const_ = new gsl_vector_const_view(*other.view_const_);
-        proxy_v_ = &(view_const_->vector);
-        v_=NULL;
+      }
-      else if (other.v_)
-        proxy_v_ = v_ = other.create_gsl_vector_copy();
+    }
-    return *this;
+  }
-  gsl_vector* vector::create_gsl_vector_copy(void) const
+  {
-    gsl_vector* vec = gsl_vector_alloc(size());
-    if (!vec)
-      throw utility::GSL_error("vector::create_gsl_vector_copy failed to allocate memory");
-    if (gsl_vector_memcpy(vec, proxy_v_))
-      throw utility::GSL_error("vector::create_gsl_matrix_copy dimension mis-match");
-    return vec;
+  }
-  void vector::delete_allocated_memory(void)
+  {
-    if (view_)
-      delete view_;
-    else if (view_const_)
-      delete view_const_;
-    else if (v_)
-      gsl_vector_free(v_);
-    proxy_v_=v_=NULL;
+  }
-  void vector::div(const vector& other)
+  {
-    assert(v_);
-    int status=gsl_vector_div(v_,other.gsl_vector_p());
-    if (status)
-      throw utility::GSL_error(std::string("vector::div",status));
+  }
-  vector::iterator vector::end(void)
+  {
-    return vector::iterator(*this, size());
+  }
-  vector::const_iterator vector::end(void) const
+  {
-    return vector::const_iterator(*this, size());
+  }
-  bool vector::equal(const vector& other, const double d) const
+  {
-    if (this==&other)
-      return true;
-    if (size()!=other.size())
-      return false;
-    // if gsl error handler disabled, out of bounds index will not
-    // abort the program.
-    for (size_t i=0; i<size(); ++i)
-      // The two last condition checks are needed for NaN detection
-      if (fabs( (*this)(i)-other(i) ) > d ||
-          (*this)(i)!=(*this)(i) || other(i)!=other(i))
-        return false;
-    return true;
+  }
-  const gsl_vector* vector::gsl_vector_p(void) const
+  {
-    return proxy_v_;
+  }
-  gsl_vector* vector::gsl_vector_p(void)
+  {
-    return v_;
+  }
-  bool vector::isview(void) const
+  {
-    return view_ || view_const_;
+  }
-  void vector::mul(const vector& other)
+  {
-    assert(v_);
-    int status=gsl_vector_mul(v_,other.gsl_vector_p());
-    if (status)
-      throw utility::GSL_error(std::string("vector::mul",status));
+  }
-  void vector::resize(size_t n, double init_value)
+  {
-    delete_allocated_memory();
-    proxy_v_ = v_ = gsl_vector_alloc(n);
-    if (!v_)
-      throw utility::GSL_error("vector::vector failed to allocate memory");
-    all(init_value);
+  }
-  void vector::reverse(void)
+  {
-    assert(v_);
-    gsl_vector_reverse(v_);
+  }
-  void vector::all(const double& value)
+  {
-    assert(v_);
-    gsl_vector_set_all(v_,value);
+  }
-  size_t vector::size(void) const
+  {
-    if (!proxy_v_)
-      return 0;
-    return proxy_v_->size;
+  }
-  void vector::swap(size_t i, size_t j)
+  {
-    assert(v_);
-    int status=gsl_vector_swap_elements(v_, i, j);
-    if (status)
-      throw utility::GSL_error(std::string("vector::swap_elements",status));
+  }
-  double& vector::operator()(size_t i)
+  {
-    assert(v_);
-    double* d=gsl_vector_ptr(v_, i);
-    if (!d)
-      throw utility::GSL_error("vector::operator()",GSL_EINVAL);
-    return *d;
+  }
-  const double& vector::operator()(size_t i) const
+  {
-    const double* d=gsl_vector_const_ptr(proxy_v_, i);
-    if (!d)
-      throw utility::GSL_error("vector::operator()",GSL_EINVAL);
-    return *d;
+  }
-  double& vector::operator[](size_t i)
+  {
-    return this->operator()(i);
+  }
-  const double& vector::operator[](size_t i) const
+  {
-    return this->operator()(i);
+  }
-  bool vector::operator==(const vector& other) const
+  {
-    return equal(other);
+  }
-  bool vector::operator!=(const vector& other) const
+  {
-    return !equal(other);
+  }
-  double vector::operator*( const vector &other ) const
+  {
-    double res = 0.0;;
-    for ( size_t i = 0; i < size(); ++i )
-      res += other(i) * (*this)(i);
-    return res;
+  }
-  const vector& vector::operator=( const vector& other )
+  {
-    assert(v_);
-    if (gsl_vector_memcpy(v_,other.gsl_vector_p()))
-      throw utility::GSL_error("vector::set dimension mis-match");
-    return *this;
+  }
-  const vector& vector::operator+=(const vector& other)
+  {
-    assert(v_);
-    int status=gsl_vector_add(v_, other.gsl_vector_p());
-    if (status)
-      throw utility::GSL_error(std::string("vector::add", status));
-    return *this;
+  }
-  const vector& vector::operator+=(double d)
+  {
-    assert(v_);
-    gsl_vector_add_constant(v_, d);
-    return *this;
+  }
-  const vector& vector::operator-=(const vector& other)
+  {
-    assert(v_);
-    int status=gsl_vector_sub(v_, other.gsl_vector_p());
-    if (status)
-      throw utility::GSL_error(std::string("vector::sub", status));
-    return *this;
+  }
-  const vector& vector::operator-=(const double d)
+  {
-    assert(v_);
-    gsl_vector_add_constant(v_, -d);
-    return *this;
+  }
-  const vector& vector::operator*=(const double d)
+  {
-    assert(v_);
-    gsl_vector_scale(v_, d);
-    return *this;
+  }
-  bool isnull(const vector& v)
+  {
-    return gsl_vector_isnull(v.gsl_vector_p());
+  }
-  double max(const vector& v)
+  {
-    return gsl_vector_max(v.gsl_vector_p());
+  }
-  size_t max_index(const vector& v)
+  {
-    return gsl_vector_max_index(v.gsl_vector_p());
+  }
-  double min(const vector& v)
+  {
-    return gsl_vector_min(v.gsl_vector_p());
+  }
-  size_t min_index(const vector& v)
+  {
-    return gsl_vector_min_index(v.gsl_vector_p());
+  }
-  bool nan(const vector& templat, vector& flag)
+  {
-    size_t vsize=templat.size();
-    if (vsize!=flag.size())
-      flag.clone(vector(vsize,1.0));
-    else
-      flag.all(1.0);
-    bool nan=false;
-    for (size_t i=0; i<vsize; i++)
-      if (std::isnan(templat(i))) {
-        flag(i)=0;
-        nan=true;
+      }
-    return nan;
+  }
-  void set_basis(vector& v, size_t i)
+  {
-    assert(v.gsl_vector_p());
-    gsl_vector_set_basis(v.gsl_vector_p(),i);
+  }
-  void shuffle(vector& invec)
+  {
-    random::DiscreteUniform rnd;
-    std::random_shuffle(invec.begin(), invec.end(), rnd);
+  }
-  void sort(vector& v)
+  {
-    assert(v.gsl_vector_p());
-    std::sort(v.begin(), v.end());
+  }
-  void sort_index(std::vector<size_t>& sort_index, const vector& invec)
+  {
-    assert(invec.gsl_vector_p());
-    gsl_permutation* p = gsl_permutation_alloc(invec.size());
-    int status=gsl_sort_vector_index(p,invec.gsl_vector_p());
-    if (status) {
-      gsl_permutation_free(p);
-      throw utility::GSL_error(std::string("sort_index(vector&,const vector&)",status));
+    }
-    sort_index=std::vector<size_t>(p->data,p->data+p->size);
-    gsl_permutation_free(p);
+  }
-  void sort_smallest_index(std::vector<size_t>& sort_index, size_t k,
-                            const vector& invec)
+  {
-    assert(invec.gsl_vector_p());
-    assert(k<=invec.size());
-    sort_index.resize(k);
-    gsl_sort_vector_smallest_index(&sort_index[0],k,invec.gsl_vector_p());
+  }
-  void sort_largest_index(std::vector<size_t>& sort_index, size_t k,
-                            const vector& invec)
+  {
-    assert(invec.gsl_vector_p());
-    assert(k<=invec.size());
-    sort_index.resize(k);
-    gsl_sort_vector_largest_index(&sort_index[0],k,invec.gsl_vector_p());
+  }
-  double sum(const vector& v)
+  {
-    double sum = 0;
-    size_t vsize=v.size();
-    for (size_t i=0; i<vsize; ++i)
-      sum += v[i];
-    return sum;
+  }
-  void swap(vector& v, vector& w)
+  {
-    assert(v.gsl_vector_p()); assert(w.gsl_vector_p());
-    int status=gsl_vector_swap(v.gsl_vector_p(),w.gsl_vector_p());
-    if (status)
-      throw utility::GSL_error(std::string("swap(vector&,vector&)",status));
+  }
-  std::ostream& operator<<(std::ostream& s, const vector& a)
+  {
-    s.setf(std::ios::dec);
-    s.precision(12);
-    for (size_t j = 0; j < a.size(); ++j) {
-      s << a[j];
-      if ( (j+1)<a.size() )
-        s << s.fill();
+    }
-    return s;
+  }
 }}} // of namespace utility, yat, and thep

branches/peter-dev/yat/utility/vectorConstView.h

-                      r994
+                      r995
 #ifndef _theplu_yat_utility_vector_
 #define _theplu_yat_utility_vector_
+#ifndef _theplu_yat_utility_vector_const_view_
+#define _theplu_yat_utility_vector_const_view_
 // $Id$
 …
 -1307, USA.
 */
+#include "vectorBase.h"
 #include "Exception.h"
 …
   */
   class vector
+  class vectorConstView : public vectorBase
+  {
   public:
-    /// \brief vector::iterator
-    typedef Iterator<double&, vector> iterator;
-    /// \brief vector::const_iterator
-    typedef Iterator<const double, const vector> const_iterator;
     /**
        \brief The default constructor.
+       \brief Default Constructor
     */
+    vector(void);
+    /**
+       \brief Allocates memory space for \a n elements, and sets all
+       elements to \a init_value.
+       \throw GSL_error if memory allocation fails.
+    */
+    vector(size_t n, double init_value=0);
+    vectorConstView(void);
     /**
        \brief The copy constructor.
        \note If the object to be copied is a vector view, the values
+       \note If the object to be copied is a vectorConstView view, the values
        of the view will be copied, i.e. the view is not copied.
 …
        fails.
     */
     vector(const vector& other);
+    vectorConstView(const vectorBase& other);
     /**
        \brief Vector view constructor.
+       \brief VectorConstView constructor.
        Create a view of vector \a v, with starting index \a offset,
+       Create a ConstView of vectorBase \a v, with starting index \a offset,
        size \a n, and an optional \a stride.
-       A vector view can be used as any vector with the difference
-       that changes made to the view will also change the object that
-       is viewed. Also, using the copy constructor will create a new
-       vector object that is a copy of whatever is viewed. If a copy
-       of the view is needed then you should use this constructor to
-       obtain a copy.
        \note If the object viewed by the view goes out of scope or is
 …
        \throw GSL_error if a view cannot be set up.
     */
+    vector(vector& v, size_t offset, size_t n, size_t stride=1);
+    /**
+       \brief Vector const view constructor.
+       Create a view of vector \a v, with starting index \a offset,
+       size \a n, and an optional \a stride.
+       A vector view can be used as any const vector. Using the copy
+       constructor will create a new vector object that is a copy of
+       whatever is viewed. If a copy of the view is needed then you
+       should use this constructor to obtain a copy.
+       \note If the object viewed by the view goes out of scope or is
+       deleted, the view becomes invalid and the result of further use
+       is undefined.
+       \throw GSL_error if a view cannot be set up.
+    */
+    vector(const vector& v, size_t offset, size_t n, size_t stride=1);
+    vectorConstView(const vectorBase& v, size_t offset, size_t n,
+                    size_t stride=1);
     ///
     /// Matrix row/column view constructor.
+    /// Matrix row/column const view constructor.
     ///
     /// Create a row/column vector view of matrix \a m, pointing at
+    /// Create a row/column vectorConstView view of matrix \a m, pointing at
     /// row/column \a i. The parameter \a row is used to set whether
     /// the view should be a row or column view. If \a row is set to
 …
     /// naturally, a column view otherwise.
     ///
-    /// A vector view can be used as any vector with the difference
-    /// that changes made to the view will also change the object that
-    /// is viewed. Also, using the copy constructor will create a new
-    /// vector object that is a copy of whatever is viewed. If a copy
-    /// of the view is needed then you should use the vector view
-    /// constructor to obtain a copy.
-    ///
     /// @note If the object viewed by the view goes out of scope or is
     /// deleted, the view becomes invalid and the result of further
     /// use is undefined.
     ///
+    vector(matrix& m, size_t i, bool row=true);
+    ///
+    /// Matrix row/column const view constructor.
+    ///
+    /// Create a row/column vector view of matrix \a m, pointing at
+    /// row/column \a i. The parameter \a row is used to set whether
+    /// the view should be a row or column view. If \a row is set to
+    /// true, the view will be a row view (default behaviour), and,
+    /// naturally, a column view otherwise.
+    ///
+    /// A const vector view can be used as any const vector. Using the
+    /// copy constructor will create a new vector object that is a
+    /// copy of whatever is viewed. If a copy of the view is needed
+    /// then you should use the vector view constructor to obtain a
+    /// copy.
+    ///
+    /// @note If the object viewed by the view goes out of scope or is
+    /// deleted, the view becomes invalid and the result of further
+    /// use is undefined.
+    ///
+    vector(const matrix& m, size_t i, bool row=true);
+    /**
+       \brief The istream constructor.
+       Either elements should be separated with white space characters
+       (default), or elements should be separated by the delimiter \a
+       sep. When delimiter \a sep is used empty elements are stored as
+       NaN's (except that empty lines are ignored). The end of input
+       to the vector is at end of file marker.
+       \throw GSL_error if memory allocation fails, IO_error if
+       unexpected input is found in the input stream.
+    */
+    explicit vector(std::istream &, char sep='\0')
+      throw(utility::IO_error, std::exception);
+    vectorConstView(const matrix& m, size_t i, bool row=true);
     ///
     /// The destructor.
     ///
+    ~vector(void);
+    /**
+       \return mutable iterator to start of vector
+     */
+    iterator begin(void);
+    /**
+       \return read-only iterator to start of vector
+     */
+    const_iterator begin(void) const;
+    /**
+       \brief Make a copy of \a other.
+       This function will make a deep copy of \a other. Memory is
+       resized and view state is changed if needed.
+    */
+    const vector& clone(const vector& other);
+    /**
+       \brief This function performs element-wise division, \f$ this_i =
+       this_i/other_i \; \forall i \f$.
+       \throw GSL_error if dimensions mis-match.
+    */
+    void div(const vector& other);
+    /**
+       \return mutable iterator to end of vector
+     */
+    iterator end(void);
+    /**
+       \return read-only iterator to end of vector
+     */
+    const_iterator end(void) const;
+    /**
+       \brief Check whether vectors are equal within a user defined
+       precision, set by \a precision.
+       \return True if each element deviates less or equal than \a
+       d. If any vector contain a NaN, false is always returned.
+       \see operator== and operator!=
+    */
+    bool equal(const vector&, const double precision=0) const;
+    ///
+    /// @return A const pointer to the internal GSL vector,
+    ///
+    const gsl_vector* gsl_vector_p(void) const;
+    ///
+    /// @return A pointer to the internal GSL vector,
+    ///
+    gsl_vector* gsl_vector_p(void);
+    ///
+    /// Check if the vector object is a view (sub-vector) to another
+    /// vector.
+    ///
+    /// @return True if the object is a view, false othwerwise.
+    ///
+    bool isview(void) const;
+    /**
+       \brief This function performs element-wise multiplication, \f$
+       this_i = this_i * other_i \; \forall i \f$.
+       \throw GSL_error if dimensions mis-match.
+    */
+    void mul(const vector& other);
+    /**
+       @brief Resize vector
+       All elements are set to @a init_value.
+       \note Underlying GSL vector is destroyed and a view into this
+       vector becomes invalid.
+    */
+    void resize(size_t, double init_value=0);
+    /**
+       \brief Reverse the order of elements in the vector.
+    */
+    void reverse(void);
+    ///
+    /// Set all elements to \a value.
+    ///
+    void all(const double& value);
+    ///
+    /// @return the number of elements in the vector.
+    ///
+    size_t size(void) const;
+    /**
+       \brief Exchange elements \a i and \a j.
+       \throw GSL_error if vector lengths differs.
+    */
+    void swap(size_t i, size_t j);
+    /**
+       \brief Element access operator.
+       \return Reference to element \a i.
+       \throw If GSL range checks are enabled in the underlying GSL
+       library a GSL_error exception is thrown if either index is out
+       of range.
+    */
+    double& operator()(size_t i);
+    /**
+       \brief Element access operator.
+       \return Const reference to element \a i.
+       \throw If GSL range checks are enabled in the underlying GSL
+       library a GSL_error exception is thrown if either index is out
+       of range.
+    */
+    const double& operator()(size_t i) const;
+    ///
+    /// Element access operator.
+    ///
+    /// @return Reference to element \a i.
+    ///
+    double& operator[](size_t i);
+    ///
+    /// Const element access operator.
+    ///
+    /// @return The value of element \a i.
+    ///
+    const double& operator[](size_t i) const;
+    /**
+       \brief Comparison operator. Takes linear time.
+       Checks are performed with exact matching, i.e., rounding off
+       effects may destroy comparison. Use the equal function for
+       comparing elements within a user defined precision.
+       \return True if all elements are equal otherwise false.
+       \see equal(const vector&, const double precision=0)
+    */
+    bool operator==(const vector&) const;
+    /**
+       \brief Comparison operator. Takes linear time.
+       Checks are performed with exact matching, i.e., rounding off
+       effects may destroy comparison. Use the equal function for
+       comparing elements within a user defined precision.
+       \return False if all elements are equal otherwise true.
+       \see equal(const vector&, const double precision=0)
+    */
+    bool operator!=(const vector&) const;
+    ///
+    /// @return The dot product.
+    ///
+    double operator*(const vector&) const;
+    /**
+       \brief The assignment operator.
+       Dimensions of the vectors must match. If the LHS vector is a
+       view, the underlying data will be changed.
+       \return A const reference to the resulting vector.
+       \see void set(const vector&).
+       \throw GSL_error if dimensions mis-match.
+    */
+    const vector& operator=(const vector&);
+    /**
+       \brief Addition and assign operator. Vector addition, \f$
+       this_i = this_i + other_i \; \forall i \f$.
+       \return A const reference to the resulting vector.
+       \throw GSL_error if dimensions mis-match.
+    */
+    const vector& operator+=(const vector&);
+    /**
+       \brief Add a constant to a vector, \f$ this_i = this_i + d \;
+       \forall i \f$.
+       \return A const reference to the resulting vector.
+    */
+    const vector& operator+=(double d);
+    /**
+       \brief Subtract and assign operator. Vector subtraction, \f$
+       this_i = this_i - other_i \; \forall i \f$.
+       \return A const reference to the resulting vector.
+       \throw GSL_error if dimensions mis-match.
+    */
+    const vector& operator-=(const vector&);
+    /**
+       \brief Subtract a constant to a vector, \f$ this_i = this_i - d
+       \; \forall i \f$.
+       \return A const reference to the resulting vector.
+    */
+    const vector& operator-=(double d);
+    /**
+       \brief Multiply with scalar and assign operator, \f$ this_i =
+       this_i * d \; \forall i \f$.
+       \return A const reference to the resulting vector.
+    */
+    const vector& operator*=(double d);
+    ~vectorConstView(void);
   private:
+    /**
+       \brief Create a new copy of the internal GSL vector.
+       Necessary memory for the new GSL vector is allocated and the
+       caller is responsible for freeing the allocated memory.
+       \return A pointer to a copy of the internal GSL vector.
+       \throw GSL_error if memory cannot be allocated for the new
+       copy, or if dimensions mis-match.
+    */
+    gsl_vector* create_gsl_vector_copy(void) const;
+    /**
+       \brief Clear all dynamically allocated memory.
+       Internal utility function.
+    */
+    void delete_allocated_memory(void);
+    gsl_vector* v_;
+    gsl_vector_view* view_;
+    const gsl_vector_const_view* view_const_;
+    // proxy_v_ is used to access the proper underlying gsl_vector
+    // in all const member functions. It is not used by design for
+    // non-const vector functions and operators. This is to make sure
+    // that runtime errors occur if a const vector is used in an
+    // inappropriate manner such as on left hand side in assignment
+    // (remember, v_ is null for const vector views).
+    const gsl_vector* proxy_v_;
+    const gsl_vector_const_view* const_view_;
   };
+  /**
+     \brief Check if all elements of the vector are zero.
+     \return True if all elements in the vector is zero, false
+     othwerwise.
+  */
+  bool isnull(const vector&);
+  /**
+     \brief Get the maximum value of the vector.
+     \return The maximum value of the vector.
+  */
+  double max(const vector&);
+  /**
+     \brief Locate the maximum value in the vector.
+     \return The index to the maximum value of the vector.
+     \note Lower index has precedence.
+  */
+  size_t max_index(const vector&);
+  /**
+     \brief Get the minimum value of the vector.
+     \return The minimum value of the vector.
+  */
+  double min(const vector&);
+  /**
+     \brief Locate the minimum value in the vector.
+     \return The index to the minimum value of the vector.
+     \note Lower index has precedence.
+  */
+  size_t min_index(const vector&);
+  /**
+     \brief Create a vector \a flag indicating NaN's in another vector
+     \a templat.
+     The \a flag vector is changed to contain 1's and 0's only. A 1
+     means that the corresponding element in the \a templat vector is
+     valid and a zero means that the corresponding element is a NaN.
+     \note Space for vector \a flag is reallocated to fit the size of
+     vector \a templat if sizes mismatch.
+     \return True if the \a templat vector contains at least one NaN.
+  */
+  bool nan(const vector& templat, vector& flag);
+  /**
+     \brief Transforms a vector to a basis vector.
+     All elements are set to zero except the \a i-th element which is
+     set to one.
+  */
+  void set_basis(vector&, size_t i);
+  /**
+     Randomly shuffles the elements in vector \a invec
+  */
+  void shuffle(vector& invec);
+  /**
+     Sort the elements in the vector.
+  */
+  void sort(vector&);
+  /**
+     Create a vector \a sort_index containing the indeces of
+     elements in a another vector \a invec.  The elements of \a
+     sort_index give the index of the vector element which would
+     have been stored in that position if the vector had been sorted
+     in place. The first element of \a sort_index gives the index of the least
+     element in \a invec, and the last element of \a sort_index gives the index of the
+     greatest element in \a invec . The vector \a invec is not changed.
+  */
+  void sort_index(std::vector<size_t>& sort_index, const vector& invec);
+  /** Similar to sort_index but creates a vector with indices to the \a k
+  smallest elements in \a invec.
+  */
+  void sort_smallest_index(std::vector<size_t>& sort_index, size_t k, const
+  vector& invec);
+  /** Similar to sort_index but creates a vector with indices to the \a k
+  largest elements in \a invec.
+  */
+  void sort_largest_index(std::vector<size_t>& sort_index, size_t k, const
+  vector& invec);
+  /**
+     \brief Calculate the sum of all vector elements.
+     \return The sum.
+  */
+  double sum(const vector&);
+  /**
+     \brief Swap vector elements by copying.
+     The two vectors must have the same length.
+     \throw GSL_error if vector lengths differs.
+  */
+  void swap(vector&, vector&);
+  /**
+     \brief The output operator for the vector class.
+  */
+  std::ostream& operator<<(std::ostream&, const vector&);
+}}} // of namespace utility, yat, and theplu
+}}}  // of namespace utility, yat, and theplu
 #endif

branches/peter-dev/yat/utility/vectorMutable.cc

-                      r994
+                      r995
 */
 #include "vector.h"
+#include "vectorMutable.h"
 #include "matrix.h"
 #include "utility.h"
 …
   vector::vector(void)
     : v_(NULL), view_(NULL), view_const_(NULL), proxy_v_(NULL)
+  vectorMutable::vectorMutable(gsl_vector* v)
+    : vectorBase(v)
+  {
+  }
+  vector::vector(size_t n, double init_value)
+    : v_(gsl_vector_alloc(n)), view_(NULL), view_const_(NULL),
+      proxy_v_(v_)
+  {
+    if (!v_)
+      throw utility::GSL_error("vector::vector failed to allocate memory");
+    all(init_value);
+  }
+  vector::vector(const vector& other)
+    : v_(other.create_gsl_vector_copy()), view_(NULL),
+      view_const_(NULL), proxy_v_(v_)
+  vectorMutable::~vectorMutable(void)
+  {
+  }
+  vector::vector(vector& v, size_t offset, size_t n, size_t stride)
+    : view_const_(NULL)
+  vectorMutable::iterator vectorMutable::begin(void)
+  {
+    view_ = new gsl_vector_view(gsl_vector_subvector_with_stride(v.v_,offset,
+                                                                 stride,n));
+    if (!view_)
+      throw utility::GSL_error("vector::vector failed to setup view");
+    proxy_v_ = v_ = &(view_->vector);
+    return vectorMutable::iterator(*this, 0);
+  }
+  vector::vector(const vector& v, size_t offset, size_t n, size_t stride)
+    : v_(NULL), view_(NULL)
+  void vectorMutable::div(const vectorBase& other)
+  {
+    view_const_ = new gsl_vector_const_view(
+                   gsl_vector_const_subvector_with_stride(v.v_,offset,stride,n));
+    if (!view_const_)
+      throw utility::GSL_error("vector::vector failed to setup view");
+    proxy_v_ = &(view_const_->vector);
+    assert(vec_);
+    int status=gsl_vector_div(vec_,other.gsl_vector_p());
+    if (status)
+      throw utility::GSL_error(std::string("vectorMutable::div",status));
+  }
+  vector::vector(matrix& m, size_t i, bool row)
+    : view_const_(NULL)
+  vectorMutable::iterator vectorMutable::end(void)
+  {
+    view_=new gsl_vector_view(row ?
+                              gsl_matrix_row   (m.gsl_matrix_p(),i) :
+                              gsl_matrix_column(m.gsl_matrix_p(),i)  );
+    if (!view_)
+      throw utility::GSL_error("vector::vector failed to setup view");
+    proxy_v_ = v_ = &(view_->vector);
+    return vectorMutable::iterator(*this, size());
+  }
+  vector::vector(const matrix& m, size_t i, bool row)
+    : v_(NULL), view_(NULL)
+  gsl_vector* vectorMutable::gsl_vector_p(void)
+  {
+    view_const_ = new gsl_vector_const_view(row ?
+                                   gsl_matrix_const_row   (m.gsl_matrix_p(),i) :
+                                   gsl_matrix_const_column(m.gsl_matrix_p(),i) );
+    if (!view_const_)
+      throw utility::GSL_error("vector::vector failed to setup view");
+    proxy_v_ = &(view_const_->vector);
+    return vec_;
+  }
+  vector::vector(std::istream& is, char sep)
+    throw (utility::IO_error, std::exception)
+    : view_(NULL), view_const_(NULL)
+  void vectorMutable::mul(const vectorBase& other)
+  {
+    // read the data file and store in stl vectors (dynamically
+    // expandable)
+    std::vector<std::vector<double> > data_matrix;
+    u_int nof_columns=0;
+    u_int nof_rows=0;
+    std::string line;
+    while(getline(is, line, '\n')) {
+      // Empty lines
+      if (!line.size())
+          continue;
+      nof_rows++;
+      std::vector<double> v;
+      std::string element;
+      std::stringstream ss(line);
+      bool ok=true;
+      while(ok) {
+        if(sep=='\0')
+          ok=(ss>>element);
+        else
+          ok=getline(ss, element, sep);
+        if(!ok)
+          break;
+        if(utility::is_double(element)) {
+          v.push_back(atof(element.c_str()));
+        }
+        else if (!element.size() || utility::is_nan(element)) {
+          v.push_back(std::numeric_limits<double>::quiet_NaN());
+        }
+        else {
+          std::stringstream ss("Warning: '");
+          ss << element << "' is not a double.";
+          throw IO_error(ss.str());
+        }
+      }
+      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
+          v.push_back(std::numeric_limits<double>::quiet_NaN());
+      if (!nof_columns)
+        nof_columns=v.size();
+      else if (nof_rows && (nof_columns>1)) {
+        std::ostringstream s;
+        s << "vector::vector(std::istream&) data file error:\n"
+          << "    File has inconsistent number of rows (" << nof_rows
+          << ") and columns (" << nof_columns
+          << ").\n    Expected a row or a column vector.";
+        throw utility::IO_error(s.str());
+      }
+      else if (v.size()!=nof_columns) {
+        std::ostringstream s;
+        s << "vector::vector(std::istream&) data file error:\n"
+          << "    Line " << nof_rows << " has " << v.size()
+          << " columns; expected " << nof_columns << " column.";
+        throw utility::IO_error(s.str());
+      }
+      data_matrix.push_back(v);
+    }
+    // manipulate the state of the stream to be good
+    is.clear(std::ios::goodbit);
+    // convert the data to a gsl vector
+    proxy_v_ = v_ = gsl_vector_alloc(nof_rows*nof_columns);
+    if (!v_)
+      throw utility::GSL_error("vector::vector failed to allocate memory");
+    size_t n=0;
+    // if gsl error handler disabled, out of bounds index will not
+    // abort the program.
+    for (size_t i=0; i<nof_rows; i++)
+      for (size_t j=0; j<nof_columns; j++)
+        gsl_vector_set( v_, n++, data_matrix[i][j] );
+    assert(vec_);
+    int status=gsl_vector_mul(vec_,other.gsl_vector_p());
+    if (status)
+      throw utility::GSL_error(std::string("vectorMutable::mul",status));
+  }
   vector::~vector(void)
+  void vectorMutable::reverse(void)
+  {
+    delete_allocated_memory();
+    assert(vec_);
+    gsl_vector_reverse(vec_);
+  }
   vector::iterator vector::begin(void)
+  void vectorMutable::all(const double& value)
+  {
+    return vector::iterator(*this, 0);
+    assert(vec_);
+    gsl_vector_set_all(vec_,value);
+  }
   vector::const_iterator vector::begin(void) const
+  void vectorMutable::swap(size_t i, size_t j)
+  {
+    return vector::const_iterator(*this, 0);
+    assert(vec_);
+    int status=gsl_vector_swap_elements(vec_, i, j);
+    if (status)
+      throw utility::GSL_error(std::string("vectorMutable::swap_elements",status));
+  }
   const vector& vector::clone(const vector& other)
+  double& vectorMutable::operator()(size_t i)
+  {
+    if (this!=&other) {
+      delete_allocated_memory();
+      if (other.view_) {
+        view_ = new gsl_vector_view(*other.view_);
+        proxy_v_ = v_ = &(view_->vector);
+      }
+      else if (other.view_const_) {
+        view_const_ = new gsl_vector_const_view(*other.view_const_);
+        proxy_v_ = &(view_const_->vector);
+        v_=NULL;
+      }
+      else if (other.v_)
+        proxy_v_ = v_ = other.create_gsl_vector_copy();
+    }
+    return *this;
+  }
+  gsl_vector* vector::create_gsl_vector_copy(void) const
+  {
+    gsl_vector* vec = gsl_vector_alloc(size());
+    if (!vec)
+      throw utility::GSL_error("vector::create_gsl_vector_copy failed to allocate memory");
+    if (gsl_vector_memcpy(vec, proxy_v_))
+      throw utility::GSL_error("vector::create_gsl_matrix_copy dimension mis-match");
+    return vec;
+  }
+  void vector::delete_allocated_memory(void)
+  {
+    if (view_)
+      delete view_;
+    else if (view_const_)
+      delete view_const_;
+    else if (v_)
+      gsl_vector_free(v_);
+    proxy_v_=v_=NULL;
+  }
+  void vector::div(const vector& other)
+  {
+    assert(v_);
+    int status=gsl_vector_div(v_,other.gsl_vector_p());
+    if (status)
+      throw utility::GSL_error(std::string("vector::div",status));
+  }
+  vector::iterator vector::end(void)
+  {
+    return vector::iterator(*this, size());
+  }
+  vector::const_iterator vector::end(void) const
+  {
+    return vector::const_iterator(*this, size());
+  }
+  bool vector::equal(const vector& other, const double d) const
+  {
+    if (this==&other)
+      return true;
+    if (size()!=other.size())
+      return false;
+    // if gsl error handler disabled, out of bounds index will not
+    // abort the program.
+    for (size_t i=0; i<size(); ++i)
+      // The two last condition checks are needed for NaN detection
+      if (fabs( (*this)(i)-other(i) ) > d ||
+          (*this)(i)!=(*this)(i) || other(i)!=other(i))
+        return false;
+    return true;
+  }
+  const gsl_vector* vector::gsl_vector_p(void) const
+  {
+    return proxy_v_;
+  }
+  gsl_vector* vector::gsl_vector_p(void)
+  {
+    return v_;
+  }
+  bool vector::isview(void) const
+  {
+    return view_ || view_const_;
+  }
+  void vector::mul(const vector& other)
+  {
+    assert(v_);
+    int status=gsl_vector_mul(v_,other.gsl_vector_p());
+    if (status)
+      throw utility::GSL_error(std::string("vector::mul",status));
+  }
+  void vector::resize(size_t n, double init_value)
+  {
+    delete_allocated_memory();
+    proxy_v_ = v_ = gsl_vector_alloc(n);
+    if (!v_)
+      throw utility::GSL_error("vector::vector failed to allocate memory");
+    all(init_value);
+  }
+  void vector::reverse(void)
+  {
+    assert(v_);
+    gsl_vector_reverse(v_);
+  }
+  void vector::all(const double& value)
+  {
+    assert(v_);
+    gsl_vector_set_all(v_,value);
+  }
+  size_t vector::size(void) const
+  {
+    if (!proxy_v_)
+      return 0;
+    return proxy_v_->size;
+  }
+  void vector::swap(size_t i, size_t j)
+  {
+    assert(v_);
+    int status=gsl_vector_swap_elements(v_, i, j);
+    if (status)
+      throw utility::GSL_error(std::string("vector::swap_elements",status));
+  }
+  double& vector::operator()(size_t i)
+  {
+    assert(v_);
+    double* d=gsl_vector_ptr(v_, i);
+    assert(vec_);
+    double* d=gsl_vector_ptr(vec_, i);
     if (!d)
       throw utility::GSL_error("vector::operator()",GSL_EINVAL);
+      throw utility::GSL_error("vectorMutable::operator()",GSL_EINVAL);
     return *d;
+  }
   const double& vector::operator()(size_t i) const
+  const vectorMutable& vectorMutable::operator+=(const vectorBase& other)
+  {
+    const double* d=gsl_vector_const_ptr(proxy_v_, i);
+    if (!d)
+      throw utility::GSL_error("vector::operator()",GSL_EINVAL);
+    return *d;
+  }
+  double& vector::operator[](size_t i)
+  {
+    return this->operator()(i);
+  }
+  const double& vector::operator[](size_t i) const
+  {
+    return this->operator()(i);
+  }
+  bool vector::operator==(const vector& other) const
+  {
+    return equal(other);
+  }
+  bool vector::operator!=(const vector& other) const
+  {
+    return !equal(other);
+  }
+  double vector::operator*( const vector &other ) const
+  {
+    double res = 0.0;;
+    for ( size_t i = 0; i < size(); ++i )
+      res += other(i) * (*this)(i);
+    return res;
+  }
+  const vector& vector::operator=( const vector& other )
+  {
+    assert(v_);
+    if (gsl_vector_memcpy(v_,other.gsl_vector_p()))
+      throw utility::GSL_error("vector::set dimension mis-match");
+    return *this;
+  }
+  const vector& vector::operator+=(const vector& other)
+  {
+    assert(v_);
+    int status=gsl_vector_add(v_, other.gsl_vector_p());
+    assert(vec_);
+    int status=gsl_vector_add(vec_, other.gsl_vector_p());
     if (status)
       throw utility::GSL_error(std::string("vector::add", status));
+      throw utility::GSL_error(std::string("vectorMutable::add", status));
     return *this;
+  }
   const vector& vector::operator+=(double d)
+  const vectorMutable& vectorMutable::operator+=(double d)
+  {
     assert(v_);
     gsl_vector_add_constant(v_, d);
+    assert(vec_);
+    gsl_vector_add_constant(vec_, d);
     return *this;
+  }
   const vector& vector::operator-=(const vector& other)
+  const vectorMutable& vectorMutable::operator-=(const vectorBase& other)
+  {
     assert(v_);
     int status=gsl_vector_sub(v_, other.gsl_vector_p());
+    assert(vec_);
+    int status=gsl_vector_sub(vec_, other.gsl_vector_p());
     if (status)
       throw utility::GSL_error(std::string("vector::sub", status));
+      throw utility::GSL_error(std::string("vectorMutable::sub", status));
     return *this;
+  }
   const vector& vector::operator-=(const double d)
+  const vectorMutable& vectorMutable::operator-=(const double d)
+  {
     assert(v_);
     gsl_vector_add_constant(v_, -d);
+    assert(vec_);
+    gsl_vector_add_constant(vec_, -d);
     return *this;
+  }
   const vector& vector::operator*=(const double d)
+  const vectorMutable& vectorMutable::operator*=(const double d)
+  {
     assert(v_);
     gsl_vector_scale(v_, d);
+    assert(vec_);
+    gsl_vector_scale(vec_, d);
     return *this;
+  }
+  bool isnull(const vector& v)
+  {
+    return gsl_vector_isnull(v.gsl_vector_p());
+  }
+  double max(const vector& v)
+  {
+    return gsl_vector_max(v.gsl_vector_p());
+  }
+  size_t max_index(const vector& v)
+  {
+    return gsl_vector_max_index(v.gsl_vector_p());
+  }
+  double min(const vector& v)
+  {
+    return gsl_vector_min(v.gsl_vector_p());
+  }
+  size_t min_index(const vector& v)
+  {
+    return gsl_vector_min_index(v.gsl_vector_p());
+  }
+  bool nan(const vector& templat, vector& flag)
+  {
+    size_t vsize=templat.size();
+    if (vsize!=flag.size())
+      flag.clone(vector(vsize,1.0));
+    else
+      flag.all(1.0);
+    bool nan=false;
+    for (size_t i=0; i<vsize; i++)
+      if (std::isnan(templat(i))) {
+        flag(i)=0;
+        nan=true;
+      }
+    return nan;
+  }
+  void set_basis(vector& v, size_t i)
+  void set_basis(vectorMutable& v, size_t i)
+  {
     assert(v.gsl_vector_p());
 …
   void shuffle(vector& invec)
+  void shuffle(vectorMutable& invec)
+  {
     random::DiscreteUniform rnd;
 …
   void sort(vector& v)
+  void sort(vectorMutable& v)
+  {
     assert(v.gsl_vector_p());
 …
+  }
+  void sort_index(std::vector<size_t>& sort_index, const vector& invec)
+  vectorMutable::operator proxy ()
+  {
+    assert(invec.gsl_vector_p());
+    gsl_permutation* p = gsl_permutation_alloc(invec.size());
+    int status=gsl_sort_vector_index(p,invec.gsl_vector_p());
+    if (status) {
+      gsl_permutation_free(p);
+      throw utility::GSL_error(std::string("sort_index(vector&,const vector&)",status));
+    }
+    sort_index=std::vector<size_t>(p->data,p->data+p->size);
+    gsl_permutation_free(p);
+    proxy p;
+    p.vec_ = this->vec_;
+    return p;
+  }
-  void sort_smallest_index(std::vector<size_t>& sort_index, size_t k,
-                            const vector& invec)
+  {
-    assert(invec.gsl_vector_p());
-    assert(k<=invec.size());
-    sort_index.resize(k);
-    gsl_sort_vector_smallest_index(&sort_index[0],k,invec.gsl_vector_p());
+  }
-  void sort_largest_index(std::vector<size_t>& sort_index, size_t k,
-                            const vector& invec)
+  {
-    assert(invec.gsl_vector_p());
-    assert(k<=invec.size());
-    sort_index.resize(k);
-    gsl_sort_vector_largest_index(&sort_index[0],k,invec.gsl_vector_p());
+  }
-  double sum(const vector& v)
+  {
-    double sum = 0;
-    size_t vsize=v.size();
-    for (size_t i=0; i<vsize; ++i)
-      sum += v[i];
-    return sum;
+  }
-  void swap(vector& v, vector& w)
+  {
-    assert(v.gsl_vector_p()); assert(w.gsl_vector_p());
-    int status=gsl_vector_swap(v.gsl_vector_p(),w.gsl_vector_p());
-    if (status)
-      throw utility::GSL_error(std::string("swap(vector&,vector&)",status));
+  }
-  std::ostream& operator<<(std::ostream& s, const vector& a)
+  {
-    s.setf(std::ios::dec);
-    s.precision(12);
-    for (size_t j = 0; j < a.size(); ++j) {
-      s << a[j];
-      if ( (j+1)<a.size() )
-        s << s.fill();
+    }
-    return s;
+  }
 }}} // of namespace utility, yat, and thep

branches/peter-dev/yat/utility/vectorMutable.h

-                      r994
+                      r995
 #ifndef _theplu_yat_utility_vector_
 #define _theplu_yat_utility_vector_
+#ifndef _theplu_yat_utility_vector_mutable_
+#define _theplu_yat_utility_vector_mutable_
 // $Id$
 …
 */
+#include "vectorBase.h"
 #include "Exception.h"
 #include "Iterator.h"
 …
   class matrix;
   /**
      @brief This is the yat interface to GSL vector.
 …
   */
   class vector
+  class vectorMutable : public vectorBase
+  {
   public:
+    /// \brief vector::iterator
+    typedef Iterator<double&, vector> iterator;
+    /// \brief vector::const_iterator
+    typedef Iterator<const double, const vector> const_iterator;
+    /// \brief vectorMutable::iterator
+    typedef Iterator<double&, vectorMutable> iterator;
     /**
        \brief The default constructor.
     */
+    vector(void);
+    /**
+       \brief Allocates memory space for \a n elements, and sets all
+       elements to \a init_value.
+       \throw GSL_error if memory allocation fails.
+    */
+    vector(size_t n, double init_value=0);
+    /**
+       \brief The copy constructor.
+       \note If the object to be copied is a vector view, the values
+       of the view will be copied, i.e. the view is not copied.
+       \throw A GSL_error is indirectly thrown if memory allocation
+       fails.
+    */
+    vector(const vector& other);
+    /**
+       \brief Vector view constructor.
+       Create a view of vector \a v, with starting index \a offset,
+       size \a n, and an optional \a stride.
+       A vector view can be used as any vector with the difference
+       that changes made to the view will also change the object that
+       is viewed. Also, using the copy constructor will create a new
+       vector object that is a copy of whatever is viewed. If a copy
+       of the view is needed then you should use this constructor to
+       obtain a copy.
+       \note If the object viewed by the view goes out of scope or is
+       deleted, the view becomes invalid and the result of further use
+       is undefined.
+       \throw GSL_error if a view cannot be set up.
+    */
+    vector(vector& v, size_t offset, size_t n, size_t stride=1);
+    /**
+       \brief Vector const view constructor.
+       Create a view of vector \a v, with starting index \a offset,
+       size \a n, and an optional \a stride.
+       A vector view can be used as any const vector. Using the copy
+       constructor will create a new vector object that is a copy of
+       whatever is viewed. If a copy of the view is needed then you
+       should use this constructor to obtain a copy.
+       \note If the object viewed by the view goes out of scope or is
+       deleted, the view becomes invalid and the result of further use
+       is undefined.
+       \throw GSL_error if a view cannot be set up.
+    */
+    vector(const vector& v, size_t offset, size_t n, size_t stride=1);
+    ///
+    /// Matrix row/column view constructor.
+    ///
+    /// Create a row/column vector view of matrix \a m, pointing at
+    /// row/column \a i. The parameter \a row is used to set whether
+    /// the view should be a row or column view. If \a row is set to
+    /// true, the view will be a row view (default behaviour), and,
+    /// naturally, a column view otherwise.
+    ///
+    /// A vector view can be used as any vector with the difference
+    /// that changes made to the view will also change the object that
+    /// is viewed. Also, using the copy constructor will create a new
+    /// vector object that is a copy of whatever is viewed. If a copy
+    /// of the view is needed then you should use the vector view
+    /// constructor to obtain a copy.
+    ///
+    /// @note If the object viewed by the view goes out of scope or is
+    /// deleted, the view becomes invalid and the result of further
+    /// use is undefined.
+    ///
+    vector(matrix& m, size_t i, bool row=true);
+    ///
+    /// Matrix row/column const view constructor.
+    ///
+    /// Create a row/column vector view of matrix \a m, pointing at
+    /// row/column \a i. The parameter \a row is used to set whether
+    /// the view should be a row or column view. If \a row is set to
+    /// true, the view will be a row view (default behaviour), and,
+    /// naturally, a column view otherwise.
+    ///
+    /// A const vector view can be used as any const vector. Using the
+    /// copy constructor will create a new vector object that is a
+    /// copy of whatever is viewed. If a copy of the view is needed
+    /// then you should use the vector view constructor to obtain a
+    /// copy.
+    ///
+    /// @note If the object viewed by the view goes out of scope or is
+    /// deleted, the view becomes invalid and the result of further
+    /// use is undefined.
+    ///
+    vector(const matrix& m, size_t i, bool row=true);
+    /**
+       \brief The istream constructor.
+       Either elements should be separated with white space characters
+       (default), or elements should be separated by the delimiter \a
+       sep. When delimiter \a sep is used empty elements are stored as
+       NaN's (except that empty lines are ignored). The end of input
+       to the vector is at end of file marker.
+       \throw GSL_error if memory allocation fails, IO_error if
+       unexpected input is found in the input stream.
+    */
+    explicit vector(std::istream &, char sep='\0')
+      throw(utility::IO_error, std::exception);
+    vectorMutable(gsl_vector*);
     ///
     /// The destructor.
     ///
+    ~vector(void);
+    /**
+       \return mutable iterator to start of vector
+    virtual ~vectorMutable(void);
+    // to enable overloading from base class
+    using vectorBase::begin;
+    /**
+       \return mutable iterator to start of vectorMutable
      */
     iterator begin(void);
     /**
-       \return read-only iterator to start of vector
-     */
-    const_iterator begin(void) const;
-    /**
-       \brief Make a copy of \a other.
-       This function will make a deep copy of \a other. Memory is
-       resized and view state is changed if needed.
-    */
-    const vector& clone(const vector& other);
-    /**
        \brief This function performs element-wise division, \f$ this_i =
        this_i/other_i \; \forall i \f$.
 …
        \throw GSL_error if dimensions mis-match.
     */
+    void div(const vector& other);
+    /**
+       \return mutable iterator to end of vector
+    void div(const vectorBase& other);
+    // to enable overloading from base class
+    using vectorBase::end;
+    /**
+       \return mutable iterator to end of vectorMutable
      */
     iterator end(void);
+    /**
+       \return read-only iterator to end of vector
+     */
+    const_iterator end(void) const;
+    /**
+       \brief Check whether vectors are equal within a user defined
+       precision, set by \a precision.
+       \return True if each element deviates less or equal than \a
+       d. If any vector contain a NaN, false is always returned.
+       \see operator== and operator!=
+    */
+    bool equal(const vector&, const double precision=0) const;
+    ///
+    /// @return A const pointer to the internal GSL vector,
+    ///
+    const gsl_vector* gsl_vector_p(void) const;
+    // to enable overloading from base class
+    using vectorBase::gsl_vector_p;
     ///
 …
     ///
     gsl_vector* gsl_vector_p(void);
-    ///
-    /// Check if the vector object is a view (sub-vector) to another
-    /// vector.
-    ///
-    /// @return True if the object is a view, false othwerwise.
-    ///
-    bool isview(void) const;
     /**
 …
        \throw GSL_error if dimensions mis-match.
     */
+    void mul(const vector& other);
+    /**
+       @brief Resize vector
+       All elements are set to @a init_value.
+       \note Underlying GSL vector is destroyed and a view into this
+       vector becomes invalid.
+    */
+    void resize(size_t, double init_value=0);
+    /**
+       \brief Reverse the order of elements in the vector.
+    void mul(const vectorBase& other);
+    /**
+       \brief Reverse the order of elements in the vectorMutable.
     */
     void reverse(void);
 …
     void all(const double& value);
-    ///
-    /// @return the number of elements in the vector.
-    ///
-    size_t size(void) const;
     /**
        \brief Exchange elements \a i and \a j.
        \throw GSL_error if vector lengths differs.
+       \throw GSL_error if vectorMutable lengths differs.
     */
     void swap(size_t i, size_t j);
+    // to allow overloading from base class
+    using vectorBase::operator();
     /**
 …
     /**
-       \brief Element access operator.
-       \return Const reference to element \a i.
-       \throw If GSL range checks are enabled in the underlying GSL
-       library a GSL_error exception is thrown if either index is out
-       of range.
-    */
-    const double& operator()(size_t i) const;
-    ///
-    /// Element access operator.
-    ///
-    /// @return Reference to element \a i.
-    ///
-    double& operator[](size_t i);
-    ///
-    /// Const element access operator.
-    ///
-    /// @return The value of element \a i.
-    ///
-    const double& operator[](size_t i) const;
-    /**
-       \brief Comparison operator. Takes linear time.
-       Checks are performed with exact matching, i.e., rounding off
-       effects may destroy comparison. Use the equal function for
-       comparing elements within a user defined precision.
-       \return True if all elements are equal otherwise false.
-       \see equal(const vector&, const double precision=0)
-    */
-    bool operator==(const vector&) const;
-    /**
-       \brief Comparison operator. Takes linear time.
-       Checks are performed with exact matching, i.e., rounding off
-       effects may destroy comparison. Use the equal function for
-       comparing elements within a user defined precision.
-       \return False if all elements are equal otherwise true.
-       \see equal(const vector&, const double precision=0)
-    */
-    bool operator!=(const vector&) const;
-    ///
-    /// @return The dot product.
-    ///
-    double operator*(const vector&) const;
-    /**
        \brief The assignment operator.
        Dimensions of the vectors must match. If the LHS vector is a
+       Dimensions of the vectorMutables must match. If the LHS vectorMutable is a
        view, the underlying data will be changed.
        \return A const reference to the resulting vector.
        \see void set(const vector&).
        \throw GSL_error if dimensions mis-match.
     */
     const vector& operator=(const vector&);
     /**
        \brief Addition and assign operator. Vector addition, \f$
+       \return A const reference to the resulting vectorMutable.
+       \see void set(const vectorMutable&).
+       \throw GSL_error if dimensions mis-match.
+    */
+    virtual const vectorMutable& operator=(vectorMutable&)=0;
+    /**
+       \brief Addition and assign operator. VectorMutable addition, \f$
        this_i = this_i + other_i \; \forall i \f$.
        \return A const reference to the resulting vector.
        \throw GSL_error if dimensions mis-match.
     */
     const vector& operator+=(const vector&);
     /**
        \brief Add a constant to a vector, \f$ this_i = this_i + d \;
+       \return A const reference to the resulting vectorMutable.
+       \throw GSL_error if dimensions mis-match.
+    */
+    const vectorMutable& operator+=(const vectorBase&);
+    /**
+       \brief Add a constant to a vectorMutable, \f$ this_i = this_i + d \;
        \forall i \f$.
        \return A const reference to the resulting vector.
     */
     const vector& operator+=(double d);
     /**
        \brief Subtract and assign operator. Vector subtraction, \f$
+       \return A const reference to the resulting vectorMutable.
+    */
+    const vectorMutable& operator+=(double d);
+    /**
+       \brief Subtract and assign operator. VectorMutable subtraction, \f$
        this_i = this_i - other_i \; \forall i \f$.
        \return A const reference to the resulting vector.
        \throw GSL_error if dimensions mis-match.
     */
     const vector& operator-=(const vector&);
     /**
        \brief Subtract a constant to a vector, \f$ this_i = this_i - d
+       \return A const reference to the resulting vectorMutable.
+       \throw GSL_error if dimensions mis-match.
+    */
+    const vectorMutable& operator-=(const vectorBase&);
+    /**
+       \brief Subtract a constant to a vectorMutable, \f$ this_i = this_i - d
        \; \forall i \f$.
        \return A const reference to the resulting vector.
     */
     const vector& operator-=(double d);
+       \return A const reference to the resulting vectorMutable.
+    */
+    const vectorMutable& operator-=(double d);
     /**
 …
        this_i * d \; \forall i \f$.
+       \return A const reference to the resulting vector.
+    */
+    const vector& operator*=(double d);
+  private:
+    /**
+       \brief Create a new copy of the internal GSL vector.
+       Necessary memory for the new GSL vector is allocated and the
+       caller is responsible for freeing the allocated memory.
+       \return A pointer to a copy of the internal GSL vector.
+       \throw GSL_error if memory cannot be allocated for the new
+       copy, or if dimensions mis-match.
+    */
+    gsl_vector* create_gsl_vector_copy(void) const;
+    /**
+       \brief Clear all dynamically allocated memory.
+       Internal utility function.
+    */
+    void delete_allocated_memory(void);
+    gsl_vector* v_;
+    gsl_vector_view* view_;
+    const gsl_vector_const_view* view_const_;
+    // proxy_v_ is used to access the proper underlying gsl_vector
+    // in all const member functions. It is not used by design for
+    // non-const vector functions and operators. This is to make sure
+    // that runtime errors occur if a const vector is used in an
+    // inappropriate manner such as on left hand side in assignment
+    // (remember, v_ is null for const vector views).
+    const gsl_vector* proxy_v_;
+       \return A const reference to the resulting vectorMutable.
+    */
+    const vectorMutable& operator*=(double d);
+  protected:
+    gsl_vector* vec_;
+    struct proxy
+    {
+      gsl_vector* vec_;
+    };
+  public:
+    /**
+     */
+    operator proxy();
   };
   /**
+     \brief Check if all elements of the vector are zero.
+     \return True if all elements in the vector is zero, false
+     othwerwise.
+  */
+  bool isnull(const vector&);
+  /**
+     \brief Get the maximum value of the vector.
+     \return The maximum value of the vector.
+  */
+  double max(const vector&);
+  /**
+     \brief Locate the maximum value in the vector.
+     \return The index to the maximum value of the vector.
+     \note Lower index has precedence.
+  */
+  size_t max_index(const vector&);
+  /**
+     \brief Get the minimum value of the vector.
+     \return The minimum value of the vector.
+  */
+  double min(const vector&);
+  /**
+     \brief Locate the minimum value in the vector.
+     \return The index to the minimum value of the vector.
+     \note Lower index has precedence.
+  */
+  size_t min_index(const vector&);
+  /**
+     \brief Create a vector \a flag indicating NaN's in another vector
+     \a templat.
+     The \a flag vector is changed to contain 1's and 0's only. A 1
+     means that the corresponding element in the \a templat vector is
+     valid and a zero means that the corresponding element is a NaN.
+     \note Space for vector \a flag is reallocated to fit the size of
+     vector \a templat if sizes mismatch.
+     \return True if the \a templat vector contains at least one NaN.
+  */
+  bool nan(const vector& templat, vector& flag);
+  /**
+     \brief Transforms a vector to a basis vector.
+     \brief Transforms a vectorMutable to a basis vector.
      All elements are set to zero except the \a i-th element which is
      set to one.
   */
+  void set_basis(vector&, size_t i);
+  /**
+     Randomly shuffles the elements in vector \a invec
+  */
+  void shuffle(vector& invec);
+  /**
+     Sort the elements in the vector.
+  */
+  void sort(vector&);
+  /**
+     Create a vector \a sort_index containing the indeces of
+     elements in a another vector \a invec.  The elements of \a
+     sort_index give the index of the vector element which would
+     have been stored in that position if the vector had been sorted
+     in place. The first element of \a sort_index gives the index of the least
+     element in \a invec, and the last element of \a sort_index gives the index of the
+     greatest element in \a invec . The vector \a invec is not changed.
+  */
+  void sort_index(std::vector<size_t>& sort_index, const vector& invec);
+  /** Similar to sort_index but creates a vector with indices to the \a k
+  smallest elements in \a invec.
+  */
+  void sort_smallest_index(std::vector<size_t>& sort_index, size_t k, const
+  vector& invec);
+  /** Similar to sort_index but creates a vector with indices to the \a k
+  largest elements in \a invec.
+  */
+  void sort_largest_index(std::vector<size_t>& sort_index, size_t k, const
+  vector& invec);
+  /**
+     \brief Calculate the sum of all vector elements.
+     \return The sum.
+  */
+  double sum(const vector&);
+  void set_basis(vectorMutable&, size_t i);
+  /**
+     Randomly shuffles the elements in vectorMutable \a invec
+  */
+  void shuffle(vectorMutable& invec);
+  /**
+     Sort the elements in the vectorMutable.
+  */
+  void sort(vectorMutable&);
   /**
      \brief Swap vector elements by copying.
+     The two vectors must have the same length.
+     \throw GSL_error if vector lengths differs.
+  */
+  void swap(vector&, vector&);
+  /**
+     \brief The output operator for the vector class.
+  */
+  std::ostream& operator<<(std::ostream&, const vector&);
+     The two vectorMutables must have the same length.
+     \throw GSL_error if vectorMutable lengths differs.
+  */
+  //void swap(vectorMutable&, vectorMutable&);
 }}} // of namespace utility, yat, and theplu

branches/peter-dev/yat/utility/vectorView.cc

-                      r994
+                      r995
 */
 #include "vector.h"
+#include "vectorView.h"
 #include "matrix.h"
 #include "utility.h"
 …
   vector::vector(void)
     : v_(NULL), view_(NULL), view_const_(NULL), proxy_v_(NULL)
+  vectorView::vectorView(void)
+    : vectorMutable(NULL)
+  {
+  }
+  vector::vector(size_t n, double init_value)
+    : v_(gsl_vector_alloc(n)), view_(NULL), view_const_(NULL),
+      proxy_v_(v_)
+  vectorView::vectorView(vectorMutable& other)
+    : vectorMutable(NULL), view_(NULL)
+  {
+    if (!v_)
+      throw utility::GSL_error("vector::vector failed to allocate memory");
+    all(init_value);
+    view_ = new gsl_vector_view(gsl_vector_subvector(other.gsl_vector_p(),0,
+                                                     other.size()));
+    const_vec_ = vec_ = &(view_->vector);
+  }
+  vector::vector(const vector& other)
+    : v_(other.create_gsl_vector_copy()), view_(NULL),
+      view_const_(NULL), proxy_v_(v_)
+  vectorView::vectorView(vectorMutable& v, size_t offset,size_t n,size_t stride)
+    : vectorMutable(NULL)
+  {
+    view_ =
+      new gsl_vector_view(gsl_vector_subvector_with_stride(v.gsl_vector_p(),
+                                                           offset, stride,n));
+    const_vec_ = vec_ = &(view_->vector);
+  }
+  vector::vector(vector& v, size_t offset, size_t n, size_t stride)
+    : view_const_(NULL)
+  {
+    view_ = new gsl_vector_view(gsl_vector_subvector_with_stride(v.v_,offset,
+                                                                 stride,n));
+    if (!view_)
+      throw utility::GSL_error("vector::vector failed to setup view");
+    proxy_v_ = v_ = &(view_->vector);
+  }
+  vector::vector(const vector& v, size_t offset, size_t n, size_t stride)
+    : v_(NULL), view_(NULL)
+  {
+    view_const_ = new gsl_vector_const_view(
+                   gsl_vector_const_subvector_with_stride(v.v_,offset,stride,n));
+    if (!view_const_)
+      throw utility::GSL_error("vector::vector failed to setup view");
+    proxy_v_ = &(view_const_->vector);
+  }
+  vector::vector(matrix& m, size_t i, bool row)
+    : view_const_(NULL)
+  vectorView::vectorView(matrix& m, size_t i, bool row)
+    : vectorMutable(NULL)
+  {
     view_=new gsl_vector_view(row ?
                               gsl_matrix_row   (m.gsl_matrix_p(),i) :
                               gsl_matrix_column(m.gsl_matrix_p(),i)  );
+    if (!view_)
+      throw utility::GSL_error("vector::vector failed to setup view");
+    proxy_v_ = v_ = &(view_->vector);
+    const_vec_ = vec_ = &(view_->vector);
+  }
   vector::vector(const matrix& m, size_t i, bool row)
     : v_(NULL), view_(NULL)
+  vectorView::vectorView(proxy p)
+    : vectorMutable(NULL)
+  {
+    view_const_ = new gsl_vector_const_view(row ?
+                                   gsl_matrix_const_row   (m.gsl_matrix_p(),i) :
+                                   gsl_matrix_const_column(m.gsl_matrix_p(),i) );
+    if (!view_const_)
+      throw utility::GSL_error("vector::vector failed to setup view");
+    proxy_v_ = &(view_const_->vector);
+    view_ = new gsl_vector_view(gsl_vector_subvector(p.vec_, 0, p.vec_->size));
+    const_vec_ = vec_ = &(view_->vector);
+  }
+  vector::vector(std::istream& is, char sep)
+    throw (utility::IO_error, std::exception)
+    : view_(NULL), view_const_(NULL)
+  vectorView::~vectorView(void)
+  {
+    // read the data file and store in stl vectors (dynamically
+    // expandable)
+    std::vector<std::vector<double> > data_matrix;
+    u_int nof_columns=0;
+    u_int nof_rows=0;
+    std::string line;
+    while(getline(is, line, '\n')) {
+      // Empty lines
+      if (!line.size())
+          continue;
+      nof_rows++;
+      std::vector<double> v;
+      std::string element;
+      std::stringstream ss(line);
+      bool ok=true;
+      while(ok) {
+        if(sep=='\0')
+          ok=(ss>>element);
+        else
+          ok=getline(ss, element, sep);
+        if(!ok)
+          break;
+        if(utility::is_double(element)) {
+          v.push_back(atof(element.c_str()));
+        }
+        else if (!element.size() || utility::is_nan(element)) {
+          v.push_back(std::numeric_limits<double>::quiet_NaN());
+        }
+        else {
+          std::stringstream ss("Warning: '");
+          ss << element << "' is not a double.";
+          throw IO_error(ss.str());
+        }
+      }
+      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
+          v.push_back(std::numeric_limits<double>::quiet_NaN());
+      if (!nof_columns)
+        nof_columns=v.size();
+      else if (nof_rows && (nof_columns>1)) {
+        std::ostringstream s;
+        s << "vector::vector(std::istream&) data file error:\n"
+          << "    File has inconsistent number of rows (" << nof_rows
+          << ") and columns (" << nof_columns
+          << ").\n    Expected a row or a column vector.";
+        throw utility::IO_error(s.str());
+      }
+      else if (v.size()!=nof_columns) {
+        std::ostringstream s;
+        s << "vector::vector(std::istream&) data file error:\n"
+          << "    Line " << nof_rows << " has " << v.size()
+          << " columns; expected " << nof_columns << " column.";
+        throw utility::IO_error(s.str());
+      }
+      data_matrix.push_back(v);
+    }
+    // manipulate the state of the stream to be good
+    is.clear(std::ios::goodbit);
+    // convert the data to a gsl vector
+    proxy_v_ = v_ = gsl_vector_alloc(nof_rows*nof_columns);
+    if (!v_)
+      throw utility::GSL_error("vector::vector failed to allocate memory");
+    size_t n=0;
+    // if gsl error handler disabled, out of bounds index will not
+    // abort the program.
+    for (size_t i=0; i<nof_rows; i++)
+      for (size_t j=0; j<nof_columns; j++)
+        gsl_vector_set( v_, n++, data_matrix[i][j] );
+    clean_up();
+  }
   vector::~vector(void)
+  const vectorMutable& vectorView::operator=(vectorMutable& other )
+  {
+    delete_allocated_memory();
+  }
+  vector::iterator vector::begin(void)
+  {
+    return vector::iterator(*this, 0);
+  }
+  vector::const_iterator vector::begin(void) const
+  {
+    return vector::const_iterator(*this, 0);
+  }
+  const vector& vector::clone(const vector& other)
+  {
+    if (this!=&other) {
+      delete_allocated_memory();
+      if (other.view_) {
+        view_ = new gsl_vector_view(*other.view_);
+        proxy_v_ = v_ = &(view_->vector);
+      }
+      else if (other.view_const_) {
+        view_const_ = new gsl_vector_const_view(*other.view_const_);
+        proxy_v_ = &(view_const_->vector);
+        v_=NULL;
+      }
+      else if (other.v_)
+        proxy_v_ = v_ = other.create_gsl_vector_copy();
+    }
+    if (size()!=other.size())
+      throw utility::GSL_error("vectorView::dimension mis-match");
+    if (!other.size())
+      return *this;
+    clean_up();
+    view_ = new gsl_vector_view(gsl_vector_subvector(other.gsl_vector_p(),0,
+                                                     other.size()));
+    const_vec_ = vec_ = &view_->vector;
     return *this;
+  }
   gsl_vector* vector::create_gsl_vector_copy(void) const
+  const vectorMutable& vectorView::operator=(proxy p)
+  {
+    gsl_vector* vec = gsl_vector_alloc(size());
+    if (!vec)
+      throw utility::GSL_error("vector::create_gsl_vector_copy failed to allocate memory");
+    if (gsl_vector_memcpy(vec, proxy_v_))
+      throw utility::GSL_error("vector::create_gsl_matrix_copy dimension mis-match");
+    return vec;
+  }
+  void vector::delete_allocated_memory(void)
+  {
+    if (view_)
+      delete view_;
+    else if (view_const_)
+      delete view_const_;
+    else if (v_)
+      gsl_vector_free(v_);
+    proxy_v_=v_=NULL;
+  }
+  void vector::div(const vector& other)
+  {
+    assert(v_);
+    int status=gsl_vector_div(v_,other.gsl_vector_p());
+    if (status)
+      throw utility::GSL_error(std::string("vector::div",status));
+  }
+  vector::iterator vector::end(void)
+  {
+    return vector::iterator(*this, size());
+  }
+  vector::const_iterator vector::end(void) const
+  {
+    return vector::const_iterator(*this, size());
+  }
+  bool vector::equal(const vector& other, const double d) const
+  {
+    if (this==&other)
+      return true;
+    if (size()!=other.size())
+      return false;
+    // if gsl error handler disabled, out of bounds index will not
+    // abort the program.
+    for (size_t i=0; i<size(); ++i)
+      // The two last condition checks are needed for NaN detection
+      if (fabs( (*this)(i)-other(i) ) > d ||
+          (*this)(i)!=(*this)(i) || other(i)!=other(i))
+        return false;
+    return true;
+  }
+  const gsl_vector* vector::gsl_vector_p(void) const
+  {
+    return proxy_v_;
+  }
+  gsl_vector* vector::gsl_vector_p(void)
+  {
+    return v_;
+  }
+  bool vector::isview(void) const
+  {
+    return view_ || view_const_;
+  }
+  void vector::mul(const vector& other)
+  {
+    assert(v_);
+    int status=gsl_vector_mul(v_,other.gsl_vector_p());
+    if (status)
+      throw utility::GSL_error(std::string("vector::mul",status));
+  }
+  void vector::resize(size_t n, double init_value)
+  {
+    delete_allocated_memory();
+    proxy_v_ = v_ = gsl_vector_alloc(n);
+    if (!v_)
+      throw utility::GSL_error("vector::vector failed to allocate memory");
+    all(init_value);
+  }
+  void vector::reverse(void)
+  {
+    assert(v_);
+    gsl_vector_reverse(v_);
+  }
+  void vector::all(const double& value)
+  {
+    assert(v_);
+    gsl_vector_set_all(v_,value);
+  }
+  size_t vector::size(void) const
+  {
+    if (!proxy_v_)
+      return 0;
+    return proxy_v_->size;
+  }
+  void vector::swap(size_t i, size_t j)
+  {
+    assert(v_);
+    int status=gsl_vector_swap_elements(v_, i, j);
+    if (status)
+      throw utility::GSL_error(std::string("vector::swap_elements",status));
+  }
+  double& vector::operator()(size_t i)
+  {
+    assert(v_);
+    double* d=gsl_vector_ptr(v_, i);
+    if (!d)
+      throw utility::GSL_error("vector::operator()",GSL_EINVAL);
+    return *d;
+  }
+  const double& vector::operator()(size_t i) const
+  {
+    const double* d=gsl_vector_const_ptr(proxy_v_, i);
+    if (!d)
+      throw utility::GSL_error("vector::operator()",GSL_EINVAL);
+    return *d;
+  }
+  double& vector::operator[](size_t i)
+  {
+    return this->operator()(i);
+  }
+  const double& vector::operator[](size_t i) const
+  {
+    return this->operator()(i);
+  }
+  bool vector::operator==(const vector& other) const
+  {
+    return equal(other);
+  }
+  bool vector::operator!=(const vector& other) const
+  {
+    return !equal(other);
+  }
+  double vector::operator*( const vector &other ) const
+  {
+    double res = 0.0;;
+    for ( size_t i = 0; i < size(); ++i )
+      res += other(i) * (*this)(i);
+    return res;
+  }
+  const vector& vector::operator=( const vector& other )
+  {
+    assert(v_);
+    if (gsl_vector_memcpy(v_,other.gsl_vector_p()))
+      throw utility::GSL_error("vector::set dimension mis-match");
+    if (size()!=p.vec_->size)
+      throw utility::GSL_error("vectorView::dimension mis-match");
+    if (!size())
+      return *this;
+    clean_up();
+    view_ = new gsl_vector_view(gsl_vector_subvector(p.vec_,0, p.vec_->size));
+    const_vec_ = vec_ = &view_->vector;
     return *this;
+  }
   const vector& vector::operator+=(const vector& other)
+  void vectorView::clean_up(void)
+  {
+    assert(v_);
+    int status=gsl_vector_add(v_, other.gsl_vector_p());
+    if (status)
+      throw utility::GSL_error(std::string("vector::add", status));
+    return *this;
+    if (view_)
+      delete view_;
+  }
-  const vector& vector::operator+=(double d)
+  {
-    assert(v_);
-    gsl_vector_add_constant(v_, d);
-    return *this;
+  }
-  const vector& vector::operator-=(const vector& other)
+  {
-    assert(v_);
-    int status=gsl_vector_sub(v_, other.gsl_vector_p());
-    if (status)
-      throw utility::GSL_error(std::string("vector::sub", status));
-    return *this;
+  }
-  const vector& vector::operator-=(const double d)
+  {
-    assert(v_);
-    gsl_vector_add_constant(v_, -d);
-    return *this;
+  }
-  const vector& vector::operator*=(const double d)
+  {
-    assert(v_);
-    gsl_vector_scale(v_, d);
-    return *this;
+  }
-  bool isnull(const vector& v)
+  {
-    return gsl_vector_isnull(v.gsl_vector_p());
+  }
-  double max(const vector& v)
+  {
-    return gsl_vector_max(v.gsl_vector_p());
+  }
-  size_t max_index(const vector& v)
+  {
-    return gsl_vector_max_index(v.gsl_vector_p());
+  }
-  double min(const vector& v)
+  {
-    return gsl_vector_min(v.gsl_vector_p());
+  }
-  size_t min_index(const vector& v)
+  {
-    return gsl_vector_min_index(v.gsl_vector_p());
+  }
-  bool nan(const vector& templat, vector& flag)
+  {
-    size_t vsize=templat.size();
-    if (vsize!=flag.size())
-      flag.clone(vector(vsize,1.0));
-    else
-      flag.all(1.0);
-    bool nan=false;
-    for (size_t i=0; i<vsize; i++)
-      if (std::isnan(templat(i))) {
-        flag(i)=0;
-        nan=true;
+      }
-    return nan;
+  }
-  void set_basis(vector& v, size_t i)
+  {
-    assert(v.gsl_vector_p());
-    gsl_vector_set_basis(v.gsl_vector_p(),i);
+  }
-  void shuffle(vector& invec)
+  {
-    random::DiscreteUniform rnd;
-    std::random_shuffle(invec.begin(), invec.end(), rnd);
+  }
-  void sort(vector& v)
+  {
-    assert(v.gsl_vector_p());
-    std::sort(v.begin(), v.end());
+  }
-  void sort_index(std::vector<size_t>& sort_index, const vector& invec)
+  {
-    assert(invec.gsl_vector_p());
-    gsl_permutation* p = gsl_permutation_alloc(invec.size());
-    int status=gsl_sort_vector_index(p,invec.gsl_vector_p());
-    if (status) {
-      gsl_permutation_free(p);
-      throw utility::GSL_error(std::string("sort_index(vector&,const vector&)",status));
+    }
-    sort_index=std::vector<size_t>(p->data,p->data+p->size);
-    gsl_permutation_free(p);
+  }
-  void sort_smallest_index(std::vector<size_t>& sort_index, size_t k,
-                            const vector& invec)
+  {
-    assert(invec.gsl_vector_p());
-    assert(k<=invec.size());
-    sort_index.resize(k);
-    gsl_sort_vector_smallest_index(&sort_index[0],k,invec.gsl_vector_p());
+  }
-  void sort_largest_index(std::vector<size_t>& sort_index, size_t k,
-                            const vector& invec)
+  {
-    assert(invec.gsl_vector_p());
-    assert(k<=invec.size());
-    sort_index.resize(k);
-    gsl_sort_vector_largest_index(&sort_index[0],k,invec.gsl_vector_p());
+  }
-  double sum(const vector& v)
+  {
-    double sum = 0;
-    size_t vsize=v.size();
-    for (size_t i=0; i<vsize; ++i)
-      sum += v[i];
-    return sum;
+  }
-  void swap(vector& v, vector& w)
+  {
-    assert(v.gsl_vector_p()); assert(w.gsl_vector_p());
-    int status=gsl_vector_swap(v.gsl_vector_p(),w.gsl_vector_p());
-    if (status)
-      throw utility::GSL_error(std::string("swap(vector&,vector&)",status));
+  }
-  std::ostream& operator<<(std::ostream& s, const vector& a)
+  {
-    s.setf(std::ios::dec);
-    s.precision(12);
-    for (size_t j = 0; j < a.size(); ++j) {
-      s << a[j];
-      if ( (j+1)<a.size() )
-        s << s.fill();
+    }
-    return s;
+  }
 }}} // of namespace utility, yat, and thep

branches/peter-dev/yat/utility/vectorView.h

-                      r994
+                      r995
 #ifndef _theplu_yat_utility_vector_
 #define _theplu_yat_utility_vector_
+#ifndef _theplu_yat_utility_vector_view_
+#define _theplu_yat_utility_vector_view_
 // $Id$
 …
 */
+#include "vectorMutable.h"
 #include "Exception.h"
 #include "Iterator.h"
 …
   class matrix;
+  class vector;
   /**
 …
   */
   class vector
+  class vectorView : public vectorMutable
+  {
   public:
+    /// \brief vector::iterator
+    typedef Iterator<double&, vector> iterator;
+    /// \brief vector::const_iterator
+    typedef Iterator<const double, const vector> const_iterator;
+    /// \brief vectorView::iterator
+    typedef Iterator<double&, vectorView> iterator;
     /**
        \brief The default constructor.
+       \brief Default constructor.
     */
+    vector(void);
+    /**
+       \brief Allocates memory space for \a n elements, and sets all
+       elements to \a init_value.
+       \throw GSL_error if memory allocation fails.
+    */
+    vector(size_t n, double init_value=0);
+    vectorView(void);
     /**
        \brief The copy constructor.
-       \note If the object to be copied is a vector view, the values
-       of the view will be copied, i.e. the view is not copied.
-       \throw A GSL_error is indirectly thrown if memory allocation
-       fails.
     */
     vector(const vector& other);
+    vectorView(vectorMutable& other);
     /**
        \brief Vector view constructor.
+       \brief VectorView constructor.
        Create a view of vector \a v, with starting index \a offset,
+       Create a view of vectorView \a v, with starting index \a offset,
        size \a n, and an optional \a stride.
        A vector view can be used as any vector with the difference
+       A vectorView view can be used as any vectorView with the difference
        that changes made to the view will also change the object that
        is viewed. Also, using the copy constructor will create a new
        vector object that is a copy of whatever is viewed. If a copy
+       vectorView object that is a copy of whatever is viewed. If a copy
        of the view is needed then you should use this constructor to
        obtain a copy.
 …
        \throw GSL_error if a view cannot be set up.
     */
+    vector(vector& v, size_t offset, size_t n, size_t stride=1);
+    /**
+       \brief Vector const view constructor.
+       Create a view of vector \a v, with starting index \a offset,
+       size \a n, and an optional \a stride.
+       A vector view can be used as any const vector. Using the copy
+       constructor will create a new vector object that is a copy of
+       whatever is viewed. If a copy of the view is needed then you
+       should use this constructor to obtain a copy.
+       \note If the object viewed by the view goes out of scope or is
+       deleted, the view becomes invalid and the result of further use
+       is undefined.
+       \throw GSL_error if a view cannot be set up.
+    */
+    vector(const vector& v, size_t offset, size_t n, size_t stride=1);
+    vectorView(vectorMutable& v, size_t offset, size_t n, size_t stride=1);
     ///
     /// Matrix row/column view constructor.
     ///
     /// Create a row/column vector view of matrix \a m, pointing at
+    /// Create a row/column vectorView view of matrix \a m, pointing at
     /// row/column \a i. The parameter \a row is used to set whether
     /// the view should be a row or column view. If \a row is set to
 …
     /// naturally, a column view otherwise.
     ///
     /// A vector view can be used as any vector with the difference
+    /// A vectorView view can be used as any vectorView with the difference
     /// that changes made to the view will also change the object that
     /// is viewed. Also, using the copy constructor will create a new
     /// vector object that is a copy of whatever is viewed. If a copy
     /// of the view is needed then you should use the vector view
+    /// vectorView object that is a copy of whatever is viewed. If a copy
+    /// of the view is needed then you should use the vectorView view
     /// constructor to obtain a copy.
     ///
 …
     /// use is undefined.
     ///
+    vector(matrix& m, size_t i, bool row=true);
+    ///
+    /// Matrix row/column const view constructor.
+    ///
+    /// Create a row/column vector view of matrix \a m, pointing at
+    /// row/column \a i. The parameter \a row is used to set whether
+    /// the view should be a row or column view. If \a row is set to
+    /// true, the view will be a row view (default behaviour), and,
+    /// naturally, a column view otherwise.
+    ///
+    /// A const vector view can be used as any const vector. Using the
+    /// copy constructor will create a new vector object that is a
+    /// copy of whatever is viewed. If a copy of the view is needed
+    /// then you should use the vector view constructor to obtain a
+    /// copy.
+    ///
+    /// @note If the object viewed by the view goes out of scope or is
+    /// deleted, the view becomes invalid and the result of further
+    /// use is undefined.
+    ///
+    vector(const matrix& m, size_t i, bool row=true);
+    vectorView(matrix& m, size_t i, bool row=true);
     /**
+       \brief The istream constructor.
+       Either elements should be separated with white space characters
+       (default), or elements should be separated by the delimiter \a
+       sep. When delimiter \a sep is used empty elements are stored as
+       NaN's (except that empty lines are ignored). The end of input
+       to the vector is at end of file marker.
+       \throw GSL_error if memory allocation fails, IO_error if
+       unexpected input is found in the input stream.
+    */
+    explicit vector(std::istream &, char sep='\0')
+      throw(utility::IO_error, std::exception);
+     */
+    vectorView(proxy p);
     ///
     /// The destructor.
     ///
+    ~vector(void);
+    /**
+       \return mutable iterator to start of vector
+     */
+    iterator begin(void);
+    /**
+       \return read-only iterator to start of vector
+     */
+    const_iterator begin(void) const;
+    /**
+       \brief Make a copy of \a other.
+       This function will make a deep copy of \a other. Memory is
+       resized and view state is changed if needed.
+    */
+    const vector& clone(const vector& other);
+    /**
+       \brief This function performs element-wise division, \f$ this_i =
+       this_i/other_i \; \forall i \f$.
+       \throw GSL_error if dimensions mis-match.
+    */
+    void div(const vector& other);
+    /**
+       \return mutable iterator to end of vector
+     */
+    iterator end(void);
+    /**
+       \return read-only iterator to end of vector
+     */
+    const_iterator end(void) const;
+    /**
+       \brief Check whether vectors are equal within a user defined
+       precision, set by \a precision.
+       \return True if each element deviates less or equal than \a
+       d. If any vector contain a NaN, false is always returned.
+       \see operator== and operator!=
+    */
+    bool equal(const vector&, const double precision=0) const;
+    ///
+    /// @return A const pointer to the internal GSL vector,
+    ///
+    const gsl_vector* gsl_vector_p(void) const;
+    ///
+    /// @return A pointer to the internal GSL vector,
+    ///
+    gsl_vector* gsl_vector_p(void);
+    ///
+    /// Check if the vector object is a view (sub-vector) to another
+    /// vector.
+    ///
+    /// @return True if the object is a view, false othwerwise.
+    ///
+    bool isview(void) const;
+    /**
+       \brief This function performs element-wise multiplication, \f$
+       this_i = this_i * other_i \; \forall i \f$.
+       \throw GSL_error if dimensions mis-match.
+    */
+    void mul(const vector& other);
+    /**
+       @brief Resize vector
+       All elements are set to @a init_value.
+       \note Underlying GSL vector is destroyed and a view into this
+       vector becomes invalid.
+    */
+    void resize(size_t, double init_value=0);
+    /**
+       \brief Reverse the order of elements in the vector.
+    */
+    void reverse(void);
+    ///
+    /// Set all elements to \a value.
+    ///
+    void all(const double& value);
+    ///
+    /// @return the number of elements in the vector.
+    ///
+    size_t size(void) const;
+    /**
+       \brief Exchange elements \a i and \a j.
+       \throw GSL_error if vector lengths differs.
+    */
+    void swap(size_t i, size_t j);
+    /**
+       \brief Element access operator.
+       \return Reference to element \a i.
+       \throw If GSL range checks are enabled in the underlying GSL
+       library a GSL_error exception is thrown if either index is out
+       of range.
+    */
+    double& operator()(size_t i);
+    /**
+       \brief Element access operator.
+       \return Const reference to element \a i.
+       \throw If GSL range checks are enabled in the underlying GSL
+       library a GSL_error exception is thrown if either index is out
+       of range.
+    */
+    const double& operator()(size_t i) const;
+    ///
+    /// Element access operator.
+    ///
+    /// @return Reference to element \a i.
+    ///
+    double& operator[](size_t i);
+    ///
+    /// Const element access operator.
+    ///
+    /// @return The value of element \a i.
+    ///
+    const double& operator[](size_t i) const;
+    /**
+       \brief Comparison operator. Takes linear time.
+       Checks are performed with exact matching, i.e., rounding off
+       effects may destroy comparison. Use the equal function for
+       comparing elements within a user defined precision.
+       \return True if all elements are equal otherwise false.
+       \see equal(const vector&, const double precision=0)
+    */
+    bool operator==(const vector&) const;
+    /**
+       \brief Comparison operator. Takes linear time.
+       Checks are performed with exact matching, i.e., rounding off
+       effects may destroy comparison. Use the equal function for
+       comparing elements within a user defined precision.
+       \return False if all elements are equal otherwise true.
+       \see equal(const vector&, const double precision=0)
+    */
+    bool operator!=(const vector&) const;
+    ///
+    /// @return The dot product.
+    ///
+    double operator*(const vector&) const;
+    ~vectorView(void);
     /**
        \brief The assignment operator.
-       Dimensions of the vectors must match. If the LHS vector is a
-       view, the underlying data will be changed.
-       \return A const reference to the resulting vector.
-       \see void set(const vector&).
-       \throw GSL_error if dimensions mis-match.
-    */
-    const vector& operator=(const vector&);
-    /**
-       \brief Addition and assign operator. Vector addition, \f$
-       this_i = this_i + other_i \; \forall i \f$.
        \return A const reference to the resulting vector.
 …
        \throw GSL_error if dimensions mis-match.
     */
     const vector& operator+=(const vector&);
+    const vectorMutable& operator=(vectorMutable&);
+    /**
+       \brief Add a constant to a vector, \f$ this_i = this_i + d \;
+       \forall i \f$.
+       \return A const reference to the resulting vector.
+    */
+    const vector& operator+=(double d);
+    /**
+       \brief Subtract and assign operator. Vector subtraction, \f$
+       this_i = this_i - other_i \; \forall i \f$.
+       \return A const reference to the resulting vector.
+       \throw GSL_error if dimensions mis-match.
+    */
+    const vector& operator-=(const vector&);
+    /**
+       \brief Subtract a constant to a vector, \f$ this_i = this_i - d
+       \; \forall i \f$.
+       \return A const reference to the resulting vector.
+    */
+    const vector& operator-=(double d);
+    /**
+       \brief Multiply with scalar and assign operator, \f$ this_i =
+       this_i * d \; \forall i \f$.
+       \return A const reference to the resulting vector.
+    */
+    const vector& operator*=(double d);
+    const vectorMutable& operator=(proxy);
   private:
+    void clean_up(void);
-    /**
-       \brief Create a new copy of the internal GSL vector.
-       Necessary memory for the new GSL vector is allocated and the
-       caller is responsible for freeing the allocated memory.
-       \return A pointer to a copy of the internal GSL vector.
-       \throw GSL_error if memory cannot be allocated for the new
-       copy, or if dimensions mis-match.
-    */
-    gsl_vector* create_gsl_vector_copy(void) const;
-    /**
-       \brief Clear all dynamically allocated memory.
-       Internal utility function.
-    */
-    void delete_allocated_memory(void);
-    gsl_vector* v_;
     gsl_vector_view* view_;
+    const gsl_vector_const_view* view_const_;
+    // proxy_v_ is used to access the proper underlying gsl_vector
+    // in all const member functions. It is not used by design for
+    // non-const vector functions and operators. This is to make sure
+    // that runtime errors occur if a const vector is used in an
+    // inappropriate manner such as on left hand side in assignment
+    // (remember, v_ is null for const vector views).
+    const gsl_vector* proxy_v_;
   };
-  /**
-     \brief Check if all elements of the vector are zero.
-     \return True if all elements in the vector is zero, false
-     othwerwise.
-  */
-  bool isnull(const vector&);
-  /**
-     \brief Get the maximum value of the vector.
-     \return The maximum value of the vector.
-  */
-  double max(const vector&);
-  /**
-     \brief Locate the maximum value in the vector.
-     \return The index to the maximum value of the vector.
-     \note Lower index has precedence.
-  */
-  size_t max_index(const vector&);
-  /**
-     \brief Get the minimum value of the vector.
-     \return The minimum value of the vector.
-  */
-  double min(const vector&);
-  /**
-     \brief Locate the minimum value in the vector.
-     \return The index to the minimum value of the vector.
-     \note Lower index has precedence.
-  */
-  size_t min_index(const vector&);
-  /**
-     \brief Create a vector \a flag indicating NaN's in another vector
-     \a templat.
-     The \a flag vector is changed to contain 1's and 0's only. A 1
-     means that the corresponding element in the \a templat vector is
-     valid and a zero means that the corresponding element is a NaN.
-     \note Space for vector \a flag is reallocated to fit the size of
-     vector \a templat if sizes mismatch.
-     \return True if the \a templat vector contains at least one NaN.
-  */
-  bool nan(const vector& templat, vector& flag);
-  /**
-     \brief Transforms a vector to a basis vector.
-     All elements are set to zero except the \a i-th element which is
-     set to one.
-  */
-  void set_basis(vector&, size_t i);
-  /**
-     Randomly shuffles the elements in vector \a invec
-  */
-  void shuffle(vector& invec);
-  /**
-     Sort the elements in the vector.
-  */
-  void sort(vector&);
-  /**
-     Create a vector \a sort_index containing the indeces of
-     elements in a another vector \a invec.  The elements of \a
-     sort_index give the index of the vector element which would
-     have been stored in that position if the vector had been sorted
-     in place. The first element of \a sort_index gives the index of the least
-     element in \a invec, and the last element of \a sort_index gives the index of the
-     greatest element in \a invec . The vector \a invec is not changed.
-  */
-  void sort_index(std::vector<size_t>& sort_index, const vector& invec);
-  /** Similar to sort_index but creates a vector with indices to the \a k
-  smallest elements in \a invec.
-  */
-  void sort_smallest_index(std::vector<size_t>& sort_index, size_t k, const
-  vector& invec);
-  /** Similar to sort_index but creates a vector with indices to the \a k
-  largest elements in \a invec.
-  */
-  void sort_largest_index(std::vector<size_t>& sort_index, size_t k, const
-  vector& invec);
-  /**
-     \brief Calculate the sum of all vector elements.
-     \return The sum.
-  */
-  double sum(const vector&);
-  /**
-     \brief Swap vector elements by copying.
-     The two vectors must have the same length.
-     \throw GSL_error if vector lengths differs.
-  */
-  void swap(vector&, vector&);
-  /**
-     \brief The output operator for the vector class.
-  */
-  std::ostream& operator<<(std::ostream&, const vector&);
 }}} // of namespace utility, yat, and theplu

Context Navigation

Legend:

branches/peter-dev/test/vector_test.cc

branches/peter-dev/yat/classifier/NBC.cc

branches/peter-dev/yat/classifier/SVM.cc

branches/peter-dev/yat/classifier/utility.cc

branches/peter-dev/yat/regression/Local.cc

branches/peter-dev/yat/regression/MultiDimensional.cc

branches/peter-dev/yat/regression/MultiDimensionalWeighted.cc

branches/peter-dev/yat/statistics/AveragerPairWeighted.h

branches/peter-dev/yat/utility/Makefile.am

branches/peter-dev/yat/utility/PCA.cc

branches/peter-dev/yat/utility/matrix.cc

branches/peter-dev/yat/utility/matrix.h

branches/peter-dev/yat/utility/vector.cc

branches/peter-dev/yat/utility/vector.h

branches/peter-dev/yat/utility/vectorBase.cc

branches/peter-dev/yat/utility/vectorBase.h

branches/peter-dev/yat/utility/vectorConstView.cc

branches/peter-dev/yat/utility/vectorConstView.h

branches/peter-dev/yat/utility/vectorMutable.cc

branches/peter-dev/yat/utility/vectorMutable.h

branches/peter-dev/yat/utility/vectorView.cc

branches/peter-dev/yat/utility/vectorView.h

Download in other formats: