Changeset 995


Ignore:
Timestamp:
Nov 30, 2007, 11:55:30 PM (14 years ago)
Author:
Peter
Message:

Splitting vector class into three different classes for vector, view and const_view, and adding abstract classes appropriately. NOTE, the tests are not going through, not even compiling; sorry about that

Location:
branches/peter-dev
Files:
14 edited
8 copied

Legend:

Unmodified
Added
Removed
  • branches/peter-dev/test/vector_test.cc

    r988 r995  
    3434
    3535using namespace theplu::yat;
     36
    3637
    3738void check_file_access(std::string& str)
  • branches/peter-dev/yat/classifier/NBC.cc

    r963 r995  
    205205
    206206    // normalize each row (label) to sum up to unity (probability)
    207     for (size_t i=0; i<prediction.rows(); ++i)
    208       utility::vector(prediction,i) *=
    209         1.0/utility::sum(utility::vector(prediction,i));
    210 
     207    for (size_t i=0; i<prediction.rows(); ++i){
     208      prediction.row_vec(i) *=
     209        1.0/sum(prediction.row_vec(i));
     210    }
    211211  }
    212212
  • branches/peter-dev/yat/classifier/SVM.cc

    r964 r995  
    198198  {
    199199    trained_=false;
    200     alpha_.clone(utility::vector(target_.size(), 0));
     200    alpha_ = utility::vector(target_.size(), 0);
    201201  }
    202202
  • branches/peter-dev/yat/classifier/utility.cc

    r865 r995  
    3535  void convert(const DataLookup1D& lookup, utility::vector& vector)
    3636  {
    37     vector.clone(utility::vector(lookup.size()));
     37    vector = utility::vector(lookup.size());
    3838    for(u_int i=0; i<lookup.size(); i++)
    3939      vector(i)=lookup(i);
     
    4444  {
    4545   
    46     value.clone(utility::vector(lookup.size()));
    47     weight.clone(utility::vector(lookup.size()));
     46    value = utility::vector(lookup.size());
     47    weight = utility::vector(lookup.size());
    4848    for(u_int i=0; i<lookup.size(); i++){
    4949      value(i)=lookup.data(i);
  • branches/peter-dev/yat/regression/Local.cc

    r865 r995  
    2727#include "OneDimensionalWeighted.h"
    2828#include "yat/utility/vector.h"
     29#include "yat/utility/vectorView.h"
    2930
    3031#include <algorithm>
     
    5960
    6061    size_t nof_fits=data_.size()/step_size;
    61     x_.clone(utility::vector(nof_fits));
    62     y_predicted_.clone(utility::vector(x_.size()));
    63     y_err_.clone(utility::vector(x_.size()));
     62    x_ = utility::vector(nof_fits);
     63    y_predicted_ = utility::vector(x_.size());
     64    y_err_ = utility::vector(x_.size());
    6465    sort(data_.begin(), data_.end());
    6566
     
    102103      assert(max_index<data_.size());
    103104                               
    104       utility::vector x_local(x, min_index, max_index-min_index+1);
    105       utility::vector y_local(y, min_index, max_index-min_index+1);
     105      utility::vectorView x_local(x, min_index, max_index-min_index+1);
     106      utility::vectorView y_local(y, min_index, max_index-min_index+1);
    106107
    107108      // calculating weights
  • branches/peter-dev/yat/regression/MultiDimensional.cc

    r865 r995  
    5858    assert(x.rows()==y.size());
    5959    covariance_.clone(utility::matrix(x.columns(),x.columns()));
    60     fit_parameters_.clone(utility::vector(x.columns()));
     60    fit_parameters_ = utility::vector(x.columns());
    6161    if (work_)
    6262      gsl_multifit_linear_free(work_);
  • branches/peter-dev/yat/regression/MultiDimensionalWeighted.cc

    r916 r995  
    5959
    6060    covariance_.clone(utility::matrix(x.columns(),x.columns()));
    61     fit_parameters_.clone(utility::vector(x.columns()));
     61    fit_parameters_ = utility::vector(x.columns());
    6262    if (work_)
    6363      gsl_multifit_linear_free(work_);
  • branches/peter-dev/yat/statistics/AveragerPairWeighted.h

    r937 r995  
    201201  {
    202202    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)));
    204204      add(x[i],y[i],wx[i],wy[i]);
    205205    }
  • branches/peter-dev/yat/utility/Makefile.am

    r981 r995  
    2828  Alignment.cc ColumnStream.cc CommandLine.cc FileUtil.cc kNNI.cc \
    2929  matrix.cc NNI.cc Option.cc OptionFile.cc OptionHelp.cc OptionSwitch.cc \
    30   PCA.cc stl_utility.cc SVD.cc TypeInfo.cc utility.cc vector.cc WeNNI.cc
     30  PCA.cc stl_utility.cc SVD.cc TypeInfo.cc utility.cc vector.cc \
     31  vectorBase.cc vectorConstView.cc vectorMutable.cc vectorView.cc WeNNI.cc Peter.cc
    3132
    3233include_utilitydir = $(includedir)/yat/utility
     
    3738  IteratorWeighted.h kNNI.h matrix.h NNI.h \
    3839  Option.h OptionArg.h OptionFile.h OptionHelp.h OptionSwitch.h \
    39   PCA.h stl_utility.h SVD.h TypeInfo.h utility.h vector.h WeNNI.h \
    40   yat_assert.h
     40  PCA.h stl_utility.h SVD.h TypeInfo.h utility.h vector.h
     41  vectorBase.h vectorMutable.h vectorView.h vectorConstView.h
     42  WeNNI.h yat_assert.h
    4143
  • branches/peter-dev/yat/utility/PCA.cc

    r865 r995  
    2929#include "SVD.h"
    3030#include "utility.h"
     31#include "vectorConstView.h"
    3132
    3233#include <iostream>
     
    7374    eigenvectors_.clone(U);
    7475    eigenvectors_ .transpose();
     76    eigenvalues_ = pSVD->s();
     77
     78    // T
     79    for( size_t i = 0; i < eigenvalues_.size(); ++i )
     80      eigenvalues_(i) = eigenvalues_(i)*eigenvalues_(i);
     81    eigenvalues_ *= 1.0/A_center.rows();
     82
     83    // Sort the eigenvectors in order of eigenvalues
     84    // Simple (not efficient) algorithm that always
     85    // make sure that the i:th element is in its correct
     86    // position (N element --> Ordo( N*N ))
     87    for ( size_t i = 0; i < eigenvalues_.size(); ++i )
     88      for ( size_t j = i + 1; j < eigenvalues_.size(); ++j )
     89        if ( eigenvalues_(j) > eigenvalues_(i) ) {
     90          std::swap( eigenvalues_(i), eigenvalues_(j) );
     91          eigenvectors_.swap_rows(i,j);
     92        }
     93  }
     94
     95
     96  /*
     97  void PCA::process_transposed_problem(void)
     98  {
     99    // Row-center the data matrix
     100    utility::matrix A_center( A_.rows(), A_.columns() );
     101    this->row_center( A_center );
     102
     103    // Transform into SVD friendly dimension
     104    A_.transpose();
     105    A_center.transpose();
     106
     107    // Single value decompose the data matrix
     108    std::auto_ptr<SVD> pSVD( new SVD( A_center ) );
     109    pSVD->decompose();
     110    utility::matrix U(pSVD->U());
     111    utility::matrix V(pSVD->V());
     112
     113    // Read the eigenvectors and eigenvalues
     114    eigenvectors_.clone(V);
     115    eigenvectors_.transpose();
    75116    eigenvalues_.clone(pSVD->s());
     117
     118    // Transform back when done with SVD!
     119    // (used V insted of U now for eigenvectors)
     120    A_.transpose();
     121    A_center.transpose();
    76122
    77123    // T
     
    91137        }
    92138  }
    93 
    94 
    95   /*
    96   void PCA::process_transposed_problem(void)
    97   {
    98     // Row-center the data matrix
    99     utility::matrix A_center( A_.rows(), A_.columns() );
    100     this->row_center( A_center );
    101 
    102     // Transform into SVD friendly dimension
    103     A_.transpose();
    104     A_center.transpose();
    105 
    106     // Single value decompose the data matrix
    107     std::auto_ptr<SVD> pSVD( new SVD( A_center ) );
    108     pSVD->decompose();
    109     utility::matrix U(pSVD->U());
    110     utility::matrix V(pSVD->V());
    111 
    112     // Read the eigenvectors and eigenvalues
    113     eigenvectors_.clone(V);
    114     eigenvectors_.transpose();
    115     eigenvalues_.clone(pSVD->s());
    116 
    117     // Transform back when done with SVD!
    118     // (used V insted of U now for eigenvectors)
    119     A_.transpose();
    120     A_center.transpose();
    121 
    122     // T
    123     for( size_t i = 0; i < eigenvalues_.size(); ++i )
    124       eigenvalues_[ i ] = eigenvalues_[ i ]*eigenvalues_[ i ];
    125     eigenvalues_ *= 1.0/A_center.rows();
    126 
    127     // Sort the eigenvectors in order of eigenvalues
    128     // Simple (not efficient) algorithm that always
    129     // make sure that the i:th element is in its correct
    130     // position (N element --> Ordo( N*N ))
    131     for ( size_t i = 0; i < eigenvalues_.size(); ++i )
    132       for ( size_t j = i + 1; j < eigenvalues_.size(); ++j )
    133         if ( eigenvalues_[ j ] > eigenvalues_[ i ] ) {
    134           std::swap( eigenvalues_[ i ], eigenvalues_[ j ] );
    135           eigenvectors_.swap_rows(i,j);
    136         }
    137   }
    138139  */
    139140
     
    189190  void PCA::row_center(utility::matrix& A_center)
    190191  {
    191     meanvalues_.clone(utility::vector(A_.rows()));
     192    meanvalues_ = vector(A_.rows());
    192193    utility::vector A_row_sum(A_.rows());
    193     for (size_t i=0; i<A_row_sum.size(); ++i)
    194       A_row_sum(i)=utility::sum(utility::vector(A_,i));
     194    for (size_t i=0; i<A_row_sum.size(); ++i){
     195      vectorConstView tmp;
     196      tmp = A_.const_row(i);
     197      A_row_sum(i) = sum(tmp);
     198    }
    195199    for( size_t i = 0; i < A_center.rows(); ++i ) {
    196       meanvalues_[i] = A_row_sum(i) / static_cast<double>(A_.columns());
     200      meanvalues_(i) = A_row_sum(i) / static_cast<double>(A_.columns());
    197201      for( size_t j = 0; j < A_center.columns(); ++j )
    198202        A_center(i,j) = A_(i,j) - meanvalues_(i);
  • branches/peter-dev/yat/utility/matrix.cc

    r916 r995  
    2727#include "matrix.h"
    2828#include "vector.h"
     29#include "vectorView.h"
     30#include "vectorConstView.h"
    2931#include "utility.h"
    3032
     
    302304
    303305
    304   void matrix::column(const size_t column, const vector& vec)
    305   {
    306     assert(m_);
    307     int status=gsl_matrix_set_col(m_, column, vec.gsl_vector_p());
    308     if (status)
    309       throw utility::GSL_error(std::string("matrix::set_column",status));
    310   }
    311 
    312 
    313   void matrix::row(const size_t row, const vector& vec)
    314   {
    315     assert(m_);
    316     int status=gsl_matrix_set_row(m_, row, vec.gsl_vector_p());
    317     if (status)
    318       throw utility::GSL_error(std::string("matrix::set_row",status));
     306  vectorView matrix::column_vec(size_t col)
     307  {
     308    vectorView res(*this, col, false);
     309    return res;
     310  }
     311
     312
     313  vectorConstView matrix::const_column(size_t col) const
     314  {
     315    return vectorConstView(*this, col, false);
     316  }
     317
     318
     319  vectorConstView matrix::const_row(size_t col) const
     320  {
     321    return vectorConstView(*this, col, true);
     322  }
     323
     324
     325  vectorView matrix::row_vec(size_t row)
     326  {
     327    vectorView res(*this, row, true);
     328    return res;
    319329  }
    320330
     
    549559
    550560
    551   vector operator*(const matrix& m, const vector& v)
     561  vector operator*(const matrix& m, const vectorBase& v)
    552562  {
    553563    utility::vector res(m.rows());
    554564    for (size_t i=0; i<res.size(); ++i)
    555       res(i) = vector(m,i) * v;
     565      res(i) = vectorConstView(m,i) * v;
    556566    return res;
    557567  }
     
    562572    utility::vector res(m.columns());
    563573    for (size_t i=0; i<res.size(); ++i)
    564       res(i) = v * vector(m,i,false);
     574      res(i) = v * vectorConstView(m,i,false);
    565575    return res;
    566576  }
  • branches/peter-dev/yat/utility/matrix.h

    r865 r995  
    2929
    3030#include "vector.h"
     31#include "vectorView.h"
     32#include "vectorConstView.h"
    3133#include "Exception.h"
    3234
     
    141143    ~matrix(void);
    142144
     145    ///
     146    /// Set all elements to \a value.
     147    ///
     148    void all(const double value);
     149
    143150    /**
    144151       \brief Make a copy of \a other.
     
    149156    const matrix& clone(const matrix& other);
    150157
     158    /**
     159     */
     160    vectorView column_vec(size_t);
     161
    151162    ///
    152163    /// @return The number of columns in the matrix.
    153164    ///
    154165    size_t columns(void) const;
     166
     167    /**
     168     */
     169    vectorConstView const_column(size_t) const;
     170
     171    /**
     172     */
     173    vectorConstView const_row(size_t) const;
    155174
    156175    /**
     
    216235    size_t rows(void) const;
    217236
    218     ///
    219     /// Set all elements to \a value.
    220     ///
    221     void all(const double value);
    222 
    223     /**
    224        \brief Set \a column values to values in \a vec.
    225 
    226        \note No check on size is done.
    227 
    228        \throw GSL_error if index is out of range or mis-match in
    229        sizes.
    230     */
    231     void column(const size_t column, const vector& vec);
    232 
    233     /**
    234        \brief Set \a row values to values in \a vec.
    235 
    236        \note No check on size is done.
    237 
    238        \throw GSL_error if index is out of range or mis-match in
    239        sizes.
    240     */
    241     void row(const size_t row, const vector& vec);
     237    /**
     238     */
     239    vectorView row_vec(size_t);
    242240
    243241    /**
  • branches/peter-dev/yat/utility/vector.cc

    r978 r995  
    4444
    4545  vector::vector(void)
    46     : v_(NULL), view_(NULL), view_const_(NULL), proxy_v_(NULL)
     46    : vectorMutable(NULL)
    4747  {
    4848  }
     
    5050
    5151  vector::vector(size_t n, double init_value)
    52     : v_(gsl_vector_alloc(n)), view_(NULL), view_const_(NULL),
    53       proxy_v_(v_)
    54   {
    55     if (!v_)
     52    : vectorMutable(gsl_vector_alloc(n))
     53  {
     54    if (!vec_)
    5655      throw utility::GSL_error("vector::vector failed to allocate memory");
    5756
     
    6160
    6261  vector::vector(const vector& other)
    63     : v_(other.create_gsl_vector_copy()), view_(NULL),
    64       view_const_(NULL), proxy_v_(v_)
    65   {
    66   }
    67 
    68 
    69   vector::vector(vector& v, size_t offset, size_t n, size_t stride)
    70     : view_const_(NULL)
    71   {
    72     view_ = new gsl_vector_view(gsl_vector_subvector_with_stride(v.v_,offset,
    73                                                                  stride,n));
    74     if (!view_)
    75       throw utility::GSL_error("vector::vector failed to setup view");
    76     proxy_v_ = v_ = &(view_->vector);
    77   }
    78 
    79 
    80   vector::vector(const vector& v, size_t offset, size_t n, size_t stride)
    81     : v_(NULL), view_(NULL)
    82   {
    83     view_const_ = new gsl_vector_const_view(
    84                    gsl_vector_const_subvector_with_stride(v.v_,offset,stride,n));
    85     if (!view_const_)
    86       throw utility::GSL_error("vector::vector failed to setup view");
    87     proxy_v_ = &(view_const_->vector);
    88   }
    89 
    90 
    91   vector::vector(matrix& m, size_t i, bool row)
    92     : view_const_(NULL)
    93   {
    94     view_=new gsl_vector_view(row ?
    95                               gsl_matrix_row   (m.gsl_matrix_p(),i) :
    96                               gsl_matrix_column(m.gsl_matrix_p(),i)  );
    97     if (!view_)
    98       throw utility::GSL_error("vector::vector failed to setup view");
    99     proxy_v_ = v_ = &(view_->vector);
    100   }
    101 
    102 
    103   vector::vector(const matrix& m, size_t i, bool row)
    104     : v_(NULL), view_(NULL)
    105   {
    106     view_const_ = new gsl_vector_const_view(row ?
    107                                    gsl_matrix_const_row   (m.gsl_matrix_p(),i) :
    108                                    gsl_matrix_const_column(m.gsl_matrix_p(),i) );
    109     if (!view_const_)
    110       throw utility::GSL_error("vector::vector failed to setup view");
    111     proxy_v_ = &(view_const_->vector);
     62    : vectorMutable(create_gsl_vector_copy(other))
     63  {
     64  }
     65
     66
     67  vector::vector(const vectorBase& other)
     68    : vectorMutable(create_gsl_vector_copy(other))
     69  {
    11270  }
    11371
     
    11573  vector::vector(std::istream& is, char sep)
    11674    throw (utility::IO_error, std::exception)
    117     : view_(NULL), view_const_(NULL)
     75    : vectorMutable(NULL)
    11876  {
    11977    // read the data file and store in stl vectors (dynamically
     
    178136    is.clear(std::ios::goodbit);
    179137    // convert the data to a gsl vector
    180     proxy_v_ = v_ = gsl_vector_alloc(nof_rows*nof_columns);
    181     if (!v_)
     138    vec_ = gsl_vector_alloc(nof_rows*nof_columns);
     139    if (!vec_)
    182140      throw utility::GSL_error("vector::vector failed to allocate memory");
    183141    size_t n=0;
     
    186144    for (size_t i=0; i<nof_rows; i++)
    187145      for (size_t j=0; j<nof_columns; j++)
    188         gsl_vector_set( v_, n++, data_matrix[i][j] );
     146        gsl_vector_set( vec_, n++, data_matrix[i][j] );
    189147  }
    190148
     
    192150  vector::~vector(void)
    193151  {
    194     delete_allocated_memory();
    195   }
    196 
    197 
    198   vector::iterator vector::begin(void)
    199   {
    200     return vector::iterator(*this, 0);
    201   }
    202 
    203 
    204   vector::const_iterator vector::begin(void) const
    205   {
    206     return vector::const_iterator(*this, 0);
    207   }
    208 
    209 
    210   const vector& vector::clone(const vector& other)
    211   {
    212     if (this!=&other) {
    213 
    214       delete_allocated_memory();
    215 
    216       if (other.view_) {
    217         view_ = new gsl_vector_view(*other.view_);
    218         proxy_v_ = v_ = &(view_->vector);
    219       }
    220       else if (other.view_const_) {
    221         view_const_ = new gsl_vector_const_view(*other.view_const_);
    222         proxy_v_ = &(view_const_->vector);
    223         v_=NULL;
    224       }
    225       else if (other.v_)
    226         proxy_v_ = v_ = other.create_gsl_vector_copy();
    227 
     152    gsl_vector_free(vec_);
     153  }
     154
     155
     156  gsl_vector* vector::create_gsl_vector_copy(const vectorBase& other) const
     157  {
     158    gsl_vector* vec = gsl_vector_alloc(other.size());
     159    if (!vec)
     160      throw utility::GSL_error("vector::create_gsl_vector_copy failed to allocate memory");
     161    if (gsl_vector_memcpy(vec, other.gsl_vector_p()))
     162      throw utility::GSL_error("vector::create_gsl_matrix_copy dimension mis-match");
     163    return vec;
     164  }
     165
     166
     167  void vector::resize(size_t n, double init_value)
     168  {
     169    gsl_vector_free(vec_);
     170    vec_ = gsl_vector_alloc(n);
     171    if (!vec_)
     172      throw utility::GSL_error("vector::vector failed to allocate memory");
     173    all(init_value);
     174  }
     175
     176
     177  //Peter use swap idiom
     178  const vector& vector::operator=( const vectorBase& other )
     179  {
     180    if (!other.size())
     181      vec_=NULL;
     182    else  {   
     183      if (size()!=other.size())
     184        resize(other.size(),0.0);
     185      gsl_vector_memcpy(vec_,other.gsl_vector_p());
    228186    }
    229187    return *this;
     
    231189
    232190
    233   gsl_vector* vector::create_gsl_vector_copy(void) const
    234   {
    235     gsl_vector* vec = gsl_vector_alloc(size());
    236     if (!vec)
    237       throw utility::GSL_error("vector::create_gsl_vector_copy failed to allocate memory");
    238     if (gsl_vector_memcpy(vec, proxy_v_))
    239       throw utility::GSL_error("vector::create_gsl_matrix_copy dimension mis-match");
    240     return vec;
    241   }
    242 
    243 
    244   void vector::delete_allocated_memory(void)
    245   {
    246     if (view_)
    247       delete view_;
    248     else if (view_const_)
    249       delete view_const_;
    250     else if (v_)
    251       gsl_vector_free(v_);
    252     proxy_v_=v_=NULL;
    253   }
    254 
    255 
    256   void vector::div(const vector& other)
    257   {
    258     assert(v_);
    259     int status=gsl_vector_div(v_,other.gsl_vector_p());
    260     if (status)
    261       throw utility::GSL_error(std::string("vector::div",status));
    262   }
    263 
    264 
    265   vector::iterator vector::end(void)
    266   {
    267     return vector::iterator(*this, size());
    268   }
    269 
    270 
    271   vector::const_iterator vector::end(void) const
    272   {
    273     return vector::const_iterator(*this, size());
    274   }
    275 
    276 
    277   bool vector::equal(const vector& other, const double d) const
    278   {
    279     if (this==&other)
    280       return true;
    281     if (size()!=other.size())
    282       return false;
    283     // if gsl error handler disabled, out of bounds index will not
    284     // abort the program.
    285     for (size_t i=0; i<size(); ++i)
    286       // The two last condition checks are needed for NaN detection
    287       if (fabs( (*this)(i)-other(i) ) > d ||
    288           (*this)(i)!=(*this)(i) || other(i)!=other(i))
    289         return false;
    290     return true;
    291   }
    292 
    293 
    294   const gsl_vector* vector::gsl_vector_p(void) const
    295   {
    296     return proxy_v_;
    297   }
    298 
    299 
    300   gsl_vector* vector::gsl_vector_p(void)
    301   {
    302     return v_;
    303   }
    304 
    305 
    306   bool vector::isview(void) const
    307   {
    308     return view_ || view_const_;
    309   }
    310 
    311 
    312   void vector::mul(const vector& other)
    313   {
    314     assert(v_);
    315     int status=gsl_vector_mul(v_,other.gsl_vector_p());
    316     if (status)
    317       throw utility::GSL_error(std::string("vector::mul",status));
    318   }
    319 
    320 
    321   void vector::resize(size_t n, double init_value)
    322   {
    323     delete_allocated_memory();
    324     proxy_v_ = v_ = gsl_vector_alloc(n);
    325     if (!v_)
    326       throw utility::GSL_error("vector::vector failed to allocate memory");
    327     all(init_value);
    328   }
    329 
    330 
    331   void vector::reverse(void)
    332   {
    333     assert(v_);
    334     gsl_vector_reverse(v_);
    335   }
    336 
    337 
    338   void vector::all(const double& value)
    339   {
    340     assert(v_);
    341     gsl_vector_set_all(v_,value);
    342   }
    343 
    344 
    345   size_t vector::size(void) const
    346   {
    347     if (!proxy_v_)
    348       return 0;
    349     return proxy_v_->size;
    350   }
    351 
    352 
    353   void vector::swap(size_t i, size_t j)
    354   {
    355     assert(v_);
    356     int status=gsl_vector_swap_elements(v_, i, j);
    357     if (status)
    358       throw utility::GSL_error(std::string("vector::swap_elements",status));
    359   }
    360 
    361 
    362   double& vector::operator()(size_t i)
    363   {
    364     assert(v_);
    365     double* d=gsl_vector_ptr(v_, i);
    366     if (!d)
    367       throw utility::GSL_error("vector::operator()",GSL_EINVAL);
    368     return *d;
    369   }
    370 
    371 
    372   const double& vector::operator()(size_t i) const
    373   {
    374     const double* d=gsl_vector_const_ptr(proxy_v_, i);
    375     if (!d)
    376       throw utility::GSL_error("vector::operator()",GSL_EINVAL);
    377     return *d;
    378   }
    379 
    380 
    381   double& vector::operator[](size_t i)
    382   {
    383     return this->operator()(i);
    384   }
    385 
    386 
    387   const double& vector::operator[](size_t i) const
    388   {
    389     return this->operator()(i);
    390   }
    391 
    392 
    393   bool vector::operator==(const vector& other) const
    394   {
    395     return equal(other);
    396   }
    397 
    398 
    399   bool vector::operator!=(const vector& other) const
    400   {
    401     return !equal(other);
    402   }
    403 
    404 
    405   double vector::operator*( const vector &other ) const
    406   {
    407     double res = 0.0;;
    408     for ( size_t i = 0; i < size(); ++i )
    409       res += other(i) * (*this)(i);
    410     return res;
    411   }
    412 
    413 
    414   const vector& vector::operator=( const vector& other )
    415   {
    416     assert(v_);
    417     if (gsl_vector_memcpy(v_,other.gsl_vector_p()))
    418       throw utility::GSL_error("vector::set dimension mis-match");
     191  //Peter use swap idiom
     192  const vectorMutable& vector::operator=( vectorMutable& other )
     193  {
     194    if (!other.size())
     195      vec_=NULL;
     196    else  {   
     197      if (size()!=other.size())
     198        resize(other.size(),0.0);
     199      gsl_vector_memcpy(vec_,other.gsl_vector_p());
     200    }
    419201    return *this;
    420202  }
    421 
    422 
    423   const vector& vector::operator+=(const vector& other)
    424   {
    425     assert(v_);
    426     int status=gsl_vector_add(v_, other.gsl_vector_p());
    427     if (status)
    428       throw utility::GSL_error(std::string("vector::add", status));
    429     return *this;
    430   }
    431 
    432 
    433   const vector& vector::operator+=(double d)
    434   {
    435     assert(v_);
    436     gsl_vector_add_constant(v_, d);
    437     return *this;
    438   }
    439 
    440 
    441   const vector& vector::operator-=(const vector& other)
    442   {
    443     assert(v_);
    444     int status=gsl_vector_sub(v_, other.gsl_vector_p());
    445     if (status)
    446       throw utility::GSL_error(std::string("vector::sub", status));
    447     return *this;
    448   }
    449 
    450 
    451   const vector& vector::operator-=(const double d)
    452   {
    453     assert(v_);
    454     gsl_vector_add_constant(v_, -d);
    455     return *this;
    456   }
    457 
    458 
    459   const vector& vector::operator*=(const double d)
    460   {
    461     assert(v_);
    462     gsl_vector_scale(v_, d);
    463     return *this;
    464   }
    465 
    466 
    467   bool isnull(const vector& v)
    468   {
    469     return gsl_vector_isnull(v.gsl_vector_p());
    470   }
    471 
    472 
    473   double max(const vector& v)
    474   {
    475     return gsl_vector_max(v.gsl_vector_p());
    476   }
    477 
    478 
    479   size_t max_index(const vector& v)
    480   {
    481     return gsl_vector_max_index(v.gsl_vector_p());
    482   }
    483 
    484 
    485   double min(const vector& v)
    486   {
    487     return gsl_vector_min(v.gsl_vector_p());
    488   }
    489 
    490 
    491   size_t min_index(const vector& v)
    492   {
    493     return gsl_vector_min_index(v.gsl_vector_p());
    494   }
    495 
    496 
    497   bool nan(const vector& templat, vector& flag)
    498   {
    499     size_t vsize=templat.size();
    500     if (vsize!=flag.size())
    501       flag.clone(vector(vsize,1.0));
    502     else
    503       flag.all(1.0);
    504     bool nan=false;
    505     for (size_t i=0; i<vsize; i++)
    506       if (std::isnan(templat(i))) {
    507         flag(i)=0;
    508         nan=true;
    509       }
    510     return nan;
    511   }
    512 
    513 
    514   void set_basis(vector& v, size_t i)
    515   {
    516     assert(v.gsl_vector_p());
    517     gsl_vector_set_basis(v.gsl_vector_p(),i);
    518   }
    519 
    520 
    521   void shuffle(vector& invec)
    522   {
    523     random::DiscreteUniform rnd;
    524     std::random_shuffle(invec.begin(), invec.end(), rnd);
    525   }
    526 
    527 
    528   void sort(vector& v)
    529   {
    530     assert(v.gsl_vector_p());
    531     std::sort(v.begin(), v.end());
    532   }
    533 
    534 
    535   void sort_index(std::vector<size_t>& sort_index, const vector& invec)
    536   {
    537     assert(invec.gsl_vector_p());
    538     gsl_permutation* p = gsl_permutation_alloc(invec.size());
    539     int status=gsl_sort_vector_index(p,invec.gsl_vector_p());
    540     if (status) {
    541       gsl_permutation_free(p);
    542       throw utility::GSL_error(std::string("sort_index(vector&,const vector&)",status));     
    543     }
    544     sort_index=std::vector<size_t>(p->data,p->data+p->size);
    545     gsl_permutation_free(p);
    546   }
    547 
    548 
    549   void sort_smallest_index(std::vector<size_t>& sort_index, size_t k,
    550                             const vector& invec)
    551   {
    552     assert(invec.gsl_vector_p());
    553     assert(k<=invec.size());
    554     sort_index.resize(k);
    555     gsl_sort_vector_smallest_index(&sort_index[0],k,invec.gsl_vector_p());
    556   }
    557  
    558   void sort_largest_index(std::vector<size_t>& sort_index, size_t k,
    559                             const vector& invec)
    560   {
    561     assert(invec.gsl_vector_p());
    562     assert(k<=invec.size());
    563     sort_index.resize(k);
    564     gsl_sort_vector_largest_index(&sort_index[0],k,invec.gsl_vector_p());
    565   }
    566 
    567 
    568   double sum(const vector& v)
    569   {
    570     double sum = 0;
    571     size_t vsize=v.size();
    572     for (size_t i=0; i<vsize; ++i)
    573       sum += v[i];
    574     return sum;
    575   }
    576203
    577204
     
    585212
    586213
    587   std::ostream& operator<<(std::ostream& s, const vector& a)
    588   {
    589     s.setf(std::ios::dec);
    590     s.precision(12);
    591     for (size_t j = 0; j < a.size(); ++j) {
    592       s << a[j];
    593       if ( (j+1)<a.size() )
    594         s << s.fill();
    595     }
    596     return s;
    597   }
    598 
    599214}}} // of namespace utility, yat, and thep
  • branches/peter-dev/yat/utility/vector.h

    r966 r995  
    2929*/
    3030
     31#include "vectorMutable.h"
    3132#include "Exception.h"
    3233#include "Iterator.h"
     
    8485  */
    8586
    86   class vector
     87  class vector : public vectorMutable
    8788  {
    8889  public:
    89 
    90     /// \brief vector::iterator
    91     typedef Iterator<double&, vector> iterator;
    92     /// \brief vector::const_iterator
    93     typedef Iterator<const double, const vector> const_iterator;
    9490
    9591    /**
     
    118114
    119115    /**
    120        \brief Vector view constructor.
     116       \brief The copy constructor.
    121117
    122        Create a view of vector \a v, with starting index \a offset,
    123        size \a n, and an optional \a stride.
     118       \note If the object to be copied is a vector view, the values
     119       of the view will be copied, i.e. the view is not copied.
    124120
    125        A vector view can be used as any vector with the difference
    126        that changes made to the view will also change the object that
    127        is viewed. Also, using the copy constructor will create a new
    128        vector object that is a copy of whatever is viewed. If a copy
    129        of the view is needed then you should use this constructor to
    130        obtain a copy.
    131 
    132        \note If the object viewed by the view goes out of scope or is
    133        deleted, the view becomes invalid and the result of further use
    134        is undefined.
    135 
    136        \throw GSL_error if a view cannot be set up.
     121       \throw A GSL_error is indirectly thrown if memory allocation
     122       fails.
    137123    */
    138     vector(vector& v, size_t offset, size_t n, size_t stride=1);
    139 
    140     /**
    141        \brief Vector const view constructor.
    142 
    143        Create a view of vector \a v, with starting index \a offset,
    144        size \a n, and an optional \a stride.
    145 
    146        A vector view can be used as any const vector. Using the copy
    147        constructor will create a new vector object that is a copy of
    148        whatever is viewed. If a copy of the view is needed then you
    149        should use this constructor to obtain a copy.
    150 
    151        \note If the object viewed by the view goes out of scope or is
    152        deleted, the view becomes invalid and the result of further use
    153        is undefined.
    154 
    155        \throw GSL_error if a view cannot be set up.
    156     */
    157     vector(const vector& v, size_t offset, size_t n, size_t stride=1);
    158 
    159     ///
    160     /// Matrix row/column view constructor.
    161     ///
    162     /// Create a row/column vector view of matrix \a m, pointing at
    163     /// row/column \a i. The parameter \a row is used to set whether
    164     /// the view should be a row or column view. If \a row is set to
    165     /// true, the view will be a row view (default behaviour), and,
    166     /// naturally, a column view otherwise.
    167     ///
    168     /// A vector view can be used as any vector with the difference
    169     /// that changes made to the view will also change the object that
    170     /// is viewed. Also, using the copy constructor will create a new
    171     /// vector object that is a copy of whatever is viewed. If a copy
    172     /// of the view is needed then you should use the vector view
    173     /// constructor to obtain a copy.
    174     ///
    175     /// @note If the object viewed by the view goes out of scope or is
    176     /// deleted, the view becomes invalid and the result of further
    177     /// use is undefined.
    178     ///
    179     vector(matrix& m, size_t i, bool row=true);
    180 
    181     ///
    182     /// Matrix row/column const view constructor.
    183     ///
    184     /// Create a row/column vector view of matrix \a m, pointing at
    185     /// row/column \a i. The parameter \a row is used to set whether
    186     /// the view should be a row or column view. If \a row is set to
    187     /// true, the view will be a row view (default behaviour), and,
    188     /// naturally, a column view otherwise.
    189     ///
    190     /// A const vector view can be used as any const vector. Using the
    191     /// copy constructor will create a new vector object that is a
    192     /// copy of whatever is viewed. If a copy of the view is needed
    193     /// then you should use the vector view constructor to obtain a
    194     /// copy.
    195     ///
    196     /// @note If the object viewed by the view goes out of scope or is
    197     /// deleted, the view becomes invalid and the result of further
    198     /// use is undefined.
    199     ///
    200     vector(const matrix& m, size_t i, bool row=true);
     124    vector(const vectorBase& other);
    201125
    202126    /**
     
    221145
    222146    /**
    223        \return mutable iterator to start of vector
    224      */
    225     iterator begin(void);
    226 
    227     /**
    228        \return read-only iterator to start of vector
    229      */
    230     const_iterator begin(void) const;
    231 
    232     /**
    233        \brief Make a copy of \a other.
    234 
    235        This function will make a deep copy of \a other. Memory is
    236        resized and view state is changed if needed.
    237     */
    238     const vector& clone(const vector& other);
    239 
    240     /**
    241        \brief This function performs element-wise division, \f$ this_i =
    242        this_i/other_i \; \forall i \f$.
    243 
    244        \throw GSL_error if dimensions mis-match.
    245     */
    246     void div(const vector& other);
    247 
    248     /**
    249        \return mutable iterator to end of vector
    250      */
    251     iterator end(void);
    252 
    253     /**
    254        \return read-only iterator to end of vector
    255      */
    256     const_iterator end(void) const;
    257 
    258     /**
    259        \brief Check whether vectors are equal within a user defined
    260        precision, set by \a precision.
    261 
    262        \return True if each element deviates less or equal than \a
    263        d. If any vector contain a NaN, false is always returned.
    264 
    265        \see operator== and operator!=
    266     */
    267     bool equal(const vector&, const double precision=0) const;
    268 
    269     ///
    270     /// @return A const pointer to the internal GSL vector,
    271     ///
    272     const gsl_vector* gsl_vector_p(void) const;
    273 
    274     ///
    275     /// @return A pointer to the internal GSL vector,
    276     ///
    277     gsl_vector* gsl_vector_p(void);
    278 
    279     ///
    280     /// Check if the vector object is a view (sub-vector) to another
    281     /// vector.
    282     ///
    283     /// @return True if the object is a view, false othwerwise.
    284     ///
    285     bool isview(void) const;
    286 
    287     /**
    288        \brief This function performs element-wise multiplication, \f$
    289        this_i = this_i * other_i \; \forall i \f$.
    290 
    291        \throw GSL_error if dimensions mis-match.
    292     */
    293     void mul(const vector& other);
    294 
    295     /**
    296147       @brief Resize vector
    297148       
     
    301152       vector becomes invalid.
    302153    */
    303     void resize(size_t, double init_value=0);
    304 
    305     /**
    306        \brief Reverse the order of elements in the vector.
    307     */
    308     void reverse(void);
    309 
    310     ///
    311     /// Set all elements to \a value.
    312     ///
    313     void all(const double& value);
    314 
    315     ///
    316     /// @return the number of elements in the vector.
    317     ///
    318     size_t size(void) const;
    319 
    320     /**
    321        \brief Exchange elements \a i and \a j.
    322 
    323        \throw GSL_error if vector lengths differs.
    324     */
    325     void swap(size_t i, size_t j);
    326 
    327     /**
    328        \brief Element access operator.
    329 
    330        \return Reference to element \a i.
    331 
    332        \throw If GSL range checks are enabled in the underlying GSL
    333        library a GSL_error exception is thrown if either index is out
    334        of range.
    335     */
    336     double& operator()(size_t i);
    337 
    338     /**
    339        \brief Element access operator.
    340 
    341        \return Const reference to element \a i.
    342 
    343        \throw If GSL range checks are enabled in the underlying GSL
    344        library a GSL_error exception is thrown if either index is out
    345        of range.
    346     */
    347     const double& operator()(size_t i) const;
    348 
    349     ///
    350     /// Element access operator.
    351     ///
    352     /// @return Reference to element \a i.
    353     ///
    354     double& operator[](size_t i);
    355 
    356     ///
    357     /// Const element access operator.
    358     ///
    359     /// @return The value of element \a i.
    360     ///
    361     const double& operator[](size_t i) const;
    362 
    363     /**
    364        \brief Comparison operator. Takes linear time.
    365 
    366        Checks are performed with exact matching, i.e., rounding off
    367        effects may destroy comparison. Use the equal function for
    368        comparing elements within a user defined precision.
    369 
    370        \return True if all elements are equal otherwise false.
    371 
    372        \see equal(const vector&, const double precision=0)
    373     */
    374     bool operator==(const vector&) const;
    375 
    376     /**
    377        \brief Comparison operator. Takes linear time.
    378 
    379        Checks are performed with exact matching, i.e., rounding off
    380        effects may destroy comparison. Use the equal function for
    381        comparing elements within a user defined precision.
    382 
    383        \return False if all elements are equal otherwise true.
    384 
    385        \see equal(const vector&, const double precision=0)
    386     */
    387     bool operator!=(const vector&) const;
    388 
    389     ///
    390     /// @return The dot product.
    391     ///
    392     double operator*(const vector&) const;
     154    void resize(size_t n, double init_value);
    393155
    394156    /**
    395157       \brief The assignment operator.
    396158
    397        Dimensions of the vectors must match. If the LHS vector is a
    398        view, the underlying data will be changed.
    399 
    400        \return A const reference to the resulting vector.
    401 
    402        \see void set(const vector&).
    403 
    404        \throw GSL_error if dimensions mis-match.
    405     */
    406     const vector& operator=(const vector&);
    407 
    408     /**
    409        \brief Addition and assign operator. Vector addition, \f$
    410        this_i = this_i + other_i \; \forall i \f$.
    411 
    412        \return A const reference to the resulting vector.
    413 
    414        \throw GSL_error if dimensions mis-match.
    415     */
    416     const vector& operator+=(const vector&);
    417 
    418     /**
    419        \brief Add a constant to a vector, \f$ this_i = this_i + d \;
    420        \forall i \f$.
    421 
    422159       \return A const reference to the resulting vector.
    423160    */
    424     const vector& operator+=(double d);
    425 
    426     /**
    427        \brief Subtract and assign operator. Vector subtraction, \f$
    428        this_i = this_i - other_i \; \forall i \f$.
    429 
    430        \return A const reference to the resulting vector.
    431 
    432        \throw GSL_error if dimensions mis-match.
    433     */
    434     const vector& operator-=(const vector&);
    435 
    436     /**
    437        \brief Subtract a constant to a vector, \f$ this_i = this_i - d
    438        \; \forall i \f$.
    439 
    440        \return A const reference to the resulting vector.
    441     */
    442     const vector& operator-=(double d);
    443 
    444     /**
    445        \brief Multiply with scalar and assign operator, \f$ this_i =
    446        this_i * d \; \forall i \f$.
    447 
    448        \return A const reference to the resulting vector.
    449     */
    450     const vector& operator*=(double d);
    451 
     161    const vectorMutable& operator=(vectorMutable&);
     162    const vector& operator=(const vectorBase&);
    452163
    453164  private:
     
    464175       copy, or if dimensions mis-match.
    465176    */
    466     gsl_vector* create_gsl_vector_copy(void) const;
     177    gsl_vector* create_gsl_vector_copy(const vectorBase&) const;
    467178
    468     /**
    469        \brief Clear all dynamically allocated memory.
    470 
    471        Internal utility function.
    472     */
    473     void delete_allocated_memory(void);
    474 
    475     gsl_vector* v_;
    476     gsl_vector_view* view_;
    477     const gsl_vector_const_view* view_const_;
    478     // proxy_v_ is used to access the proper underlying gsl_vector
    479     // in all const member functions. It is not used by design for
    480     // non-const vector functions and operators. This is to make sure
    481     // that runtime errors occur if a const vector is used in an
    482     // inappropriate manner such as on left hand side in assignment
    483     // (remember, v_ is null for const vector views).
    484     const gsl_vector* proxy_v_;
    485179  };
    486 
    487   /**
    488      \brief Check if all elements of the vector are zero.
    489 
    490      \return True if all elements in the vector is zero, false
    491      othwerwise.
    492   */
    493   bool isnull(const vector&);
    494 
    495   /**
    496      \brief Get the maximum value of the vector.
    497 
    498      \return The maximum value of the vector.
    499   */
    500   double max(const vector&);
    501 
    502   /**
    503      \brief Locate the maximum value in the vector.
    504 
    505      \return The index to the maximum value of the vector.
    506 
    507      \note Lower index has precedence.
    508   */
    509   size_t max_index(const vector&);
    510 
    511   /**
    512      \brief Get the minimum value of the vector.
    513 
    514      \return The minimum value of the vector.
    515   */
    516   double min(const vector&);
    517 
    518   /**
    519      \brief Locate the minimum value in the vector.
    520 
    521      \return The index to the minimum value of the vector.
    522 
    523      \note Lower index has precedence.
    524   */
    525   size_t min_index(const vector&);
    526 
    527   /**
    528      \brief Create a vector \a flag indicating NaN's in another vector
    529      \a templat.
    530 
    531      The \a flag vector is changed to contain 1's and 0's only. A 1
    532      means that the corresponding element in the \a templat vector is
    533      valid and a zero means that the corresponding element is a NaN.
    534 
    535      \note Space for vector \a flag is reallocated to fit the size of
    536      vector \a templat if sizes mismatch.
    537 
    538      \return True if the \a templat vector contains at least one NaN.
    539   */
    540   bool nan(const vector& templat, vector& flag);
    541 
    542   /**
    543      \brief Transforms a vector to a basis vector.
    544 
    545      All elements are set to zero except the \a i-th element which is
    546      set to one.
    547   */
    548   void set_basis(vector&, size_t i);
    549 
    550   /**
    551      Randomly shuffles the elements in vector \a invec
    552   */
    553   void shuffle(vector& invec);
    554 
    555   /**
    556      Sort the elements in the vector.
    557   */
    558   void sort(vector&);
    559 
    560   /**
    561      Create a vector \a sort_index containing the indeces of
    562      elements in a another vector \a invec.  The elements of \a
    563      sort_index give the index of the vector element which would
    564      have been stored in that position if the vector had been sorted
    565      in place. The first element of \a sort_index gives the index of the least
    566      element in \a invec, and the last element of \a sort_index gives the index of the
    567      greatest element in \a invec . The vector \a invec is not changed.
    568 
    569   */
    570   void sort_index(std::vector<size_t>& sort_index, const vector& invec);
    571 
    572   /** Similar to sort_index but creates a vector with indices to the \a k
    573   smallest elements in \a invec. 
    574   */
    575   void sort_smallest_index(std::vector<size_t>& sort_index, size_t k, const
    576   vector& invec);
    577 
    578   /** Similar to sort_index but creates a vector with indices to the \a k
    579   largest elements in \a invec. 
    580   */
    581   void sort_largest_index(std::vector<size_t>& sort_index, size_t k, const
    582   vector& invec);
    583 
    584  
    585 
    586   /**
    587      \brief Calculate the sum of all vector elements.
    588 
    589      \return The sum.
    590   */
    591   double sum(const vector&);
    592180
    593181  /**
     
    600188  void swap(vector&, vector&);
    601189
    602   /**
    603      \brief The output operator for the vector class.
    604   */
    605   std::ostream& operator<<(std::ostream&, const vector&);
    606 
    607190}}} // of namespace utility, yat, and theplu
    608191
  • branches/peter-dev/yat/utility/vectorBase.cc

    r994 r995  
    2525*/
    2626
    27 #include "vector.h"
     27#include "vectorBase.h"
    2828#include "matrix.h"
    2929#include "utility.h"
     
    4343
    4444
    45   vector::vector(void)
    46     : v_(NULL), view_(NULL), view_const_(NULL), proxy_v_(NULL)
    47   {
    48   }
    49 
    50 
    51   vector::vector(size_t n, double init_value)
    52     : v_(gsl_vector_alloc(n)), view_(NULL), view_const_(NULL),
    53       proxy_v_(v_)
    54   {
    55     if (!v_)
    56       throw utility::GSL_error("vector::vector failed to allocate memory");
    57 
    58     all(init_value);
    59   }
    60 
    61 
    62   vector::vector(const vector& other)
    63     : v_(other.create_gsl_vector_copy()), view_(NULL),
    64       view_const_(NULL), proxy_v_(v_)
    65   {
    66   }
    67 
    68 
    69   vector::vector(vector& v, size_t offset, size_t n, size_t stride)
    70     : view_const_(NULL)
    71   {
    72     view_ = new gsl_vector_view(gsl_vector_subvector_with_stride(v.v_,offset,
    73                                                                  stride,n));
    74     if (!view_)
    75       throw utility::GSL_error("vector::vector failed to setup view");
    76     proxy_v_ = v_ = &(view_->vector);
    77   }
    78 
    79 
    80   vector::vector(const vector& v, size_t offset, size_t n, size_t stride)
    81     : v_(NULL), view_(NULL)
    82   {
    83     view_const_ = new gsl_vector_const_view(
    84                    gsl_vector_const_subvector_with_stride(v.v_,offset,stride,n));
    85     if (!view_const_)
    86       throw utility::GSL_error("vector::vector failed to setup view");
    87     proxy_v_ = &(view_const_->vector);
    88   }
    89 
    90 
    91   vector::vector(matrix& m, size_t i, bool row)
    92     : view_const_(NULL)
    93   {
    94     view_=new gsl_vector_view(row ?
    95                               gsl_matrix_row   (m.gsl_matrix_p(),i) :
    96                               gsl_matrix_column(m.gsl_matrix_p(),i)  );
    97     if (!view_)
    98       throw utility::GSL_error("vector::vector failed to setup view");
    99     proxy_v_ = v_ = &(view_->vector);
    100   }
    101 
    102 
    103   vector::vector(const matrix& m, size_t i, bool row)
    104     : v_(NULL), view_(NULL)
    105   {
    106     view_const_ = new gsl_vector_const_view(row ?
    107                                    gsl_matrix_const_row   (m.gsl_matrix_p(),i) :
    108                                    gsl_matrix_const_column(m.gsl_matrix_p(),i) );
    109     if (!view_const_)
    110       throw utility::GSL_error("vector::vector failed to setup view");
    111     proxy_v_ = &(view_const_->vector);
    112   }
    113 
    114 
    115   vector::vector(std::istream& is, char sep)
    116     throw (utility::IO_error, std::exception)
    117     : view_(NULL), view_const_(NULL)
    118   {
    119     // read the data file and store in stl vectors (dynamically
    120     // expandable)
    121     std::vector<std::vector<double> > data_matrix;
    122     u_int nof_columns=0;
    123     u_int nof_rows=0;
    124     std::string line;
    125     while(getline(is, line, '\n')) {
    126       // Empty lines
    127       if (!line.size())
    128           continue;
    129       nof_rows++;
    130 
    131       std::vector<double> v;
    132       std::string element;
    133       std::stringstream ss(line);
    134       bool ok=true;
    135       while(ok) {
    136         if(sep=='\0')
    137           ok=(ss>>element);
    138         else
    139           ok=getline(ss, element, sep);
    140         if(!ok)
    141           break;       
    142 
    143         if(utility::is_double(element)) {
    144           v.push_back(atof(element.c_str()));
    145         }
    146         else if (!element.size() || utility::is_nan(element)) {
    147           v.push_back(std::numeric_limits<double>::quiet_NaN());
    148         }
    149         else {
    150           std::stringstream ss("Warning: '");
    151           ss << element << "' is not a double.";
    152           throw IO_error(ss.str());
    153         }
    154       }
    155       if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
    156           v.push_back(std::numeric_limits<double>::quiet_NaN());
    157       if (!nof_columns)
    158         nof_columns=v.size();
    159       else if (nof_rows && (nof_columns>1)) {
    160         std::ostringstream s;
    161         s << "vector::vector(std::istream&) data file error:\n"
    162           << "    File has inconsistent number of rows (" << nof_rows
    163           << ") and columns (" << nof_columns
    164           << ").\n    Expected a row or a column vector.";
    165         throw utility::IO_error(s.str());
    166       }
    167       else if (v.size()!=nof_columns) {
    168         std::ostringstream s;
    169         s << "vector::vector(std::istream&) data file error:\n"
    170           << "    Line " << nof_rows << " has " << v.size()
    171           << " columns; expected " << nof_columns << " column.";
    172         throw utility::IO_error(s.str());
    173       }
    174       data_matrix.push_back(v);
    175     }
    176 
    177     // manipulate the state of the stream to be good
    178     is.clear(std::ios::goodbit);
    179     // convert the data to a gsl vector
    180     proxy_v_ = v_ = gsl_vector_alloc(nof_rows*nof_columns);
    181     if (!v_)
    182       throw utility::GSL_error("vector::vector failed to allocate memory");
    183     size_t n=0;
    184     // if gsl error handler disabled, out of bounds index will not
    185     // abort the program.
    186     for (size_t i=0; i<nof_rows; i++)
    187       for (size_t j=0; j<nof_columns; j++)
    188         gsl_vector_set( v_, n++, data_matrix[i][j] );
    189   }
    190 
    191 
    192   vector::~vector(void)
    193   {
    194     delete_allocated_memory();
    195   }
    196 
    197 
    198   vector::iterator vector::begin(void)
    199   {
    200     return vector::iterator(*this, 0);
    201   }
    202 
    203 
    204   vector::const_iterator vector::begin(void) const
    205   {
    206     return vector::const_iterator(*this, 0);
    207   }
    208 
    209 
    210   const vector& vector::clone(const vector& other)
    211   {
    212     if (this!=&other) {
    213 
    214       delete_allocated_memory();
    215 
    216       if (other.view_) {
    217         view_ = new gsl_vector_view(*other.view_);
    218         proxy_v_ = v_ = &(view_->vector);
    219       }
    220       else if (other.view_const_) {
    221         view_const_ = new gsl_vector_const_view(*other.view_const_);
    222         proxy_v_ = &(view_const_->vector);
    223         v_=NULL;
    224       }
    225       else if (other.v_)
    226         proxy_v_ = v_ = other.create_gsl_vector_copy();
    227 
    228     }
    229     return *this;
    230   }
    231 
    232 
    233   gsl_vector* vector::create_gsl_vector_copy(void) const
    234   {
    235     gsl_vector* vec = gsl_vector_alloc(size());
    236     if (!vec)
    237       throw utility::GSL_error("vector::create_gsl_vector_copy failed to allocate memory");
    238     if (gsl_vector_memcpy(vec, proxy_v_))
    239       throw utility::GSL_error("vector::create_gsl_matrix_copy dimension mis-match");
    240     return vec;
    241   }
    242 
    243 
    244   void vector::delete_allocated_memory(void)
    245   {
    246     if (view_)
    247       delete view_;
    248     else if (view_const_)
    249       delete view_const_;
    250     else if (v_)
    251       gsl_vector_free(v_);
    252     proxy_v_=v_=NULL;
    253   }
    254 
    255 
    256   void vector::div(const vector& other)
    257   {
    258     assert(v_);
    259     int status=gsl_vector_div(v_,other.gsl_vector_p());
    260     if (status)
    261       throw utility::GSL_error(std::string("vector::div",status));
    262   }
    263 
    264 
    265   vector::iterator vector::end(void)
    266   {
    267     return vector::iterator(*this, size());
    268   }
    269 
    270 
    271   vector::const_iterator vector::end(void) const
    272   {
    273     return vector::const_iterator(*this, size());
    274   }
    275 
    276 
    277   bool vector::equal(const vector& other, const double d) const
     45  vectorBase::vectorBase(const gsl_vector* v)
     46    : const_vec_(v)
     47  {
     48  }
     49
     50
     51  vectorBase::~vectorBase(void)
     52  {
     53  }
     54
     55
     56  vectorBase::const_iterator vectorBase::begin(void) const
     57  {
     58    return vectorBase::const_iterator(*this, 0);
     59  }
     60
     61
     62  vectorBase::const_iterator vectorBase::end(void) const
     63  {
     64    return vectorBase::const_iterator(*this, size());
     65  }
     66
     67
     68  bool vectorBase::equal(const vectorBase& other, const double d) const
    27869  {
    27970    if (this==&other)
     
    28475    // abort the program.
    28576    for (size_t i=0; i<size(); ++i)
    286       // The two last condition checks are needed for NaN detection
    28777      if (fabs( (*this)(i)-other(i) ) > d ||
    288           (*this)(i)!=(*this)(i) || other(i)!=other(i))
     78          std::isnan((*this)(i)) || std::isnan(other(i)) )
    28979        return false;
    29080    return true;
     
    29282
    29383
    294   const gsl_vector* vector::gsl_vector_p(void) const
    295   {
    296     return proxy_v_;
    297   }
    298 
    299 
    300   gsl_vector* vector::gsl_vector_p(void)
    301   {
    302     return v_;
    303   }
    304 
    305 
    306   bool vector::isview(void) const
    307   {
    308     return view_ || view_const_;
    309   }
    310 
    311 
    312   void vector::mul(const vector& other)
    313   {
    314     assert(v_);
    315     int status=gsl_vector_mul(v_,other.gsl_vector_p());
    316     if (status)
    317       throw utility::GSL_error(std::string("vector::mul",status));
    318   }
    319 
    320 
    321   void vector::resize(size_t n, double init_value)
    322   {
    323     delete_allocated_memory();
    324     proxy_v_ = v_ = gsl_vector_alloc(n);
    325     if (!v_)
    326       throw utility::GSL_error("vector::vector failed to allocate memory");
    327     all(init_value);
    328   }
    329 
    330 
    331   void vector::reverse(void)
    332   {
    333     assert(v_);
    334     gsl_vector_reverse(v_);
    335   }
    336 
    337 
    338   void vector::all(const double& value)
    339   {
    340     assert(v_);
    341     gsl_vector_set_all(v_,value);
    342   }
    343 
    344 
    345   size_t vector::size(void) const
    346   {
    347     if (!proxy_v_)
     84  const gsl_vector* vectorBase::gsl_vector_p(void) const
     85  {
     86    return const_vec_;
     87  }
     88
     89
     90  size_t vectorBase::size(void) const
     91  {
     92    if (const_vec_)
    34893      return 0;
    349     return proxy_v_->size;
    350   }
    351 
    352 
    353   void vector::swap(size_t i, size_t j)
    354   {
    355     assert(v_);
    356     int status=gsl_vector_swap_elements(v_, i, j);
    357     if (status)
    358       throw utility::GSL_error(std::string("vector::swap_elements",status));
    359   }
    360 
    361 
    362   double& vector::operator()(size_t i)
    363   {
    364     assert(v_);
    365     double* d=gsl_vector_ptr(v_, i);
     94    return const_vec_->size;
     95  }
     96
     97
     98  const double& vectorBase::operator()(size_t i) const
     99  {
     100    const double* d=gsl_vector_const_ptr(const_vec_, i);
    366101    if (!d)
    367       throw utility::GSL_error("vector::operator()",GSL_EINVAL);
     102      throw utility::GSL_error("vectorBase::operator()",GSL_EINVAL);
    368103    return *d;
    369104  }
    370105
    371106
    372   const double& vector::operator()(size_t i) const
    373   {
    374     const double* d=gsl_vector_const_ptr(proxy_v_, i);
     107  const double& vectorBase::operator[](size_t i) const
     108  {
     109    const double* d=gsl_vector_const_ptr(const_vec_, i);
    375110    if (!d)
    376       throw utility::GSL_error("vector::operator()",GSL_EINVAL);
     111      throw utility::GSL_error("vectorBase::operator[]",GSL_EINVAL);
    377112    return *d;
    378113  }
    379114
    380115
    381   double& vector::operator[](size_t i)
    382   {
    383     return this->operator()(i);
    384   }
    385 
    386 
    387   const double& vector::operator[](size_t i) const
    388   {
    389     return this->operator()(i);
    390   }
    391 
    392 
    393   bool vector::operator==(const vector& other) const
     116  bool vectorBase::operator==(const vectorBase& other) const
    394117  {
    395118    return equal(other);
     
    397120
    398121
    399   bool vector::operator!=(const vector& other) const
     122  bool vectorBase::operator!=(const vectorBase& other) const
    400123  {
    401124    return !equal(other);
     
    403126
    404127
    405   double vector::operator*( const vector &other ) const
     128  double vectorBase::operator*( const vectorBase &other ) const
    406129  {
    407130    double res = 0.0;;
     
    412135
    413136
    414   const vector& vector::operator=( const vector& other )
    415   {
    416     assert(v_);
    417     if (gsl_vector_memcpy(v_,other.gsl_vector_p()))
    418       throw utility::GSL_error("vector::set dimension mis-match");
    419     return *this;
    420   }
    421 
    422 
    423   const vector& vector::operator+=(const vector& other)
    424   {
    425     assert(v_);
    426     int status=gsl_vector_add(v_, other.gsl_vector_p());
    427     if (status)
    428       throw utility::GSL_error(std::string("vector::add", status));
    429     return *this;
    430   }
    431 
    432 
    433   const vector& vector::operator+=(double d)
    434   {
    435     assert(v_);
    436     gsl_vector_add_constant(v_, d);
    437     return *this;
    438   }
    439 
    440 
    441   const vector& vector::operator-=(const vector& other)
    442   {
    443     assert(v_);
    444     int status=gsl_vector_sub(v_, other.gsl_vector_p());
    445     if (status)
    446       throw utility::GSL_error(std::string("vector::sub", status));
    447     return *this;
    448   }
    449 
    450 
    451   const vector& vector::operator-=(const double d)
    452   {
    453     assert(v_);
    454     gsl_vector_add_constant(v_, -d);
    455     return *this;
    456   }
    457 
    458 
    459   const vector& vector::operator*=(const double d)
    460   {
    461     assert(v_);
    462     gsl_vector_scale(v_, d);
    463     return *this;
    464   }
    465 
    466 
    467   bool isnull(const vector& v)
     137  bool isnull(const vectorBase& v)
    468138  {
    469139    return gsl_vector_isnull(v.gsl_vector_p());
     
    471141
    472142
    473   double max(const vector& v)
     143  double max(const vectorBase& v)
    474144  {
    475145    return gsl_vector_max(v.gsl_vector_p());
     
    477147
    478148
    479   size_t max_index(const vector& v)
     149  size_t max_index(const vectorBase& v)
    480150  {
    481151    return gsl_vector_max_index(v.gsl_vector_p());
     
    483153
    484154
    485   double min(const vector& v)
     155  double min(const vectorBase& v)
    486156  {
    487157    return gsl_vector_min(v.gsl_vector_p());
     
    489159
    490160
    491   size_t min_index(const vector& v)
     161  size_t min_index(const vectorBase& v)
    492162  {
    493163    return gsl_vector_min_index(v.gsl_vector_p());
     
    495165
    496166
    497   bool nan(const vector& templat, vector& flag)
    498   {
    499     size_t vsize=templat.size();
    500     if (vsize!=flag.size())
    501       flag.clone(vector(vsize,1.0));
    502     else
    503       flag.all(1.0);
     167  bool nan(const vectorBase& templat, vector& flag)
     168  {
     169    size_t vsize(templat.size());
     170    flag = vector(vsize, 1.0);
    504171    bool nan=false;
    505172    for (size_t i=0; i<vsize; i++)
     
    512179
    513180
    514   void set_basis(vector& v, size_t i)
    515   {
    516     assert(v.gsl_vector_p());
    517     gsl_vector_set_basis(v.gsl_vector_p(),i);
    518   }
    519 
    520 
    521   void shuffle(vector& invec)
    522   {
    523     random::DiscreteUniform rnd;
    524     std::random_shuffle(invec.begin(), invec.end(), rnd);
    525   }
    526 
    527 
    528   void sort(vector& v)
    529   {
    530     assert(v.gsl_vector_p());
    531     std::sort(v.begin(), v.end());
    532   }
    533 
    534 
    535   void sort_index(std::vector<size_t>& sort_index, const vector& invec)
     181  void sort_index(std::vector<size_t>& sort_index, const vectorBase& invec)
    536182  {
    537183    assert(invec.gsl_vector_p());
     
    540186    if (status) {
    541187      gsl_permutation_free(p);
    542       throw utility::GSL_error(std::string("sort_index(vector&,const vector&)",status));     
     188      throw utility::GSL_error(std::string("sort_index(vector&,const vectorBase&)",status));     
    543189    }
    544190    sort_index=std::vector<size_t>(p->data,p->data+p->size);
     
    548194
    549195  void sort_smallest_index(std::vector<size_t>& sort_index, size_t k,
    550                             const vector& invec)
     196                            const vectorBase& invec)
    551197  {
    552198    assert(invec.gsl_vector_p());
     
    557203 
    558204  void sort_largest_index(std::vector<size_t>& sort_index, size_t k,
    559                             const vector& invec)
     205                            const vectorBase& invec)
    560206  {
    561207    assert(invec.gsl_vector_p());
     
    566212
    567213
    568   double sum(const vector& v)
     214  double sum(const vectorBase& v)
    569215  {
    570216    double sum = 0;
    571217    size_t vsize=v.size();
    572218    for (size_t i=0; i<vsize; ++i)
    573       sum += v[i];
     219      sum += v(i);
    574220    return sum;
    575221  }
    576222
    577223
    578   void swap(vector& v, vector& w)
    579   {
    580     assert(v.gsl_vector_p()); assert(w.gsl_vector_p());
    581     int status=gsl_vector_swap(v.gsl_vector_p(),w.gsl_vector_p());
    582     if (status)
    583       throw utility::GSL_error(std::string("swap(vector&,vector&)",status));
    584   }
    585 
    586 
    587   std::ostream& operator<<(std::ostream& s, const vector& a)
     224  std::ostream& operator<<(std::ostream& s, const vectorBase& a)
    588225  {
    589226    s.setf(std::ios::dec);
    590227    s.precision(12);
    591228    for (size_t j = 0; j < a.size(); ++j) {
    592       s << a[j];
     229      s << a(j);
    593230      if ( (j+1)<a.size() )
    594231        s << s.fill();
  • branches/peter-dev/yat/utility/vectorBase.h

    r994 r995  
    1 #ifndef _theplu_yat_utility_vector_
    2 #define _theplu_yat_utility_vector_
     1#ifndef _theplu_yat_utility_vector_base_
     2#define _theplu_yat_utility_vector_base_
    33
    44// $Id$
     
    4444
    4545  class matrix;
     46  class vector;
    4647
    4748  /**
     
    8485  */
    8586
    86   class vector
     87  class vectorBase
    8788  {
    8889  public:
    8990
    90     /// \brief vector::iterator
    91     typedef Iterator<double&, vector> iterator;
    92     /// \brief vector::const_iterator
    93     typedef Iterator<const double, const vector> const_iterator;
    94 
    95     /**
    96        \brief The default constructor.
    97     */
    98     vector(void);
    99 
    100     /**
    101        \brief Allocates memory space for \a n elements, and sets all
    102        elements to \a init_value.
    103        
    104        \throw GSL_error if memory allocation fails.
    105     */
    106     vector(size_t n, double init_value=0);
    107 
    108     /**
    109        \brief The copy constructor.
    110 
    111        \note If the object to be copied is a vector view, the values
    112        of the view will be copied, i.e. the view is not copied.
    113 
    114        \throw A GSL_error is indirectly thrown if memory allocation
    115        fails.
    116     */
    117     vector(const vector& other);
    118 
    119     /**
    120        \brief Vector view constructor.
    121 
    122        Create a view of vector \a v, with starting index \a offset,
    123        size \a n, and an optional \a stride.
    124 
    125        A vector view can be used as any vector with the difference
    126        that changes made to the view will also change the object that
    127        is viewed. Also, using the copy constructor will create a new
    128        vector object that is a copy of whatever is viewed. If a copy
    129        of the view is needed then you should use this constructor to
    130        obtain a copy.
    131 
    132        \note If the object viewed by the view goes out of scope or is
    133        deleted, the view becomes invalid and the result of further use
    134        is undefined.
    135 
    136        \throw GSL_error if a view cannot be set up.
    137     */
    138     vector(vector& v, size_t offset, size_t n, size_t stride=1);
    139 
    140     /**
    141        \brief Vector const view constructor.
    142 
    143        Create a view of vector \a v, with starting index \a offset,
    144        size \a n, and an optional \a stride.
    145 
    146        A vector view can be used as any const vector. Using the copy
    147        constructor will create a new vector object that is a copy of
    148        whatever is viewed. If a copy of the view is needed then you
    149        should use this constructor to obtain a copy.
    150 
    151        \note If the object viewed by the view goes out of scope or is
    152        deleted, the view becomes invalid and the result of further use
    153        is undefined.
    154 
    155        \throw GSL_error if a view cannot be set up.
    156     */
    157     vector(const vector& v, size_t offset, size_t n, size_t stride=1);
    158 
    159     ///
    160     /// Matrix row/column view constructor.
    161     ///
    162     /// Create a row/column vector view of matrix \a m, pointing at
    163     /// row/column \a i. The parameter \a row is used to set whether
    164     /// the view should be a row or column view. If \a row is set to
    165     /// true, the view will be a row view (default behaviour), and,
    166     /// naturally, a column view otherwise.
    167     ///
    168     /// A vector view can be used as any vector with the difference
    169     /// that changes made to the view will also change the object that
    170     /// is viewed. Also, using the copy constructor will create a new
    171     /// vector object that is a copy of whatever is viewed. If a copy
    172     /// of the view is needed then you should use the vector view
    173     /// constructor to obtain a copy.
    174     ///
    175     /// @note If the object viewed by the view goes out of scope or is
    176     /// deleted, the view becomes invalid and the result of further
    177     /// use is undefined.
    178     ///
    179     vector(matrix& m, size_t i, bool row=true);
    180 
    181     ///
    182     /// Matrix row/column const view constructor.
    183     ///
    184     /// Create a row/column vector view of matrix \a m, pointing at
    185     /// row/column \a i. The parameter \a row is used to set whether
    186     /// the view should be a row or column view. If \a row is set to
    187     /// true, the view will be a row view (default behaviour), and,
    188     /// naturally, a column view otherwise.
    189     ///
    190     /// A const vector view can be used as any const vector. Using the
    191     /// copy constructor will create a new vector object that is a
    192     /// copy of whatever is viewed. If a copy of the view is needed
    193     /// then you should use the vector view constructor to obtain a
    194     /// copy.
    195     ///
    196     /// @note If the object viewed by the view goes out of scope or is
    197     /// deleted, the view becomes invalid and the result of further
    198     /// use is undefined.
    199     ///
    200     vector(const matrix& m, size_t i, bool row=true);
    201 
    202     /**
    203        \brief The istream constructor.
    204 
    205        Either elements should be separated with white space characters
    206        (default), or elements should be separated by the delimiter \a
    207        sep. When delimiter \a sep is used empty elements are stored as
    208        NaN's (except that empty lines are ignored). The end of input
    209        to the vector is at end of file marker.
    210 
    211        \throw GSL_error if memory allocation fails, IO_error if
    212        unexpected input is found in the input stream.
    213     */
    214     explicit vector(std::istream &, char sep='\0')
    215       throw(utility::IO_error, std::exception);
     91    /// \brief vectorBase::const_iterator
     92    typedef Iterator<const double, const vectorBase> const_iterator;
     93
     94    /**
     95       \brief Constructor.
     96    */
     97    vectorBase(const gsl_vector*);
    21698
    21799    ///
    218100    /// The destructor.
    219101    ///
    220     ~vector(void);
    221 
    222     /**
    223        \return mutable iterator to start of vector
    224      */
    225     iterator begin(void);
    226 
    227     /**
    228        \return read-only iterator to start of vector
     102    virtual ~vectorBase(void);
     103
     104    /**
     105       \return read-only iterator to start of vectorBase
    229106     */
    230107    const_iterator begin(void) const;
    231108
    232109    /**
    233        \brief Make a copy of \a other.
    234 
    235        This function will make a deep copy of \a other. Memory is
    236        resized and view state is changed if needed.
    237     */
    238     const vector& clone(const vector& other);
    239 
    240     /**
    241        \brief This function performs element-wise division, \f$ this_i =
    242        this_i/other_i \; \forall i \f$.
    243 
    244        \throw GSL_error if dimensions mis-match.
    245     */
    246     void div(const vector& other);
    247 
    248     /**
    249        \return mutable iterator to end of vector
    250      */
    251     iterator end(void);
    252 
    253     /**
    254        \return read-only iterator to end of vector
     110       \return read-only iterator to end of vectorBase
    255111     */
    256112    const_iterator end(void) const;
    257113
    258114    /**
    259        \brief Check whether vectors are equal within a user defined
     115       \brief Check whether vectorBases are equal within a user defined
    260116       precision, set by \a precision.
    261117
    262118       \return True if each element deviates less or equal than \a
    263        d. If any vector contain a NaN, false is always returned.
     119       d. If any vectorBase contain a NaN, false is always returned.
    264120
    265121       \see operator== and operator!=
    266122    */
    267     bool equal(const vector&, const double precision=0) const;
     123    bool equal(const vectorBase&, const double precision=0) const;
    268124
    269125    ///
     
    273129
    274130    ///
    275     /// @return A pointer to the internal GSL vector,
    276     ///
    277     gsl_vector* gsl_vector_p(void);
    278 
    279     ///
    280     /// Check if the vector object is a view (sub-vector) to another
    281     /// vector.
    282     ///
    283     /// @return True if the object is a view, false othwerwise.
    284     ///
    285     bool isview(void) const;
    286 
    287     /**
    288        \brief This function performs element-wise multiplication, \f$
    289        this_i = this_i * other_i \; \forall i \f$.
    290 
    291        \throw GSL_error if dimensions mis-match.
    292     */
    293     void mul(const vector& other);
    294 
    295     /**
    296        @brief Resize vector
    297        
    298        All elements are set to @a init_value.
    299 
    300        \note Underlying GSL vector is destroyed and a view into this
    301        vector becomes invalid.
    302     */
    303     void resize(size_t, double init_value=0);
    304 
    305     /**
    306        \brief Reverse the order of elements in the vector.
    307     */
    308     void reverse(void);
    309 
    310     ///
    311     /// Set all elements to \a value.
    312     ///
    313     void all(const double& value);
    314 
    315     ///
    316     /// @return the number of elements in the vector.
     131    /// @return the number of elements in the vectorBase.
    317132    ///
    318133    size_t size(void) const;
    319134
    320135    /**
    321        \brief Exchange elements \a i and \a j.
    322 
    323        \throw GSL_error if vector lengths differs.
    324     */
    325     void swap(size_t i, size_t j);
    326 
    327     /**
    328136       \brief Element access operator.
    329137
    330        \return Reference to element \a i.
     138       \return Const reference to element \a i.
    331139
    332140       \throw If GSL range checks are enabled in the underlying GSL
     
    334142       of range.
    335143    */
    336     double& operator()(size_t i);
    337 
    338     /**
    339        \brief Element access operator.
    340 
    341        \return Const reference to element \a i.
    342 
    343        \throw If GSL range checks are enabled in the underlying GSL
    344        library a GSL_error exception is thrown if either index is out
    345        of range.
    346     */
    347144    const double& operator()(size_t i) const;
    348 
    349     ///
    350     /// Element access operator.
    351     ///
    352     /// @return Reference to element \a i.
    353     ///
    354     double& operator[](size_t i);
    355 
    356     ///
    357     /// Const element access operator.
    358     ///
    359     /// @return The value of element \a i.
    360     ///
    361145    const double& operator[](size_t i) const;
    362146
     
    370154       \return True if all elements are equal otherwise false.
    371155
    372        \see equal(const vector&, const double precision=0)
    373     */
    374     bool operator==(const vector&) const;
     156       \see equal(const vectorBase&, const double precision=0)
     157    */
     158    bool operator==(const vectorBase&) const;
    375159
    376160    /**
     
    383167       \return False if all elements are equal otherwise true.
    384168
    385        \see equal(const vector&, const double precision=0)
    386     */
    387     bool operator!=(const vector&) const;
     169       \see equal(const vectorBase&, const double precision=0)
     170    */
     171    bool operator!=(const vectorBase&) const;
    388172
    389173    ///
    390174    /// @return The dot product.
    391175    ///
    392     double operator*(const vector&) const;
    393 
    394     /**
    395        \brief The assignment operator.
    396 
    397        Dimensions of the vectors must match. If the LHS vector is a
    398        view, the underlying data will be changed.
    399 
    400        \return A const reference to the resulting vector.
    401 
    402        \see void set(const vector&).
    403 
    404        \throw GSL_error if dimensions mis-match.
    405     */
    406     const vector& operator=(const vector&);
    407 
    408     /**
    409        \brief Addition and assign operator. Vector addition, \f$
    410        this_i = this_i + other_i \; \forall i \f$.
    411 
    412        \return A const reference to the resulting vector.
    413 
    414        \throw GSL_error if dimensions mis-match.
    415     */
    416     const vector& operator+=(const vector&);
    417 
    418     /**
    419        \brief Add a constant to a vector, \f$ this_i = this_i + d \;
    420        \forall i \f$.
    421 
    422        \return A const reference to the resulting vector.
    423     */
    424     const vector& operator+=(double d);
    425 
    426     /**
    427        \brief Subtract and assign operator. Vector subtraction, \f$
    428        this_i = this_i - other_i \; \forall i \f$.
    429 
    430        \return A const reference to the resulting vector.
    431 
    432        \throw GSL_error if dimensions mis-match.
    433     */
    434     const vector& operator-=(const vector&);
    435 
    436     /**
    437        \brief Subtract a constant to a vector, \f$ this_i = this_i - d
    438        \; \forall i \f$.
    439 
    440        \return A const reference to the resulting vector.
    441     */
    442     const vector& operator-=(double d);
    443 
    444     /**
    445        \brief Multiply with scalar and assign operator, \f$ this_i =
    446        this_i * d \; \forall i \f$.
    447 
    448        \return A const reference to the resulting vector.
    449     */
    450     const vector& operator*=(double d);
    451 
    452 
    453   private:
    454 
    455     /**
    456        \brief Create a new copy of the internal GSL vector.
    457 
    458        Necessary memory for the new GSL vector is allocated and the
    459        caller is responsible for freeing the allocated memory.
    460 
    461        \return A pointer to a copy of the internal GSL vector.
    462 
    463        \throw GSL_error if memory cannot be allocated for the new
    464        copy, or if dimensions mis-match.
    465     */
    466     gsl_vector* create_gsl_vector_copy(void) const;
    467 
    468     /**
    469        \brief Clear all dynamically allocated memory.
    470 
    471        Internal utility function.
    472     */
    473     void delete_allocated_memory(void);
    474 
    475     gsl_vector* v_;
    476     gsl_vector_view* view_;
    477     const gsl_vector_const_view* view_const_;
    478     // proxy_v_ is used to access the proper underlying gsl_vector
    479     // in all const member functions. It is not used by design for
    480     // non-const vector functions and operators. This is to make sure
    481     // that runtime errors occur if a const vector is used in an
    482     // inappropriate manner such as on left hand side in assignment
    483     // (remember, v_ is null for const vector views).
    484     const gsl_vector* proxy_v_;
     176    double operator*(const vectorBase&) const;
     177
     178  protected:
     179    const gsl_vector* const_vec_;
    485180  };
    486181
    487182  /**
    488      \brief Check if all elements of the vector are zero.
    489 
    490      \return True if all elements in the vector is zero, false
     183     \brief Check if all elements of the vectorBase are zero.
     184
     185     \return True if all elements in the vectorBase is zero, false
    491186     othwerwise.
    492187  */
    493   bool isnull(const vector&);
    494 
    495   /**
    496      \brief Get the maximum value of the vector.
    497 
    498      \return The maximum value of the vector.
    499   */
    500   double max(const vector&);
    501 
    502   /**
    503      \brief Locate the maximum value in the vector.
    504 
    505      \return The index to the maximum value of the vector.
     188  bool isnull(const vectorBase&);
     189
     190  /**
     191     \brief Get the maximum value of the vectorBase.
     192
     193     \return The maximum value of the vectorBase.
     194  */
     195  double max(const vectorBase&);
     196
     197  /**
     198     \brief Locate the maximum value in the vectorBase.
     199
     200     \return The index to the maximum value of the vectorBase.
    506201
    507202     \note Lower index has precedence.
    508203  */
    509   size_t max_index(const vector&);
    510 
    511   /**
    512      \brief Get the minimum value of the vector.
    513 
    514      \return The minimum value of the vector.
    515   */
    516   double min(const vector&);
    517 
    518   /**
    519      \brief Locate the minimum value in the vector.
    520 
    521      \return The index to the minimum value of the vector.
     204  size_t max_index(const vectorBase&);
     205
     206  /**
     207     \brief Get the minimum value of the vectorBase.
     208
     209     \return The minimum value of the vectorBase.
     210  */
     211  double min(const vectorBase&);
     212
     213  /**
     214     \brief Locate the minimum value in the vectorBase.
     215
     216     \return The index to the minimum value of the vectorBase.
    522217
    523218     \note Lower index has precedence.
    524219  */
    525   size_t min_index(const vector&);
    526 
    527   /**
    528      \brief Create a vector \a flag indicating NaN's in another vector
     220  size_t min_index(const vectorBase&);
     221
     222  /**
     223     \brief Create a vectorBase \a flag indicating NaN's in another vectorBase
    529224     \a templat.
    530225
    531      The \a flag vector is changed to contain 1's and 0's only. A 1
    532      means that the corresponding element in the \a templat vector is
     226     The \a flag vectorBase is changed to contain 1's and 0's only. A 1
     227     means that the corresponding element in the \a templat vectorBase is
    533228     valid and a zero means that the corresponding element is a NaN.
    534229
    535      \note Space for vector \a flag is reallocated to fit the size of
    536      vector \a templat if sizes mismatch.
    537 
    538      \return True if the \a templat vector contains at least one NaN.
    539   */
    540   bool nan(const vector& templat, vector& flag);
    541 
    542   /**
    543      \brief Transforms a vector to a basis vector.
    544 
    545      All elements are set to zero except the \a i-th element which is
    546      set to one.
    547   */
    548   void set_basis(vector&, size_t i);
    549 
    550   /**
    551      Randomly shuffles the elements in vector \a invec
    552   */
    553   void shuffle(vector& invec);
    554 
    555   /**
    556      Sort the elements in the vector.
    557   */
    558   void sort(vector&);
    559 
    560   /**
    561      Create a vector \a sort_index containing the indeces of
    562      elements in a another vector \a invec.  The elements of \a
    563      sort_index give the index of the vector element which would
    564      have been stored in that position if the vector had been sorted
     230     \note Space for vectorBase \a flag is reallocated to fit the size of
     231     vectorBase \a templat if sizes mismatch.
     232
     233     \return True if the \a templat vectorBase contains at least one NaN.
     234  */
     235  bool nan(const vectorBase& templat, vector& flag);
     236
     237  /**
     238     Create a vectorBase \a sort_index containing the indeces of
     239     elements in a another vectorBase \a invec.  The elements of \a
     240     sort_index give the index of the vectorBase element which would
     241     have been stored in that position if the vectorBase had been sorted
    565242     in place. The first element of \a sort_index gives the index of the least
    566      element in \a invec, and the last element of \a sort_index gives the index of the
    567      greatest element in \a invec . The vector \a invec is not changed.
    568 
    569   */
    570   void sort_index(std::vector<size_t>& sort_index, const vector& invec);
    571 
    572   /** Similar to sort_index but creates a vector with indices to the \a k
     243     element in \a invec, and the last element of \a sort_index gives the
     244     index of the greatest element in \a invec . The vectorBase \a invec
     245     is not changed.
     246  */
     247  void sort_index(std::vector<size_t>& sort_index, const vectorBase& invec);
     248
     249  /**
     250      Similar to sort_index but creates a vectorBase with indices to the \a k
    573251  smallest elements in \a invec. 
    574252  */
    575253  void sort_smallest_index(std::vector<size_t>& sort_index, size_t k, const
    576   vector& invec);
    577 
    578   /** Similar to sort_index but creates a vector with indices to the \a k
     254  vectorBase& invec);
     255
     256  /**
     257      Similar to sort_index but creates a vectorBase with indices to the \a k
    579258  largest elements in \a invec. 
    580259  */
    581260  void sort_largest_index(std::vector<size_t>& sort_index, size_t k, const
    582   vector& invec);
     261  vectorBase& invec);
    583262
    584263 
    585264
    586265  /**
    587      \brief Calculate the sum of all vector elements.
     266     \brief Calculate the sum of all vectorBase elements.
    588267
    589268     \return The sum.
    590269  */
    591   double sum(const vector&);
    592 
    593   /**
    594      \brief Swap vector elements by copying.
    595 
    596      The two vectors must have the same length.
    597 
    598      \throw GSL_error if vector lengths differs.
    599   */
    600   void swap(vector&, vector&);
    601 
    602   /**
    603      \brief The output operator for the vector class.
    604   */
    605   std::ostream& operator<<(std::ostream&, const vector&);
     270  double sum(const vectorBase&);
     271
     272  /**
     273     \brief The output operator for the vectorBase class.
     274  */
     275  std::ostream& operator<<(std::ostream&, const vectorBase&);
    606276
    607277}}} // of namespace utility, yat, and theplu
  • branches/peter-dev/yat/utility/vectorConstView.cc

    r994 r995  
    2525*/
    2626
    27 #include "vector.h"
     27#include "vectorConstView.h"
    2828#include "matrix.h"
    2929#include "utility.h"
     
    4343
    4444
    45   vector::vector(void)
    46     : v_(NULL), view_(NULL), view_const_(NULL), proxy_v_(NULL)
     45  vectorConstView::vectorConstView(void)
     46    : vectorBase(NULL)
    4747  {
    4848  }
    4949
    5050
    51   vector::vector(size_t n, double init_value)
    52     : v_(gsl_vector_alloc(n)), view_(NULL), view_const_(NULL),
    53       proxy_v_(v_)
    54   {
    55     if (!v_)
    56       throw utility::GSL_error("vector::vector failed to allocate memory");
    57 
    58     all(init_value);
    59   }
    60 
    61 
    62   vector::vector(const vector& other)
    63     : v_(other.create_gsl_vector_copy()), view_(NULL),
    64       view_const_(NULL), proxy_v_(v_)
     51  vectorConstView::vectorConstView(const vectorBase& other)
     52    : vectorBase(other.gsl_vector_p())
    6553  {
    6654  }
    6755
    6856
    69   vector::vector(vector& v, size_t offset, size_t n, size_t stride)
    70     : view_const_(NULL)
     57  vectorConstView::vectorConstView(const vectorBase& v, size_t offset,
     58                                   size_t n, size_t stride)
     59    : vectorBase(NULL)
    7160  {
    72     view_ = new gsl_vector_view(gsl_vector_subvector_with_stride(v.v_,offset,
    73                                                                  stride,n));
    74     if (!view_)
    75       throw utility::GSL_error("vector::vector failed to setup view");
    76     proxy_v_ = v_ = &(view_->vector);
     61    const_view_ = new gsl_vector_const_view(
     62                   gsl_vector_const_subvector_with_stride(v.gsl_vector_p(),
     63                                                          offset, stride, n));
     64    const_vec_ = &(const_view_->vector);
    7765  }
    7866
    7967
    80   vector::vector(const vector& v, size_t offset, size_t n, size_t stride)
    81     : v_(NULL), view_(NULL)
     68  vectorConstView::vectorConstView(const matrix& m, size_t i, bool row)
     69    : vectorBase(NULL)
    8270  {
    83     view_const_ = new gsl_vector_const_view(
    84                    gsl_vector_const_subvector_with_stride(v.v_,offset,stride,n));
    85     if (!view_const_)
    86       throw utility::GSL_error("vector::vector failed to setup view");
    87     proxy_v_ = &(view_const_->vector);
     71    const_view_ =
     72      new gsl_vector_const_view(row ?
     73                                gsl_matrix_const_row(m.gsl_matrix_p(),i) :
     74                                gsl_matrix_const_column(m.gsl_matrix_p(),i));
     75    const_vec_ = &(const_view_->vector);
    8876  }
    8977
    9078
    91   vector::vector(matrix& m, size_t i, bool row)
    92     : view_const_(NULL)
     79  vectorConstView::~vectorConstView(void)
    9380  {
    94     view_=new gsl_vector_view(row ?
    95                               gsl_matrix_row   (m.gsl_matrix_p(),i) :
    96                               gsl_matrix_column(m.gsl_matrix_p(),i)  );
    97     if (!view_)
    98       throw utility::GSL_error("vector::vector failed to setup view");
    99     proxy_v_ = v_ = &(view_->vector);
    10081  }
    10182
    10283
    103   vector::vector(const matrix& m, size_t i, bool row)
    104     : v_(NULL), view_(NULL)
    105   {
    106     view_const_ = new gsl_vector_const_view(row ?
    107                                    gsl_matrix_const_row   (m.gsl_matrix_p(),i) :
    108                                    gsl_matrix_const_column(m.gsl_matrix_p(),i) );
    109     if (!view_const_)
    110       throw utility::GSL_error("vector::vector failed to setup view");
    111     proxy_v_ = &(view_const_->vector);
    112   }
    113 
    114 
    115   vector::vector(std::istream& is, char sep)
    116     throw (utility::IO_error, std::exception)
    117     : view_(NULL), view_const_(NULL)
    118   {
    119     // read the data file and store in stl vectors (dynamically
    120     // expandable)
    121     std::vector<std::vector<double> > data_matrix;
    122     u_int nof_columns=0;
    123     u_int nof_rows=0;
    124     std::string line;
    125     while(getline(is, line, '\n')) {
    126       // Empty lines
    127       if (!line.size())
    128           continue;
    129       nof_rows++;
    130 
    131       std::vector<double> v;
    132       std::string element;
    133       std::stringstream ss(line);
    134       bool ok=true;
    135       while(ok) {
    136         if(sep=='\0')
    137           ok=(ss>>element);
    138         else
    139           ok=getline(ss, element, sep);
    140         if(!ok)
    141           break;       
    142 
    143         if(utility::is_double(element)) {
    144           v.push_back(atof(element.c_str()));
    145         }
    146         else if (!element.size() || utility::is_nan(element)) {
    147           v.push_back(std::numeric_limits<double>::quiet_NaN());
    148         }
    149         else {
    150           std::stringstream ss("Warning: '");
    151           ss << element << "' is not a double.";
    152           throw IO_error(ss.str());
    153         }
    154       }
    155       if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
    156           v.push_back(std::numeric_limits<double>::quiet_NaN());
    157       if (!nof_columns)
    158         nof_columns=v.size();
    159       else if (nof_rows && (nof_columns>1)) {
    160         std::ostringstream s;
    161         s << "vector::vector(std::istream&) data file error:\n"
    162           << "    File has inconsistent number of rows (" << nof_rows
    163           << ") and columns (" << nof_columns
    164           << ").\n    Expected a row or a column vector.";
    165         throw utility::IO_error(s.str());
    166       }
    167       else if (v.size()!=nof_columns) {
    168         std::ostringstream s;
    169         s << "vector::vector(std::istream&) data file error:\n"
    170           << "    Line " << nof_rows << " has " << v.size()
    171           << " columns; expected " << nof_columns << " column.";
    172         throw utility::IO_error(s.str());
    173       }
    174       data_matrix.push_back(v);
    175     }
    176 
    177     // manipulate the state of the stream to be good
    178     is.clear(std::ios::goodbit);
    179     // convert the data to a gsl vector
    180     proxy_v_ = v_ = gsl_vector_alloc(nof_rows*nof_columns);
    181     if (!v_)
    182       throw utility::GSL_error("vector::vector failed to allocate memory");
    183     size_t n=0;
    184     // if gsl error handler disabled, out of bounds index will not
    185     // abort the program.
    186     for (size_t i=0; i<nof_rows; i++)
    187       for (size_t j=0; j<nof_columns; j++)
    188         gsl_vector_set( v_, n++, data_matrix[i][j] );
    189   }
    190 
    191 
    192   vector::~vector(void)
    193   {
    194     delete_allocated_memory();
    195   }
    196 
    197 
    198   vector::iterator vector::begin(void)
    199   {
    200     return vector::iterator(*this, 0);
    201   }
    202 
    203 
    204   vector::const_iterator vector::begin(void) const
    205   {
    206     return vector::const_iterator(*this, 0);
    207   }
    208 
    209 
    210   const vector& vector::clone(const vector& other)
    211   {
    212     if (this!=&other) {
    213 
    214       delete_allocated_memory();
    215 
    216       if (other.view_) {
    217         view_ = new gsl_vector_view(*other.view_);
    218         proxy_v_ = v_ = &(view_->vector);
    219       }
    220       else if (other.view_const_) {
    221         view_const_ = new gsl_vector_const_view(*other.view_const_);
    222         proxy_v_ = &(view_const_->vector);
    223         v_=NULL;
    224       }
    225       else if (other.v_)
    226         proxy_v_ = v_ = other.create_gsl_vector_copy();
    227 
    228     }
    229     return *this;
    230   }
    231 
    232 
    233   gsl_vector* vector::create_gsl_vector_copy(void) const
    234   {
    235     gsl_vector* vec = gsl_vector_alloc(size());
    236     if (!vec)
    237       throw utility::GSL_error("vector::create_gsl_vector_copy failed to allocate memory");
    238     if (gsl_vector_memcpy(vec, proxy_v_))
    239       throw utility::GSL_error("vector::create_gsl_matrix_copy dimension mis-match");
    240     return vec;
    241   }
    242 
    243 
    244   void vector::delete_allocated_memory(void)
    245   {
    246     if (view_)
    247       delete view_;
    248     else if (view_const_)
    249       delete view_const_;
    250     else if (v_)
    251       gsl_vector_free(v_);
    252     proxy_v_=v_=NULL;
    253   }
    254 
    255 
    256   void vector::div(const vector& other)
    257   {
    258     assert(v_);
    259     int status=gsl_vector_div(v_,other.gsl_vector_p());
    260     if (status)
    261       throw utility::GSL_error(std::string("vector::div",status));
    262   }
    263 
    264 
    265   vector::iterator vector::end(void)
    266   {
    267     return vector::iterator(*this, size());
    268   }
    269 
    270 
    271   vector::const_iterator vector::end(void) const
    272   {
    273     return vector::const_iterator(*this, size());
    274   }
    275 
    276 
    277   bool vector::equal(const vector& other, const double d) const
    278   {
    279     if (this==&other)
    280       return true;
    281     if (size()!=other.size())
    282       return false;
    283     // if gsl error handler disabled, out of bounds index will not
    284     // abort the program.
    285     for (size_t i=0; i<size(); ++i)
    286       // The two last condition checks are needed for NaN detection
    287       if (fabs( (*this)(i)-other(i) ) > d ||
    288           (*this)(i)!=(*this)(i) || other(i)!=other(i))
    289         return false;
    290     return true;
    291   }
    292 
    293 
    294   const gsl_vector* vector::gsl_vector_p(void) const
    295   {
    296     return proxy_v_;
    297   }
    298 
    299 
    300   gsl_vector* vector::gsl_vector_p(void)
    301   {
    302     return v_;
    303   }
    304 
    305 
    306   bool vector::isview(void) const
    307   {
    308     return view_ || view_const_;
    309   }
    310 
    311 
    312   void vector::mul(const vector& other)
    313   {
    314     assert(v_);
    315     int status=gsl_vector_mul(v_,other.gsl_vector_p());
    316     if (status)
    317       throw utility::GSL_error(std::string("vector::mul",status));
    318   }
    319 
    320 
    321   void vector::resize(size_t n, double init_value)
    322   {
    323     delete_allocated_memory();
    324     proxy_v_ = v_ = gsl_vector_alloc(n);
    325     if (!v_)
    326       throw utility::GSL_error("vector::vector failed to allocate memory");
    327     all(init_value);
    328   }
    329 
    330 
    331   void vector::reverse(void)
    332   {
    333     assert(v_);
    334     gsl_vector_reverse(v_);
    335   }
    336 
    337 
    338   void vector::all(const double& value)
    339   {
    340     assert(v_);
    341     gsl_vector_set_all(v_,value);
    342   }
    343 
    344 
    345   size_t vector::size(void) const
    346   {
    347     if (!proxy_v_)
    348       return 0;
    349     return proxy_v_->size;
    350   }
    351 
    352 
    353   void vector::swap(size_t i, size_t j)
    354   {
    355     assert(v_);
    356     int status=gsl_vector_swap_elements(v_, i, j);
    357     if (status)
    358       throw utility::GSL_error(std::string("vector::swap_elements",status));
    359   }
    360 
    361 
    362   double& vector::operator()(size_t i)
    363   {
    364     assert(v_);
    365     double* d=gsl_vector_ptr(v_, i);
    366     if (!d)
    367       throw utility::GSL_error("vector::operator()",GSL_EINVAL);
    368     return *d;
    369   }
    370 
    371 
    372   const double& vector::operator()(size_t i) const
    373   {
    374     const double* d=gsl_vector_const_ptr(proxy_v_, i);
    375     if (!d)
    376       throw utility::GSL_error("vector::operator()",GSL_EINVAL);
    377     return *d;
    378   }
    379 
    380 
    381   double& vector::operator[](size_t i)
    382   {
    383     return this->operator()(i);
    384   }
    385 
    386 
    387   const double& vector::operator[](size_t i) const
    388   {
    389     return this->operator()(i);
    390   }
    391 
    392 
    393   bool vector::operator==(const vector& other) const
    394   {
    395     return equal(other);
    396   }
    397 
    398 
    399   bool vector::operator!=(const vector& other) const
    400   {
    401     return !equal(other);
    402   }
    403 
    404 
    405   double vector::operator*( const vector &other ) const
    406   {
    407     double res = 0.0;;
    408     for ( size_t i = 0; i < size(); ++i )
    409       res += other(i) * (*this)(i);
    410     return res;
    411   }
    412 
    413 
    414   const vector& vector::operator=( const vector& other )
    415   {
    416     assert(v_);
    417     if (gsl_vector_memcpy(v_,other.gsl_vector_p()))
    418       throw utility::GSL_error("vector::set dimension mis-match");
    419     return *this;
    420   }
    421 
    422 
    423   const vector& vector::operator+=(const vector& other)
    424   {
    425     assert(v_);
    426     int status=gsl_vector_add(v_, other.gsl_vector_p());
    427     if (status)
    428       throw utility::GSL_error(std::string("vector::add", status));
    429     return *this;
    430   }
    431 
    432 
    433   const vector& vector::operator+=(double d)
    434   {
    435     assert(v_);
    436     gsl_vector_add_constant(v_, d);
    437     return *this;
    438   }
    439 
    440 
    441   const vector& vector::operator-=(const vector& other)
    442   {
    443     assert(v_);
    444     int status=gsl_vector_sub(v_, other.gsl_vector_p());
    445     if (status)
    446       throw utility::GSL_error(std::string("vector::sub", status));
    447     return *this;
    448   }
    449 
    450 
    451   const vector& vector::operator-=(const double d)
    452   {
    453     assert(v_);
    454     gsl_vector_add_constant(v_, -d);
    455     return *this;
    456   }
    457 
    458 
    459   const vector& vector::operator*=(const double d)
    460   {
    461     assert(v_);
    462     gsl_vector_scale(v_, d);
    463     return *this;
    464   }
    465 
    466 
    467   bool isnull(const vector& v)
    468   {
    469     return gsl_vector_isnull(v.gsl_vector_p());
    470   }
    471 
    472 
    473   double max(const vector& v)
    474   {
    475     return gsl_vector_max(v.gsl_vector_p());
    476   }
    477 
    478 
    479   size_t max_index(const vector& v)
    480   {
    481     return gsl_vector_max_index(v.gsl_vector_p());
    482   }
    483 
    484 
    485   double min(const vector& v)
    486   {
    487     return gsl_vector_min(v.gsl_vector_p());
    488   }
    489 
    490 
    491   size_t min_index(const vector& v)
    492   {
    493     return gsl_vector_min_index(v.gsl_vector_p());
    494   }
    495 
    496 
    497   bool nan(const vector& templat, vector& flag)
    498   {
    499     size_t vsize=templat.size();
    500     if (vsize!=flag.size())
    501       flag.clone(vector(vsize,1.0));
    502     else
    503       flag.all(1.0);
    504     bool nan=false;
    505     for (size_t i=0; i<vsize; i++)
    506       if (std::isnan(templat(i))) {
    507         flag(i)=0;
    508         nan=true;
    509       }
    510     return nan;
    511   }
    512 
    513 
    514   void set_basis(vector& v, size_t i)
    515   {
    516     assert(v.gsl_vector_p());
    517     gsl_vector_set_basis(v.gsl_vector_p(),i);
    518   }
    519 
    520 
    521   void shuffle(vector& invec)
    522   {
    523     random::DiscreteUniform rnd;
    524     std::random_shuffle(invec.begin(), invec.end(), rnd);
    525   }
    526 
    527 
    528   void sort(vector& v)
    529   {
    530     assert(v.gsl_vector_p());
    531     std::sort(v.begin(), v.end());
    532   }
    533 
    534 
    535   void sort_index(std::vector<size_t>& sort_index, const vector& invec)
    536   {
    537     assert(invec.gsl_vector_p());
    538     gsl_permutation* p = gsl_permutation_alloc(invec.size());
    539     int status=gsl_sort_vector_index(p,invec.gsl_vector_p());
    540     if (status) {
    541       gsl_permutation_free(p);
    542       throw utility::GSL_error(std::string("sort_index(vector&,const vector&)",status));     
    543     }
    544     sort_index=std::vector<size_t>(p->data,p->data+p->size);
    545     gsl_permutation_free(p);
    546   }
    547 
    548 
    549   void sort_smallest_index(std::vector<size_t>& sort_index, size_t k,
    550                             const vector& invec)
    551   {
    552     assert(invec.gsl_vector_p());
    553     assert(k<=invec.size());
    554     sort_index.resize(k);
    555     gsl_sort_vector_smallest_index(&sort_index[0],k,invec.gsl_vector_p());
    556   }
    557  
    558   void sort_largest_index(std::vector<size_t>& sort_index, size_t k,
    559                             const vector& invec)
    560   {
    561     assert(invec.gsl_vector_p());
    562     assert(k<=invec.size());
    563     sort_index.resize(k);
    564     gsl_sort_vector_largest_index(&sort_index[0],k,invec.gsl_vector_p());
    565   }
    566 
    567 
    568   double sum(const vector& v)
    569   {
    570     double sum = 0;
    571     size_t vsize=v.size();
    572     for (size_t i=0; i<vsize; ++i)
    573       sum += v[i];
    574     return sum;
    575   }
    576 
    577 
    578   void swap(vector& v, vector& w)
    579   {
    580     assert(v.gsl_vector_p()); assert(w.gsl_vector_p());
    581     int status=gsl_vector_swap(v.gsl_vector_p(),w.gsl_vector_p());
    582     if (status)
    583       throw utility::GSL_error(std::string("swap(vector&,vector&)",status));
    584   }
    585 
    586 
    587   std::ostream& operator<<(std::ostream& s, const vector& a)
    588   {
    589     s.setf(std::ios::dec);
    590     s.precision(12);
    591     for (size_t j = 0; j < a.size(); ++j) {
    592       s << a[j];
    593       if ( (j+1)<a.size() )
    594         s << s.fill();
    595     }
    596     return s;
    597   }
    598 
    59984}}} // of namespace utility, yat, and thep
  • branches/peter-dev/yat/utility/vectorConstView.h

    r994 r995  
    1 #ifndef _theplu_yat_utility_vector_
    2 #define _theplu_yat_utility_vector_
     1#ifndef _theplu_yat_utility_vector_const_view_
     2#define _theplu_yat_utility_vector_const_view_
    33
    44// $Id$
     
    2828  02111-1307, USA.
    2929*/
     30
     31#include "vectorBase.h"
    3032
    3133#include "Exception.h"
     
    8486  */
    8587
    86   class vector
     88  class vectorConstView : public vectorBase
    8789  {
    8890  public:
    8991
    90     /// \brief vector::iterator
    91     typedef Iterator<double&, vector> iterator;
    92     /// \brief vector::const_iterator
    93     typedef Iterator<const double, const vector> const_iterator;
    94 
    9592    /**
    96        \brief The default constructor.
     93       \brief Default Constructor
    9794    */
    98     vector(void);
    99 
    100     /**
    101        \brief Allocates memory space for \a n elements, and sets all
    102        elements to \a init_value.
    103        
    104        \throw GSL_error if memory allocation fails.
    105     */
    106     vector(size_t n, double init_value=0);
     95    vectorConstView(void);
    10796
    10897    /**
    10998       \brief The copy constructor.
    11099
    111        \note If the object to be copied is a vector view, the values
     100       \note If the object to be copied is a vectorConstView view, the values
    112101       of the view will be copied, i.e. the view is not copied.
    113102
     
    115104       fails.
    116105    */
    117     vector(const vector& other);
     106    vectorConstView(const vectorBase& other);
    118107
    119108    /**
    120        \brief Vector view constructor.
     109       \brief VectorConstView constructor.
    121110
    122        Create a view of vector \a v, with starting index \a offset,
     111       Create a ConstView of vectorBase \a v, with starting index \a offset,
    123112       size \a n, and an optional \a stride.
    124 
    125        A vector view can be used as any vector with the difference
    126        that changes made to the view will also change the object that
    127        is viewed. Also, using the copy constructor will create a new
    128        vector object that is a copy of whatever is viewed. If a copy
    129        of the view is needed then you should use this constructor to
    130        obtain a copy.
    131113
    132114       \note If the object viewed by the view goes out of scope or is
     
    136118       \throw GSL_error if a view cannot be set up.
    137119    */
    138     vector(vector& v, size_t offset, size_t n, size_t stride=1);
    139 
    140     /**
    141        \brief Vector const view constructor.
    142 
    143        Create a view of vector \a v, with starting index \a offset,
    144        size \a n, and an optional \a stride.
    145 
    146        A vector view can be used as any const vector. Using the copy
    147        constructor will create a new vector object that is a copy of
    148        whatever is viewed. If a copy of the view is needed then you
    149        should use this constructor to obtain a copy.
    150 
    151        \note If the object viewed by the view goes out of scope or is
    152        deleted, the view becomes invalid and the result of further use
    153        is undefined.
    154 
    155        \throw GSL_error if a view cannot be set up.
    156     */
    157     vector(const vector& v, size_t offset, size_t n, size_t stride=1);
     120    vectorConstView(const vectorBase& v, size_t offset, size_t n,
     121                    size_t stride=1);
    158122
    159123    ///
    160     /// Matrix row/column view constructor.
     124    /// Matrix row/column const view constructor.
    161125    ///
    162     /// Create a row/column vector view of matrix \a m, pointing at
     126    /// Create a row/column vectorConstView view of matrix \a m, pointing at
    163127    /// row/column \a i. The parameter \a row is used to set whether
    164128    /// the view should be a row or column view. If \a row is set to
     
    166130    /// naturally, a column view otherwise.
    167131    ///
    168     /// A vector view can be used as any vector with the difference
    169     /// that changes made to the view will also change the object that
    170     /// is viewed. Also, using the copy constructor will create a new
    171     /// vector object that is a copy of whatever is viewed. If a copy
    172     /// of the view is needed then you should use the vector view
    173     /// constructor to obtain a copy.
    174     ///
    175132    /// @note If the object viewed by the view goes out of scope or is
    176133    /// deleted, the view becomes invalid and the result of further
    177134    /// use is undefined.
    178135    ///
    179     vector(matrix& m, size_t i, bool row=true);
    180 
    181     ///
    182     /// Matrix row/column const view constructor.
    183     ///
    184     /// Create a row/column vector view of matrix \a m, pointing at
    185     /// row/column \a i. The parameter \a row is used to set whether
    186     /// the view should be a row or column view. If \a row is set to
    187     /// true, the view will be a row view (default behaviour), and,
    188     /// naturally, a column view otherwise.
    189     ///
    190     /// A const vector view can be used as any const vector. Using the
    191     /// copy constructor will create a new vector object that is a
    192     /// copy of whatever is viewed. If a copy of the view is needed
    193     /// then you should use the vector view constructor to obtain a
    194     /// copy.
    195     ///
    196     /// @note If the object viewed by the view goes out of scope or is
    197     /// deleted, the view becomes invalid and the result of further
    198     /// use is undefined.
    199     ///
    200     vector(const matrix& m, size_t i, bool row=true);
    201 
    202     /**
    203        \brief The istream constructor.
    204 
    205        Either elements should be separated with white space characters
    206        (default), or elements should be separated by the delimiter \a
    207        sep. When delimiter \a sep is used empty elements are stored as
    208        NaN's (except that empty lines are ignored). The end of input
    209        to the vector is at end of file marker.
    210 
    211        \throw GSL_error if memory allocation fails, IO_error if
    212        unexpected input is found in the input stream.
    213     */
    214     explicit vector(std::istream &, char sep='\0')
    215       throw(utility::IO_error, std::exception);
     136    vectorConstView(const matrix& m, size_t i, bool row=true);
    216137
    217138    ///
    218139    /// The destructor.
    219140    ///
    220     ~vector(void);
    221 
    222     /**
    223        \return mutable iterator to start of vector
    224      */
    225     iterator begin(void);
    226 
    227     /**
    228        \return read-only iterator to start of vector
    229      */
    230     const_iterator begin(void) const;
    231 
    232     /**
    233        \brief Make a copy of \a other.
    234 
    235        This function will make a deep copy of \a other. Memory is
    236        resized and view state is changed if needed.
    237     */
    238     const vector& clone(const vector& other);
    239 
    240     /**
    241        \brief This function performs element-wise division, \f$ this_i =
    242        this_i/other_i \; \forall i \f$.
    243 
    244        \throw GSL_error if dimensions mis-match.
    245     */
    246     void div(const vector& other);
    247 
    248     /**
    249        \return mutable iterator to end of vector
    250      */
    251     iterator end(void);
    252 
    253     /**
    254        \return read-only iterator to end of vector
    255      */
    256     const_iterator end(void) const;
    257 
    258     /**
    259        \brief Check whether vectors are equal within a user defined
    260        precision, set by \a precision.
    261 
    262        \return True if each element deviates less or equal than \a
    263        d. If any vector contain a NaN, false is always returned.
    264 
    265        \see operator== and operator!=
    266     */
    267     bool equal(const vector&, const double precision=0) const;
    268 
    269     ///
    270     /// @return A const pointer to the internal GSL vector,
    271     ///
    272     const gsl_vector* gsl_vector_p(void) const;
    273 
    274     ///
    275     /// @return A pointer to the internal GSL vector,
    276     ///
    277     gsl_vector* gsl_vector_p(void);
    278 
    279     ///
    280     /// Check if the vector object is a view (sub-vector) to another
    281     /// vector.
    282     ///
    283     /// @return True if the object is a view, false othwerwise.
    284     ///
    285     bool isview(void) const;
    286 
    287     /**
    288        \brief This function performs element-wise multiplication, \f$
    289        this_i = this_i * other_i \; \forall i \f$.
    290 
    291        \throw GSL_error if dimensions mis-match.
    292     */
    293     void mul(const vector& other);
    294 
    295     /**
    296        @brief Resize vector
    297        
    298        All elements are set to @a init_value.
    299 
    300        \note Underlying GSL vector is destroyed and a view into this
    301        vector becomes invalid.
    302     */
    303     void resize(size_t, double init_value=0);
    304 
    305     /**
    306        \brief Reverse the order of elements in the vector.
    307     */
    308     void reverse(void);
    309 
    310     ///
    311     /// Set all elements to \a value.
    312     ///
    313     void all(const double& value);
    314 
    315     ///
    316     /// @return the number of elements in the vector.
    317     ///
    318     size_t size(void) const;
    319 
    320     /**
    321        \brief Exchange elements \a i and \a j.
    322 
    323        \throw GSL_error if vector lengths differs.
    324     */
    325     void swap(size_t i, size_t j);
    326 
    327     /**
    328        \brief Element access operator.
    329 
    330        \return Reference to element \a i.
    331 
    332        \throw If GSL range checks are enabled in the underlying GSL
    333        library a GSL_error exception is thrown if either index is out
    334        of range.
    335     */
    336     double& operator()(size_t i);
    337 
    338     /**
    339        \brief Element access operator.
    340 
    341        \return Const reference to element \a i.
    342 
    343        \throw If GSL range checks are enabled in the underlying GSL
    344        library a GSL_error exception is thrown if either index is out
    345        of range.
    346     */
    347     const double& operator()(size_t i) const;
    348 
    349     ///
    350     /// Element access operator.
    351     ///
    352     /// @return Reference to element \a i.
    353     ///
    354     double& operator[](size_t i);
    355 
    356     ///
    357     /// Const element access operator.
    358     ///
    359     /// @return The value of element \a i.
    360     ///
    361     const double& operator[](size_t i) const;
    362 
    363     /**
    364        \brief Comparison operator. Takes linear time.
    365 
    366        Checks are performed with exact matching, i.e., rounding off
    367        effects may destroy comparison. Use the equal function for
    368        comparing elements within a user defined precision.
    369 
    370        \return True if all elements are equal otherwise false.
    371 
    372        \see equal(const vector&, const double precision=0)
    373     */
    374     bool operator==(const vector&) const;
    375 
    376     /**
    377        \brief Comparison operator. Takes linear time.
    378 
    379        Checks are performed with exact matching, i.e., rounding off
    380        effects may destroy comparison. Use the equal function for
    381        comparing elements within a user defined precision.
    382 
    383        \return False if all elements are equal otherwise true.
    384 
    385        \see equal(const vector&, const double precision=0)
    386     */
    387     bool operator!=(const vector&) const;
    388 
    389     ///
    390     /// @return The dot product.
    391     ///
    392     double operator*(const vector&) const;
    393 
    394     /**
    395        \brief The assignment operator.
    396 
    397        Dimensions of the vectors must match. If the LHS vector is a
    398        view, the underlying data will be changed.
    399 
    400        \return A const reference to the resulting vector.
    401 
    402        \see void set(const vector&).
    403 
    404        \throw GSL_error if dimensions mis-match.
    405     */
    406     const vector& operator=(const vector&);
    407 
    408     /**
    409        \brief Addition and assign operator. Vector addition, \f$
    410        this_i = this_i + other_i \; \forall i \f$.
    411 
    412        \return A const reference to the resulting vector.
    413 
    414        \throw GSL_error if dimensions mis-match.
    415     */
    416     const vector& operator+=(const vector&);
    417 
    418     /**
    419        \brief Add a constant to a vector, \f$ this_i = this_i + d \;
    420        \forall i \f$.
    421 
    422        \return A const reference to the resulting vector.
    423     */
    424     const vector& operator+=(double d);
    425 
    426     /**
    427        \brief Subtract and assign operator. Vector subtraction, \f$
    428        this_i = this_i - other_i \; \forall i \f$.
    429 
    430        \return A const reference to the resulting vector.
    431 
    432        \throw GSL_error if dimensions mis-match.
    433     */
    434     const vector& operator-=(const vector&);
    435 
    436     /**
    437        \brief Subtract a constant to a vector, \f$ this_i = this_i - d
    438        \; \forall i \f$.
    439 
    440        \return A const reference to the resulting vector.
    441     */
    442     const vector& operator-=(double d);
    443 
    444     /**
    445        \brief Multiply with scalar and assign operator, \f$ this_i =
    446        this_i * d \; \forall i \f$.
    447 
    448        \return A const reference to the resulting vector.
    449     */
    450     const vector& operator*=(double d);
    451 
     141    ~vectorConstView(void);
    452142
    453143  private:
    454 
    455     /**
    456        \brief Create a new copy of the internal GSL vector.
    457 
    458        Necessary memory for the new GSL vector is allocated and the
    459        caller is responsible for freeing the allocated memory.
    460 
    461        \return A pointer to a copy of the internal GSL vector.
    462 
    463        \throw GSL_error if memory cannot be allocated for the new
    464        copy, or if dimensions mis-match.
    465     */
    466     gsl_vector* create_gsl_vector_copy(void) const;
    467 
    468     /**
    469        \brief Clear all dynamically allocated memory.
    470 
    471        Internal utility function.
    472     */
    473     void delete_allocated_memory(void);
    474 
    475     gsl_vector* v_;
    476     gsl_vector_view* view_;
    477     const gsl_vector_const_view* view_const_;
    478     // proxy_v_ is used to access the proper underlying gsl_vector
    479     // in all const member functions. It is not used by design for
    480     // non-const vector functions and operators. This is to make sure
    481     // that runtime errors occur if a const vector is used in an
    482     // inappropriate manner such as on left hand side in assignment
    483     // (remember, v_ is null for const vector views).
    484     const gsl_vector* proxy_v_;
     144    const gsl_vector_const_view* const_view_;
    485145  };
    486146
    487   /**
    488      \brief Check if all elements of the vector are zero.
    489 
    490      \return True if all elements in the vector is zero, false
    491      othwerwise.
    492   */
    493   bool isnull(const vector&);
    494 
    495   /**
    496      \brief Get the maximum value of the vector.
    497 
    498      \return The maximum value of the vector.
    499   */
    500   double max(const vector&);
    501 
    502   /**
    503      \brief Locate the maximum value in the vector.
    504 
    505      \return The index to the maximum value of the vector.
    506 
    507      \note Lower index has precedence.
    508   */
    509   size_t max_index(const vector&);
    510 
    511   /**
    512      \brief Get the minimum value of the vector.
    513 
    514      \return The minimum value of the vector.
    515   */
    516   double min(const vector&);
    517 
    518   /**
    519      \brief Locate the minimum value in the vector.
    520 
    521      \return The index to the minimum value of the vector.
    522 
    523      \note Lower index has precedence.
    524   */
    525   size_t min_index(const vector&);
    526 
    527   /**
    528      \brief Create a vector \a flag indicating NaN's in another vector
    529      \a templat.
    530 
    531      The \a flag vector is changed to contain 1's and 0's only. A 1
    532      means that the corresponding element in the \a templat vector is
    533      valid and a zero means that the corresponding element is a NaN.
    534 
    535      \note Space for vector \a flag is reallocated to fit the size of
    536      vector \a templat if sizes mismatch.
    537 
    538      \return True if the \a templat vector contains at least one NaN.
    539   */
    540   bool nan(const vector& templat, vector& flag);
    541 
    542   /**
    543      \brief Transforms a vector to a basis vector.
    544 
    545      All elements are set to zero except the \a i-th element which is
    546      set to one.
    547   */
    548   void set_basis(vector&, size_t i);
    549 
    550   /**
    551      Randomly shuffles the elements in vector \a invec
    552   */
    553   void shuffle(vector& invec);
    554 
    555   /**
    556      Sort the elements in the vector.
    557   */
    558   void sort(vector&);
    559 
    560   /**
    561      Create a vector \a sort_index containing the indeces of
    562      elements in a another vector \a invec.  The elements of \a
    563      sort_index give the index of the vector element which would
    564      have been stored in that position if the vector had been sorted
    565      in place. The first element of \a sort_index gives the index of the least
    566      element in \a invec, and the last element of \a sort_index gives the index of the
    567      greatest element in \a invec . The vector \a invec is not changed.
    568 
    569   */
    570   void sort_index(std::vector<size_t>& sort_index, const vector& invec);
    571 
    572   /** Similar to sort_index but creates a vector with indices to the \a k
    573   smallest elements in \a invec. 
    574   */
    575   void sort_smallest_index(std::vector<size_t>& sort_index, size_t k, const
    576   vector& invec);
    577 
    578   /** Similar to sort_index but creates a vector with indices to the \a k
    579   largest elements in \a invec. 
    580   */
    581   void sort_largest_index(std::vector<size_t>& sort_index, size_t k, const
    582   vector& invec);
    583 
    584  
    585 
    586   /**
    587      \brief Calculate the sum of all vector elements.
    588 
    589      \return The sum.
    590   */
    591   double sum(const vector&);
    592 
    593   /**
    594      \brief Swap vector elements by copying.
    595 
    596      The two vectors must have the same length.
    597 
    598      \throw GSL_error if vector lengths differs.
    599   */
    600   void swap(vector&, vector&);
    601 
    602   /**
    603      \brief The output operator for the vector class.
    604   */
    605   std::ostream& operator<<(std::ostream&, const vector&);
    606 
    607 }}} // of namespace utility, yat, and theplu
     147}}}  // of namespace utility, yat, and theplu
    608148
    609149#endif
  • branches/peter-dev/yat/utility/vectorMutable.cc

    r994 r995  
    2525*/
    2626
    27 #include "vector.h"
     27#include "vectorMutable.h"
    2828#include "matrix.h"
    2929#include "utility.h"
     
    4343
    4444
    45   vector::vector(void)
    46     : v_(NULL), view_(NULL), view_const_(NULL), proxy_v_(NULL)
     45  vectorMutable::vectorMutable(gsl_vector* v)
     46    : vectorBase(v)
    4747  {
    4848  }
    4949
    5050
    51   vector::vector(size_t n, double init_value)
    52     : v_(gsl_vector_alloc(n)), view_(NULL), view_const_(NULL),
    53       proxy_v_(v_)
    54   {
    55     if (!v_)
    56       throw utility::GSL_error("vector::vector failed to allocate memory");
    57 
    58     all(init_value);
    59   }
    60 
    61 
    62   vector::vector(const vector& other)
    63     : v_(other.create_gsl_vector_copy()), view_(NULL),
    64       view_const_(NULL), proxy_v_(v_)
     51  vectorMutable::~vectorMutable(void)
    6552  {
    6653  }
    6754
    6855
    69   vector::vector(vector& v, size_t offset, size_t n, size_t stride)
    70     : view_const_(NULL)
     56  vectorMutable::iterator vectorMutable::begin(void)
    7157  {
    72     view_ = new gsl_vector_view(gsl_vector_subvector_with_stride(v.v_,offset,
    73                                                                  stride,n));
    74     if (!view_)
    75       throw utility::GSL_error("vector::vector failed to setup view");
    76     proxy_v_ = v_ = &(view_->vector);
     58    return vectorMutable::iterator(*this, 0);
    7759  }
    7860
    7961
    80   vector::vector(const vector& v, size_t offset, size_t n, size_t stride)
    81     : v_(NULL), view_(NULL)
     62  void vectorMutable::div(const vectorBase& other)
    8263  {
    83     view_const_ = new gsl_vector_const_view(
    84                    gsl_vector_const_subvector_with_stride(v.v_,offset,stride,n));
    85     if (!view_const_)
    86       throw utility::GSL_error("vector::vector failed to setup view");
    87     proxy_v_ = &(view_const_->vector);
     64    assert(vec_);
     65    int status=gsl_vector_div(vec_,other.gsl_vector_p());
     66    if (status)
     67      throw utility::GSL_error(std::string("vectorMutable::div",status));
    8868  }
    8969
    9070
    91   vector::vector(matrix& m, size_t i, bool row)
    92     : view_const_(NULL)
     71  vectorMutable::iterator vectorMutable::end(void)
    9372  {
    94     view_=new gsl_vector_view(row ?
    95                               gsl_matrix_row   (m.gsl_matrix_p(),i) :
    96                               gsl_matrix_column(m.gsl_matrix_p(),i)  );
    97     if (!view_)
    98       throw utility::GSL_error("vector::vector failed to setup view");
    99     proxy_v_ = v_ = &(view_->vector);
     73    return vectorMutable::iterator(*this, size());
    10074  }
    10175
    10276
    103   vector::vector(const matrix& m, size_t i, bool row)
    104     : v_(NULL), view_(NULL)
     77  gsl_vector* vectorMutable::gsl_vector_p(void)
    10578  {
    106     view_const_ = new gsl_vector_const_view(row ?
    107                                    gsl_matrix_const_row   (m.gsl_matrix_p(),i) :
    108                                    gsl_matrix_const_column(m.gsl_matrix_p(),i) );
    109     if (!view_const_)
    110       throw utility::GSL_error("vector::vector failed to setup view");
    111     proxy_v_ = &(view_const_->vector);
     79    return vec_;
    11280  }
    11381
    11482
    115   vector::vector(std::istream& is, char sep)
    116     throw (utility::IO_error, std::exception)
    117     : view_(NULL), view_const_(NULL)
     83  void vectorMutable::mul(const vectorBase& other)
    11884  {
    119     // read the data file and store in stl vectors (dynamically
    120     // expandable)
    121     std::vector<std::vector<double> > data_matrix;
    122     u_int nof_columns=0;
    123     u_int nof_rows=0;
    124     std::string line;
    125     while(getline(is, line, '\n')) {
    126       // Empty lines
    127       if (!line.size())
    128           continue;
    129       nof_rows++;
    130 
    131       std::vector<double> v;
    132       std::string element;
    133       std::stringstream ss(line);
    134       bool ok=true;
    135       while(ok) {
    136         if(sep=='\0')
    137           ok=(ss>>element);
    138         else
    139           ok=getline(ss, element, sep);
    140         if(!ok)
    141           break;       
    142 
    143         if(utility::is_double(element)) {
    144           v.push_back(atof(element.c_str()));
    145         }
    146         else if (!element.size() || utility::is_nan(element)) {
    147           v.push_back(std::numeric_limits<double>::quiet_NaN());
    148         }
    149         else {
    150           std::stringstream ss("Warning: '");
    151           ss << element << "' is not a double.";
    152           throw IO_error(ss.str());
    153         }
    154       }
    155       if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
    156           v.push_back(std::numeric_limits<double>::quiet_NaN());
    157       if (!nof_columns)
    158         nof_columns=v.size();
    159       else if (nof_rows && (nof_columns>1)) {
    160         std::ostringstream s;
    161         s << "vector::vector(std::istream&) data file error:\n"
    162           << "    File has inconsistent number of rows (" << nof_rows
    163           << ") and columns (" << nof_columns
    164           << ").\n    Expected a row or a column vector.";
    165         throw utility::IO_error(s.str());
    166       }
    167       else if (v.size()!=nof_columns) {
    168         std::ostringstream s;
    169         s << "vector::vector(std::istream&) data file error:\n"
    170           << "    Line " << nof_rows << " has " << v.size()
    171           << " columns; expected " << nof_columns << " column.";
    172         throw utility::IO_error(s.str());
    173       }
    174       data_matrix.push_back(v);
    175     }
    176 
    177     // manipulate the state of the stream to be good
    178     is.clear(std::ios::goodbit);
    179     // convert the data to a gsl vector
    180     proxy_v_ = v_ = gsl_vector_alloc(nof_rows*nof_columns);
    181     if (!v_)
    182       throw utility::GSL_error("vector::vector failed to allocate memory");
    183     size_t n=0;
    184     // if gsl error handler disabled, out of bounds index will not
    185     // abort the program.
    186     for (size_t i=0; i<nof_rows; i++)
    187       for (size_t j=0; j<nof_columns; j++)
    188         gsl_vector_set( v_, n++, data_matrix[i][j] );
     85    assert(vec_);
     86    int status=gsl_vector_mul(vec_,other.gsl_vector_p());
     87    if (status)
     88      throw utility::GSL_error(std::string("vectorMutable::mul",status));
    18989  }
    19090
    19191
    192   vector::~vector(void)
     92  void vectorMutable::reverse(void)
    19393  {
    194     delete_allocated_memory();
     94    assert(vec_);
     95    gsl_vector_reverse(vec_);
    19596  }
    19697
    19798
    198   vector::iterator vector::begin(void)
     99  void vectorMutable::all(const double& value)
    199100  {
    200     return vector::iterator(*this, 0);
     101    assert(vec_);
     102    gsl_vector_set_all(vec_,value);
    201103  }
    202104
    203105
    204   vector::const_iterator vector::begin(void) const
     106  void vectorMutable::swap(size_t i, size_t j)
    205107  {
    206     return vector::const_iterator(*this, 0);
     108    assert(vec_);
     109    int status=gsl_vector_swap_elements(vec_, i, j);
     110    if (status)
     111      throw utility::GSL_error(std::string("vectorMutable::swap_elements",status));
    207112  }
    208113
    209114
    210   const vector& vector::clone(const vector& other)
     115  double& vectorMutable::operator()(size_t i)
    211116  {
    212     if (this!=&other) {
    213 
    214       delete_allocated_memory();
    215 
    216       if (other.view_) {
    217         view_ = new gsl_vector_view(*other.view_);
    218         proxy_v_ = v_ = &(view_->vector);
    219       }
    220       else if (other.view_const_) {
    221         view_const_ = new gsl_vector_const_view(*other.view_const_);
    222         proxy_v_ = &(view_const_->vector);
    223         v_=NULL;
    224       }
    225       else if (other.v_)
    226         proxy_v_ = v_ = other.create_gsl_vector_copy();
    227 
    228     }
    229     return *this;
    230   }
    231 
    232 
    233   gsl_vector* vector::create_gsl_vector_copy(void) const
    234   {
    235     gsl_vector* vec = gsl_vector_alloc(size());
    236     if (!vec)
    237       throw utility::GSL_error("vector::create_gsl_vector_copy failed to allocate memory");
    238     if (gsl_vector_memcpy(vec, proxy_v_))
    239       throw utility::GSL_error("vector::create_gsl_matrix_copy dimension mis-match");
    240     return vec;
    241   }
    242 
    243 
    244   void vector::delete_allocated_memory(void)
    245   {
    246     if (view_)
    247       delete view_;
    248     else if (view_const_)
    249       delete view_const_;
    250     else if (v_)
    251       gsl_vector_free(v_);
    252     proxy_v_=v_=NULL;
    253   }
    254 
    255 
    256   void vector::div(const vector& other)
    257   {
    258     assert(v_);
    259     int status=gsl_vector_div(v_,other.gsl_vector_p());
    260     if (status)
    261       throw utility::GSL_error(std::string("vector::div",status));
    262   }
    263 
    264 
    265   vector::iterator vector::end(void)
    266   {
    267     return vector::iterator(*this, size());
    268   }
    269 
    270 
    271   vector::const_iterator vector::end(void) const
    272   {
    273     return vector::const_iterator(*this, size());
    274   }
    275 
    276 
    277   bool vector::equal(const vector& other, const double d) const
    278   {
    279     if (this==&other)
    280       return true;
    281     if (size()!=other.size())
    282       return false;
    283     // if gsl error handler disabled, out of bounds index will not
    284     // abort the program.
    285     for (size_t i=0; i<size(); ++i)
    286       // The two last condition checks are needed for NaN detection
    287       if (fabs( (*this)(i)-other(i) ) > d ||
    288           (*this)(i)!=(*this)(i) || other(i)!=other(i))
    289         return false;
    290     return true;
    291   }
    292 
    293 
    294   const gsl_vector* vector::gsl_vector_p(void) const
    295   {
    296     return proxy_v_;
    297   }
    298 
    299 
    300   gsl_vector* vector::gsl_vector_p(void)
    301   {
    302     return v_;
    303   }
    304 
    305 
    306   bool vector::isview(void) const
    307   {
    308     return view_ || view_const_;
    309   }
    310 
    311 
    312   void vector::mul(const vector& other)
    313   {
    314     assert(v_);
    315     int status=gsl_vector_mul(v_,other.gsl_vector_p());
    316     if (status)
    317       throw utility::GSL_error(std::string("vector::mul",status));
    318   }
    319 
    320 
    321   void vector::resize(size_t n, double init_value)
    322   {
    323     delete_allocated_memory();
    324     proxy_v_ = v_ = gsl_vector_alloc(n);
    325     if (!v_)
    326       throw utility::GSL_error("vector::vector failed to allocate memory");
    327     all(init_value);
    328   }
    329 
    330 
    331   void vector::reverse(void)
    332   {
    333     assert(v_);
    334     gsl_vector_reverse(v_);
    335   }
    336 
    337 
    338   void vector::all(const double& value)
    339   {
    340     assert(v_);
    341     gsl_vector_set_all(v_,value);
    342   }
    343 
    344 
    345   size_t vector::size(void) const
    346   {
    347     if (!proxy_v_)
    348       return 0;
    349     return proxy_v_->size;
    350   }
    351 
    352 
    353   void vector::swap(size_t i, size_t j)
    354   {
    355     assert(v_);
    356     int status=gsl_vector_swap_elements(v_, i, j);
    357     if (status)
    358       throw utility::GSL_error(std::string("vector::swap_elements",status));
    359   }
    360 
    361 
    362   double& vector::operator()(size_t i)
    363   {
    364     assert(v_);
    365     double* d=gsl_vector_ptr(v_, i);
     117    assert(vec_);
     118    double* d=gsl_vector_ptr(vec_, i);
    366119    if (!d)
    367       throw utility::GSL_error("vector::operator()",GSL_EINVAL);
     120      throw utility::GSL_error("vectorMutable::operator()",GSL_EINVAL);
    368121    return *d;
    369122  }
    370123
    371124
    372   const double& vector::operator()(size_t i) const
     125  const vectorMutable& vectorMutable::operator+=(const vectorBase& other)
    373126  {
    374     const double* d=gsl_vector_const_ptr(proxy_v_, i);
    375     if (!d)
    376       throw utility::GSL_error("vector::operator()",GSL_EINVAL);
    377     return *d;
    378   }
    379 
    380 
    381   double& vector::operator[](size_t i)
    382   {
    383     return this->operator()(i);
    384   }
    385 
    386 
    387   const double& vector::operator[](size_t i) const
    388   {
    389     return this->operator()(i);
    390   }
    391 
    392 
    393   bool vector::operator==(const vector& other) const
    394   {
    395     return equal(other);
    396   }
    397 
    398 
    399   bool vector::operator!=(const vector& other) const
    400   {
    401     return !equal(other);
    402   }
    403 
    404 
    405   double vector::operator*( const vector &other ) const
    406   {
    407     double res = 0.0;;
    408     for ( size_t i = 0; i < size(); ++i )
    409       res += other(i) * (*this)(i);
    410     return res;
    411   }
    412 
    413 
    414   const vector& vector::operator=( const vector& other )
    415   {
    416     assert(v_);
    417     if (gsl_vector_memcpy(v_,other.gsl_vector_p()))
    418       throw utility::GSL_error("vector::set dimension mis-match");
    419     return *this;
    420   }
    421 
    422 
    423   const vector& vector::operator+=(const vector& other)
    424   {
    425     assert(v_);
    426     int status=gsl_vector_add(v_, other.gsl_vector_p());
     127    assert(vec_);
     128    int status=gsl_vector_add(vec_, other.gsl_vector_p());
    427129    if (status)
    428       throw utility::GSL_error(std::string("vector::add", status));
     130      throw utility::GSL_error(std::string("vectorMutable::add", status));
    429131    return *this;
    430132  }
    431133
    432134
    433   const vector& vector::operator+=(double d)
     135  const vectorMutable& vectorMutable::operator+=(double d)
    434136  {
    435     assert(v_);
    436     gsl_vector_add_constant(v_, d);
     137    assert(vec_);
     138    gsl_vector_add_constant(vec_, d);
    437139    return *this;
    438140  }
    439141
    440142
    441   const vector& vector::operator-=(const vector& other)
     143  const vectorMutable& vectorMutable::operator-=(const vectorBase& other)
    442144  {
    443     assert(v_);
    444     int status=gsl_vector_sub(v_, other.gsl_vector_p());
     145    assert(vec_);
     146    int status=gsl_vector_sub(vec_, other.gsl_vector_p());
    445147    if (status)
    446       throw utility::GSL_error(std::string("vector::sub", status));
     148      throw utility::GSL_error(std::string("vectorMutable::sub", status));
    447149    return *this;
    448150  }
    449151
    450152
    451   const vector& vector::operator-=(const double d)
     153  const vectorMutable& vectorMutable::operator-=(const double d)
    452154  {
    453     assert(v_);
    454     gsl_vector_add_constant(v_, -d);
     155    assert(vec_);
     156    gsl_vector_add_constant(vec_, -d);
    455157    return *this;
    456158  }
    457159
    458160
    459   const vector& vector::operator*=(const double d)
     161  const vectorMutable& vectorMutable::operator*=(const double d)
    460162  {
    461     assert(v_);
    462     gsl_vector_scale(v_, d);
     163    assert(vec_);
     164    gsl_vector_scale(vec_, d);
    463165    return *this;
    464166  }
    465167
    466168
    467   bool isnull(const vector& v)
    468   {
    469     return gsl_vector_isnull(v.gsl_vector_p());
    470   }
    471 
    472 
    473   double max(const vector& v)
    474   {
    475     return gsl_vector_max(v.gsl_vector_p());
    476   }
    477 
    478 
    479   size_t max_index(const vector& v)
    480   {
    481     return gsl_vector_max_index(v.gsl_vector_p());
    482   }
    483 
    484 
    485   double min(const vector& v)
    486   {
    487     return gsl_vector_min(v.gsl_vector_p());
    488   }
    489 
    490 
    491   size_t min_index(const vector& v)
    492   {
    493     return gsl_vector_min_index(v.gsl_vector_p());
    494   }
    495 
    496 
    497   bool nan(const vector& templat, vector& flag)
    498   {
    499     size_t vsize=templat.size();
    500     if (vsize!=flag.size())
    501       flag.clone(vector(vsize,1.0));
    502     else
    503       flag.all(1.0);
    504     bool nan=false;
    505     for (size_t i=0; i<vsize; i++)
    506       if (std::isnan(templat(i))) {
    507         flag(i)=0;
    508         nan=true;
    509       }
    510     return nan;
    511   }
    512 
    513 
    514   void set_basis(vector& v, size_t i)
     169  void set_basis(vectorMutable& v, size_t i)
    515170  {
    516171    assert(v.gsl_vector_p());
     
    519174
    520175
    521   void shuffle(vector& invec)
     176  void shuffle(vectorMutable& invec)
    522177  {
    523178    random::DiscreteUniform rnd;
     
    526181
    527182
    528   void sort(vector& v)
     183  void sort(vectorMutable& v)
    529184  {
    530185    assert(v.gsl_vector_p());
     
    532187  }
    533188
    534 
    535   void sort_index(std::vector<size_t>& sort_index, const vector& invec)
     189  vectorMutable::operator proxy ()
    536190  {
    537     assert(invec.gsl_vector_p());
    538     gsl_permutation* p = gsl_permutation_alloc(invec.size());
    539     int status=gsl_sort_vector_index(p,invec.gsl_vector_p());
    540     if (status) {
    541       gsl_permutation_free(p);
    542       throw utility::GSL_error(std::string("sort_index(vector&,const vector&)",status));     
    543     }
    544     sort_index=std::vector<size_t>(p->data,p->data+p->size);
    545     gsl_permutation_free(p);
     191    proxy p;
     192    p.vec_ = this->vec_;
     193    return p;
    546194  }
    547195
    548196
    549   void sort_smallest_index(std::vector<size_t>& sort_index, size_t k,
    550                             const vector& invec)
    551   {
    552     assert(invec.gsl_vector_p());
    553     assert(k<=invec.size());
    554     sort_index.resize(k);
    555     gsl_sort_vector_smallest_index(&sort_index[0],k,invec.gsl_vector_p());
    556   }
    557  
    558   void sort_largest_index(std::vector<size_t>& sort_index, size_t k,
    559                             const vector& invec)
    560   {
    561     assert(invec.gsl_vector_p());
    562     assert(k<=invec.size());
    563     sort_index.resize(k);
    564     gsl_sort_vector_largest_index(&sort_index[0],k,invec.gsl_vector_p());
    565   }
    566 
    567 
    568   double sum(const vector& v)
    569   {
    570     double sum = 0;
    571     size_t vsize=v.size();
    572     for (size_t i=0; i<vsize; ++i)
    573       sum += v[i];
    574     return sum;
    575   }
    576 
    577 
    578   void swap(vector& v, vector& w)
    579   {
    580     assert(v.gsl_vector_p()); assert(w.gsl_vector_p());
    581     int status=gsl_vector_swap(v.gsl_vector_p(),w.gsl_vector_p());
    582     if (status)
    583       throw utility::GSL_error(std::string("swap(vector&,vector&)",status));
    584   }
    585 
    586 
    587   std::ostream& operator<<(std::ostream& s, const vector& a)
    588   {
    589     s.setf(std::ios::dec);
    590     s.precision(12);
    591     for (size_t j = 0; j < a.size(); ++j) {
    592       s << a[j];
    593       if ( (j+1)<a.size() )
    594         s << s.fill();
    595     }
    596     return s;
    597   }
    598197
    599198}}} // of namespace utility, yat, and thep
  • branches/peter-dev/yat/utility/vectorMutable.h

    r994 r995  
    1 #ifndef _theplu_yat_utility_vector_
    2 #define _theplu_yat_utility_vector_
     1#ifndef _theplu_yat_utility_vector_mutable_
     2#define _theplu_yat_utility_vector_mutable_
    33
    44// $Id$
     
    2929*/
    3030
     31#include "vectorBase.h"
     32
    3133#include "Exception.h"
    3234#include "Iterator.h"
     
    4446
    4547  class matrix;
    46 
     48 
    4749  /**
    4850     @brief This is the yat interface to GSL vector.
     
    8486  */
    8587
    86   class vector
     88  class vectorMutable : public vectorBase
    8789  {
    8890  public:
    8991
    90     /// \brief vector::iterator
    91     typedef Iterator<double&, vector> iterator;
    92     /// \brief vector::const_iterator
    93     typedef Iterator<const double, const vector> const_iterator;
     92    /// \brief vectorMutable::iterator
     93    typedef Iterator<double&, vectorMutable> iterator;
    9494
    9595    /**
    9696       \brief The default constructor.
    9797    */
    98     vector(void);
    99 
    100     /**
    101        \brief Allocates memory space for \a n elements, and sets all
    102        elements to \a init_value.
    103        
    104        \throw GSL_error if memory allocation fails.
    105     */
    106     vector(size_t n, double init_value=0);
    107 
    108     /**
    109        \brief The copy constructor.
    110 
    111        \note If the object to be copied is a vector view, the values
    112        of the view will be copied, i.e. the view is not copied.
    113 
    114        \throw A GSL_error is indirectly thrown if memory allocation
    115        fails.
    116     */
    117     vector(const vector& other);
    118 
    119     /**
    120        \brief Vector view constructor.
    121 
    122        Create a view of vector \a v, with starting index \a offset,
    123        size \a n, and an optional \a stride.
    124 
    125        A vector view can be used as any vector with the difference
    126        that changes made to the view will also change the object that
    127        is viewed. Also, using the copy constructor will create a new
    128        vector object that is a copy of whatever is viewed. If a copy
    129        of the view is needed then you should use this constructor to
    130        obtain a copy.
    131 
    132        \note If the object viewed by the view goes out of scope or is
    133        deleted, the view becomes invalid and the result of further use
    134        is undefined.
    135 
    136        \throw GSL_error if a view cannot be set up.
    137     */
    138     vector(vector& v, size_t offset, size_t n, size_t stride=1);
    139 
    140     /**
    141        \brief Vector const view constructor.
    142 
    143        Create a view of vector \a v, with starting index \a offset,
    144        size \a n, and an optional \a stride.
    145 
    146        A vector view can be used as any const vector. Using the copy
    147        constructor will create a new vector object that is a copy of
    148        whatever is viewed. If a copy of the view is needed then you
    149        should use this constructor to obtain a copy.
    150 
    151        \note If the object viewed by the view goes out of scope or is
    152        deleted, the view becomes invalid and the result of further use
    153        is undefined.
    154 
    155        \throw GSL_error if a view cannot be set up.
    156     */
    157     vector(const vector& v, size_t offset, size_t n, size_t stride=1);
    158 
    159     ///
    160     /// Matrix row/column view constructor.
    161     ///
    162     /// Create a row/column vector view of matrix \a m, pointing at
    163     /// row/column \a i. The parameter \a row is used to set whether
    164     /// the view should be a row or column view. If \a row is set to
    165     /// true, the view will be a row view (default behaviour), and,
    166     /// naturally, a column view otherwise.
    167     ///
    168     /// A vector view can be used as any vector with the difference
    169     /// that changes made to the view will also change the object that
    170     /// is viewed. Also, using the copy constructor will create a new
    171     /// vector object that is a copy of whatever is viewed. If a copy
    172     /// of the view is needed then you should use the vector view
    173     /// constructor to obtain a copy.
    174     ///
    175     /// @note If the object viewed by the view goes out of scope or is
    176     /// deleted, the view becomes invalid and the result of further
    177     /// use is undefined.
    178     ///
    179     vector(matrix& m, size_t i, bool row=true);
    180 
    181     ///
    182     /// Matrix row/column const view constructor.
    183     ///
    184     /// Create a row/column vector view of matrix \a m, pointing at
    185     /// row/column \a i. The parameter \a row is used to set whether
    186     /// the view should be a row or column view. If \a row is set to
    187     /// true, the view will be a row view (default behaviour), and,
    188     /// naturally, a column view otherwise.
    189     ///
    190     /// A const vector view can be used as any const vector. Using the
    191     /// copy constructor will create a new vector object that is a
    192     /// copy of whatever is viewed. If a copy of the view is needed
    193     /// then you should use the vector view constructor to obtain a
    194     /// copy.
    195     ///
    196     /// @note If the object viewed by the view goes out of scope or is
    197     /// deleted, the view becomes invalid and the result of further
    198     /// use is undefined.
    199     ///
    200     vector(const matrix& m, size_t i, bool row=true);
    201 
    202     /**
    203        \brief The istream constructor.
    204 
    205        Either elements should be separated with white space characters
    206        (default), or elements should be separated by the delimiter \a
    207        sep. When delimiter \a sep is used empty elements are stored as
    208        NaN's (except that empty lines are ignored). The end of input
    209        to the vector is at end of file marker.
    210 
    211        \throw GSL_error if memory allocation fails, IO_error if
    212        unexpected input is found in the input stream.
    213     */
    214     explicit vector(std::istream &, char sep='\0')
    215       throw(utility::IO_error, std::exception);
     98    vectorMutable(gsl_vector*);
    21699
    217100    ///
    218101    /// The destructor.
    219102    ///
    220     ~vector(void);
    221 
    222     /**
    223        \return mutable iterator to start of vector
     103    virtual ~vectorMutable(void);
     104
     105    // to enable overloading from base class
     106    using vectorBase::begin;
     107
     108    /**
     109       \return mutable iterator to start of vectorMutable
    224110     */
    225111    iterator begin(void);
    226112
    227113    /**
    228        \return read-only iterator to start of vector
    229      */
    230     const_iterator begin(void) const;
    231 
    232     /**
    233        \brief Make a copy of \a other.
    234 
    235        This function will make a deep copy of \a other. Memory is
    236        resized and view state is changed if needed.
    237     */
    238     const vector& clone(const vector& other);
    239 
    240     /**
    241114       \brief This function performs element-wise division, \f$ this_i =
    242115       this_i/other_i \; \forall i \f$.
     
    244117       \throw GSL_error if dimensions mis-match.
    245118    */
    246     void div(const vector& other);
    247 
    248     /**
    249        \return mutable iterator to end of vector
     119    void div(const vectorBase& other);
     120
     121    // to enable overloading from base class
     122    using vectorBase::end;
     123
     124    /**
     125       \return mutable iterator to end of vectorMutable
    250126     */
    251127    iterator end(void);
    252128
    253     /**
    254        \return read-only iterator to end of vector
    255      */
    256     const_iterator end(void) const;
    257 
    258     /**
    259        \brief Check whether vectors are equal within a user defined
    260        precision, set by \a precision.
    261 
    262        \return True if each element deviates less or equal than \a
    263        d. If any vector contain a NaN, false is always returned.
    264 
    265        \see operator== and operator!=
    266     */
    267     bool equal(const vector&, const double precision=0) const;
    268 
    269     ///
    270     /// @return A const pointer to the internal GSL vector,
    271     ///
    272     const gsl_vector* gsl_vector_p(void) const;
     129    // to enable overloading from base class
     130    using vectorBase::gsl_vector_p;
    273131
    274132    ///
     
    276134    ///
    277135    gsl_vector* gsl_vector_p(void);
    278 
    279     ///
    280     /// Check if the vector object is a view (sub-vector) to another
    281     /// vector.
    282     ///
    283     /// @return True if the object is a view, false othwerwise.
    284     ///
    285     bool isview(void) const;
    286136
    287137    /**
     
    291141       \throw GSL_error if dimensions mis-match.
    292142    */
    293     void mul(const vector& other);
    294 
    295     /**
    296        @brief Resize vector
    297        
    298        All elements are set to @a init_value.
    299 
    300        \note Underlying GSL vector is destroyed and a view into this
    301        vector becomes invalid.
    302     */
    303     void resize(size_t, double init_value=0);
    304 
    305     /**
    306        \brief Reverse the order of elements in the vector.
     143    void mul(const vectorBase& other);
     144
     145    /**
     146       \brief Reverse the order of elements in the vectorMutable.
    307147    */
    308148    void reverse(void);
     
    313153    void all(const double& value);
    314154
    315     ///
    316     /// @return the number of elements in the vector.
    317     ///
    318     size_t size(void) const;
    319 
    320155    /**
    321156       \brief Exchange elements \a i and \a j.
    322157
    323        \throw GSL_error if vector lengths differs.
     158       \throw GSL_error if vectorMutable lengths differs.
    324159    */
    325160    void swap(size_t i, size_t j);
     161
     162    // to allow overloading from base class
     163    using vectorBase::operator();
    326164
    327165    /**
     
    337175
    338176    /**
    339        \brief Element access operator.
    340 
    341        \return Const reference to element \a i.
    342 
    343        \throw If GSL range checks are enabled in the underlying GSL
    344        library a GSL_error exception is thrown if either index is out
    345        of range.
    346     */
    347     const double& operator()(size_t i) const;
    348 
    349     ///
    350     /// Element access operator.
    351     ///
    352     /// @return Reference to element \a i.
    353     ///
    354     double& operator[](size_t i);
    355 
    356     ///
    357     /// Const element access operator.
    358     ///
    359     /// @return The value of element \a i.
    360     ///
    361     const double& operator[](size_t i) const;
    362 
    363     /**
    364        \brief Comparison operator. Takes linear time.
    365 
    366        Checks are performed with exact matching, i.e., rounding off
    367        effects may destroy comparison. Use the equal function for
    368        comparing elements within a user defined precision.
    369 
    370        \return True if all elements are equal otherwise false.
    371 
    372        \see equal(const vector&, const double precision=0)
    373     */
    374     bool operator==(const vector&) const;
    375 
    376     /**
    377        \brief Comparison operator. Takes linear time.
    378 
    379        Checks are performed with exact matching, i.e., rounding off
    380        effects may destroy comparison. Use the equal function for
    381        comparing elements within a user defined precision.
    382 
    383        \return False if all elements are equal otherwise true.
    384 
    385        \see equal(const vector&, const double precision=0)
    386     */
    387     bool operator!=(const vector&) const;
    388 
    389     ///
    390     /// @return The dot product.
    391     ///
    392     double operator*(const vector&) const;
    393 
    394     /**
    395177       \brief The assignment operator.
    396178
    397        Dimensions of the vectors must match. If the LHS vector is a
     179       Dimensions of the vectorMutables must match. If the LHS vectorMutable is a
    398180       view, the underlying data will be changed.
    399181
    400        \return A const reference to the resulting vector.
    401 
    402        \see void set(const vector&).
    403 
    404        \throw GSL_error if dimensions mis-match.
    405     */
    406     const vector& operator=(const vector&);
    407 
    408     /**
    409        \brief Addition and assign operator. Vector addition, \f$
     182       \return A const reference to the resulting vectorMutable.
     183
     184       \see void set(const vectorMutable&).
     185
     186       \throw GSL_error if dimensions mis-match.
     187    */
     188    virtual const vectorMutable& operator=(vectorMutable&)=0;
     189
     190    /**
     191       \brief Addition and assign operator. VectorMutable addition, \f$
    410192       this_i = this_i + other_i \; \forall i \f$.
    411193
    412        \return A const reference to the resulting vector.
    413 
    414        \throw GSL_error if dimensions mis-match.
    415     */
    416     const vector& operator+=(const vector&);
    417 
    418     /**
    419        \brief Add a constant to a vector, \f$ this_i = this_i + d \;
     194       \return A const reference to the resulting vectorMutable.
     195
     196       \throw GSL_error if dimensions mis-match.
     197    */
     198    const vectorMutable& operator+=(const vectorBase&);
     199
     200    /**
     201       \brief Add a constant to a vectorMutable, \f$ this_i = this_i + d \;
    420202       \forall i \f$.
    421203
    422        \return A const reference to the resulting vector.
    423     */
    424     const vector& operator+=(double d);
    425 
    426     /**
    427        \brief Subtract and assign operator. Vector subtraction, \f$
     204       \return A const reference to the resulting vectorMutable.
     205    */
     206    const vectorMutable& operator+=(double d);
     207
     208    /**
     209       \brief Subtract and assign operator. VectorMutable subtraction, \f$
    428210       this_i = this_i - other_i \; \forall i \f$.
    429211
    430        \return A const reference to the resulting vector.
    431 
    432        \throw GSL_error if dimensions mis-match.
    433     */
    434     const vector& operator-=(const vector&);
    435 
    436     /**
    437        \brief Subtract a constant to a vector, \f$ this_i = this_i - d
     212       \return A const reference to the resulting vectorMutable.
     213
     214       \throw GSL_error if dimensions mis-match.
     215    */
     216    const vectorMutable& operator-=(const vectorBase&);
     217
     218    /**
     219       \brief Subtract a constant to a vectorMutable, \f$ this_i = this_i - d
    438220       \; \forall i \f$.
    439221
    440        \return A const reference to the resulting vector.
    441     */
    442     const vector& operator-=(double d);
     222       \return A const reference to the resulting vectorMutable.
     223    */
     224    const vectorMutable& operator-=(double d);
    443225
    444226    /**
     
    446228       this_i * d \; \forall i \f$.
    447229
    448        \return A const reference to the resulting vector.
    449     */
    450     const vector& operator*=(double d);
    451 
    452 
    453   private:
    454 
    455     /**
    456        \brief Create a new copy of the internal GSL vector.
    457 
    458        Necessary memory for the new GSL vector is allocated and the
    459        caller is responsible for freeing the allocated memory.
    460 
    461        \return A pointer to a copy of the internal GSL vector.
    462 
    463        \throw GSL_error if memory cannot be allocated for the new
    464        copy, or if dimensions mis-match.
    465     */
    466     gsl_vector* create_gsl_vector_copy(void) const;
    467 
    468     /**
    469        \brief Clear all dynamically allocated memory.
    470 
    471        Internal utility function.
    472     */
    473     void delete_allocated_memory(void);
    474 
    475     gsl_vector* v_;
    476     gsl_vector_view* view_;
    477     const gsl_vector_const_view* view_const_;
    478     // proxy_v_ is used to access the proper underlying gsl_vector
    479     // in all const member functions. It is not used by design for
    480     // non-const vector functions and operators. This is to make sure
    481     // that runtime errors occur if a const vector is used in an
    482     // inappropriate manner such as on left hand side in assignment
    483     // (remember, v_ is null for const vector views).
    484     const gsl_vector* proxy_v_;
     230       \return A const reference to the resulting vectorMutable.
     231    */
     232    const vectorMutable& operator*=(double d);
     233
     234  protected:
     235    gsl_vector* vec_;
     236
     237    struct proxy
     238    {
     239      gsl_vector* vec_;
     240    };
     241
     242  public:
     243    /**
     244     */
     245    operator proxy();
     246
    485247  };
    486248
    487249  /**
    488      \brief Check if all elements of the vector are zero.
    489 
    490      \return True if all elements in the vector is zero, false
    491      othwerwise.
    492   */
    493   bool isnull(const vector&);
    494 
    495   /**
    496      \brief Get the maximum value of the vector.
    497 
    498      \return The maximum value of the vector.
    499   */
    500   double max(const vector&);
    501 
    502   /**
    503      \brief Locate the maximum value in the vector.
    504 
    505      \return The index to the maximum value of the vector.
    506 
    507      \note Lower index has precedence.
    508   */
    509   size_t max_index(const vector&);
    510 
    511   /**
    512      \brief Get the minimum value of the vector.
    513 
    514      \return The minimum value of the vector.
    515   */
    516   double min(const vector&);
    517 
    518   /**
    519      \brief Locate the minimum value in the vector.
    520 
    521      \return The index to the minimum value of the vector.
    522 
    523      \note Lower index has precedence.
    524   */
    525   size_t min_index(const vector&);
    526 
    527   /**
    528      \brief Create a vector \a flag indicating NaN's in another vector
    529      \a templat.
    530 
    531      The \a flag vector is changed to contain 1's and 0's only. A 1
    532      means that the corresponding element in the \a templat vector is
    533      valid and a zero means that the corresponding element is a NaN.
    534 
    535      \note Space for vector \a flag is reallocated to fit the size of
    536      vector \a templat if sizes mismatch.
    537 
    538      \return True if the \a templat vector contains at least one NaN.
    539   */
    540   bool nan(const vector& templat, vector& flag);
    541 
    542   /**
    543      \brief Transforms a vector to a basis vector.
     250     \brief Transforms a vectorMutable to a basis vector.
    544251
    545252     All elements are set to zero except the \a i-th element which is
    546253     set to one.
    547254  */
    548   void set_basis(vector&, size_t i);
    549 
    550   /**
    551      Randomly shuffles the elements in vector \a invec
    552   */
    553   void shuffle(vector& invec);
    554 
    555   /**
    556      Sort the elements in the vector.
    557   */
    558   void sort(vector&);
    559 
    560   /**
    561      Create a vector \a sort_index containing the indeces of
    562      elements in a another vector \a invec.  The elements of \a
    563      sort_index give the index of the vector element which would
    564      have been stored in that position if the vector had been sorted
    565      in place. The first element of \a sort_index gives the index of the least
    566      element in \a invec, and the last element of \a sort_index gives the index of the
    567      greatest element in \a invec . The vector \a invec is not changed.
    568 
    569   */
    570   void sort_index(std::vector<size_t>& sort_index, const vector& invec);
    571 
    572   /** Similar to sort_index but creates a vector with indices to the \a k
    573   smallest elements in \a invec. 
    574   */
    575   void sort_smallest_index(std::vector<size_t>& sort_index, size_t k, const
    576   vector& invec);
    577 
    578   /** Similar to sort_index but creates a vector with indices to the \a k
    579   largest elements in \a invec. 
    580   */
    581   void sort_largest_index(std::vector<size_t>& sort_index, size_t k, const
    582   vector& invec);
    583 
    584  
    585 
    586   /**
    587      \brief Calculate the sum of all vector elements.
    588 
    589      \return The sum.
    590   */
    591   double sum(const vector&);
     255  void set_basis(vectorMutable&, size_t i);
     256
     257  /**
     258     Randomly shuffles the elements in vectorMutable \a invec
     259  */
     260  void shuffle(vectorMutable& invec);
     261
     262  /**
     263     Sort the elements in the vectorMutable.
     264  */
     265  void sort(vectorMutable&);
    592266
    593267  /**
    594268     \brief Swap vector elements by copying.
    595269
    596      The two vectors must have the same length.
    597 
    598      \throw GSL_error if vector lengths differs.
    599   */
    600   void swap(vector&, vector&);
    601 
    602   /**
    603      \brief The output operator for the vector class.
    604   */
    605   std::ostream& operator<<(std::ostream&, const vector&);
     270     The two vectorMutables must have the same length.
     271
     272     \throw GSL_error if vectorMutable lengths differs.
     273  */
     274  //void swap(vectorMutable&, vectorMutable&);
    606275
    607276}}} // of namespace utility, yat, and theplu
  • branches/peter-dev/yat/utility/vectorView.cc

    r994 r995  
    2525*/
    2626
    27 #include "vector.h"
     27#include "vectorView.h"
    2828#include "matrix.h"
    2929#include "utility.h"
     
    4343
    4444
    45   vector::vector(void)
    46     : v_(NULL), view_(NULL), view_const_(NULL), proxy_v_(NULL)
     45  vectorView::vectorView(void)
     46    : vectorMutable(NULL)
    4747  {
    4848  }
    4949
    5050
    51   vector::vector(size_t n, double init_value)
    52     : v_(gsl_vector_alloc(n)), view_(NULL), view_const_(NULL),
    53       proxy_v_(v_)
     51  vectorView::vectorView(vectorMutable& other)
     52    : vectorMutable(NULL), view_(NULL)
    5453  {
    55     if (!v_)
    56       throw utility::GSL_error("vector::vector failed to allocate memory");
    57 
    58     all(init_value);
     54    view_ = new gsl_vector_view(gsl_vector_subvector(other.gsl_vector_p(),0,
     55                                                     other.size()));
     56    const_vec_ = vec_ = &(view_->vector);
    5957  }
    6058
    6159
    62   vector::vector(const vector& other)
    63     : v_(other.create_gsl_vector_copy()), view_(NULL),
    64       view_const_(NULL), proxy_v_(v_)
     60  vectorView::vectorView(vectorMutable& v, size_t offset,size_t n,size_t stride)
     61    : vectorMutable(NULL)
    6562  {
     63    view_ =
     64      new gsl_vector_view(gsl_vector_subvector_with_stride(v.gsl_vector_p(),
     65                                                           offset, stride,n));
     66    const_vec_ = vec_ = &(view_->vector);
    6667  }
    6768
    6869
    69   vector::vector(vector& v, size_t offset, size_t n, size_t stride)
    70     : view_const_(NULL)
    71   {
    72     view_ = new gsl_vector_view(gsl_vector_subvector_with_stride(v.v_,offset,
    73                                                                  stride,n));
    74     if (!view_)
    75       throw utility::GSL_error("vector::vector failed to setup view");
    76     proxy_v_ = v_ = &(view_->vector);
    77   }
    78 
    79 
    80   vector::vector(const vector& v, size_t offset, size_t n, size_t stride)
    81     : v_(NULL), view_(NULL)
    82   {
    83     view_const_ = new gsl_vector_const_view(
    84                    gsl_vector_const_subvector_with_stride(v.v_,offset,stride,n));
    85     if (!view_const_)
    86       throw utility::GSL_error("vector::vector failed to setup view");
    87     proxy_v_ = &(view_const_->vector);
    88   }
    89 
    90 
    91   vector::vector(matrix& m, size_t i, bool row)
    92     : view_const_(NULL)
     70  vectorView::vectorView(matrix& m, size_t i, bool row)
     71    : vectorMutable(NULL)
    9372  {
    9473    view_=new gsl_vector_view(row ?
    9574                              gsl_matrix_row   (m.gsl_matrix_p(),i) :
    9675                              gsl_matrix_column(m.gsl_matrix_p(),i)  );
    97     if (!view_)
    98       throw utility::GSL_error("vector::vector failed to setup view");
    99     proxy_v_ = v_ = &(view_->vector);
     76    const_vec_ = vec_ = &(view_->vector);
    10077  }
    10178
    10279
    103   vector::vector(const matrix& m, size_t i, bool row)
    104     : v_(NULL), view_(NULL)
     80  vectorView::vectorView(proxy p)
     81    : vectorMutable(NULL)
    10582  {
    106     view_const_ = new gsl_vector_const_view(row ?
    107                                    gsl_matrix_const_row   (m.gsl_matrix_p(),i) :
    108                                    gsl_matrix_const_column(m.gsl_matrix_p(),i) );
    109     if (!view_const_)
    110       throw utility::GSL_error("vector::vector failed to setup view");
    111     proxy_v_ = &(view_const_->vector);
     83    view_ = new gsl_vector_view(gsl_vector_subvector(p.vec_, 0, p.vec_->size));
     84    const_vec_ = vec_ = &(view_->vector);
    11285  }
    11386
    11487
    115   vector::vector(std::istream& is, char sep)
    116     throw (utility::IO_error, std::exception)
    117     : view_(NULL), view_const_(NULL)
     88  vectorView::~vectorView(void)
    11889  {
    119     // read the data file and store in stl vectors (dynamically
    120     // expandable)
    121     std::vector<std::vector<double> > data_matrix;
    122     u_int nof_columns=0;
    123     u_int nof_rows=0;
    124     std::string line;
    125     while(getline(is, line, '\n')) {
    126       // Empty lines
    127       if (!line.size())
    128           continue;
    129       nof_rows++;
    130 
    131       std::vector<double> v;
    132       std::string element;
    133       std::stringstream ss(line);
    134       bool ok=true;
    135       while(ok) {
    136         if(sep=='\0')
    137           ok=(ss>>element);
    138         else
    139           ok=getline(ss, element, sep);
    140         if(!ok)
    141           break;       
    142 
    143         if(utility::is_double(element)) {
    144           v.push_back(atof(element.c_str()));
    145         }
    146         else if (!element.size() || utility::is_nan(element)) {
    147           v.push_back(std::numeric_limits<double>::quiet_NaN());
    148         }
    149         else {
    150           std::stringstream ss("Warning: '");
    151           ss << element << "' is not a double.";
    152           throw IO_error(ss.str());
    153         }
    154       }
    155       if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
    156           v.push_back(std::numeric_limits<double>::quiet_NaN());
    157       if (!nof_columns)
    158         nof_columns=v.size();
    159       else if (nof_rows && (nof_columns>1)) {
    160         std::ostringstream s;
    161         s << "vector::vector(std::istream&) data file error:\n"
    162           << "    File has inconsistent number of rows (" << nof_rows
    163           << ") and columns (" << nof_columns
    164           << ").\n    Expected a row or a column vector.";
    165         throw utility::IO_error(s.str());
    166       }
    167       else if (v.size()!=nof_columns) {
    168         std::ostringstream s;
    169         s << "vector::vector(std::istream&) data file error:\n"
    170           << "    Line " << nof_rows << " has " << v.size()
    171           << " columns; expected " << nof_columns << " column.";
    172         throw utility::IO_error(s.str());
    173       }
    174       data_matrix.push_back(v);
    175     }
    176 
    177     // manipulate the state of the stream to be good
    178     is.clear(std::ios::goodbit);
    179     // convert the data to a gsl vector
    180     proxy_v_ = v_ = gsl_vector_alloc(nof_rows*nof_columns);
    181     if (!v_)
    182       throw utility::GSL_error("vector::vector failed to allocate memory");
    183     size_t n=0;
    184     // if gsl error handler disabled, out of bounds index will not
    185     // abort the program.
    186     for (size_t i=0; i<nof_rows; i++)
    187       for (size_t j=0; j<nof_columns; j++)
    188         gsl_vector_set( v_, n++, data_matrix[i][j] );
     90    clean_up();
    18991  }
    19092
    19193
    192   vector::~vector(void)
     94  const vectorMutable& vectorView::operator=(vectorMutable& other )
    19395  {
    194     delete_allocated_memory();
    195   }
    196 
    197 
    198   vector::iterator vector::begin(void)
    199   {
    200     return vector::iterator(*this, 0);
    201   }
    202 
    203 
    204   vector::const_iterator vector::begin(void) const
    205   {
    206     return vector::const_iterator(*this, 0);
    207   }
    208 
    209 
    210   const vector& vector::clone(const vector& other)
    211   {
    212     if (this!=&other) {
    213 
    214       delete_allocated_memory();
    215 
    216       if (other.view_) {
    217         view_ = new gsl_vector_view(*other.view_);
    218         proxy_v_ = v_ = &(view_->vector);
    219       }
    220       else if (other.view_const_) {
    221         view_const_ = new gsl_vector_const_view(*other.view_const_);
    222         proxy_v_ = &(view_const_->vector);
    223         v_=NULL;
    224       }
    225       else if (other.v_)
    226         proxy_v_ = v_ = other.create_gsl_vector_copy();
    227 
    228     }
     96    if (size()!=other.size())
     97      throw utility::GSL_error("vectorView::dimension mis-match");
     98    if (!other.size())
     99      return *this;
     100    clean_up();
     101    view_ = new gsl_vector_view(gsl_vector_subvector(other.gsl_vector_p(),0,
     102                                                     other.size()));
     103    const_vec_ = vec_ = &view_->vector;
    229104    return *this;
    230105  }
    231106
    232107
    233   gsl_vector* vector::create_gsl_vector_copy(void) const
     108  const vectorMutable& vectorView::operator=(proxy p)
    234109  {
    235     gsl_vector* vec = gsl_vector_alloc(size());
    236     if (!vec)
    237       throw utility::GSL_error("vector::create_gsl_vector_copy failed to allocate memory");
    238     if (gsl_vector_memcpy(vec, proxy_v_))
    239       throw utility::GSL_error("vector::create_gsl_matrix_copy dimension mis-match");
    240     return vec;
    241   }
    242 
    243 
    244   void vector::delete_allocated_memory(void)
    245   {
    246     if (view_)
    247       delete view_;
    248     else if (view_const_)
    249       delete view_const_;
    250     else if (v_)
    251       gsl_vector_free(v_);
    252     proxy_v_=v_=NULL;
    253   }
    254 
    255 
    256   void vector::div(const vector& other)
    257   {
    258     assert(v_);
    259     int status=gsl_vector_div(v_,other.gsl_vector_p());
    260     if (status)
    261       throw utility::GSL_error(std::string("vector::div",status));
    262   }
    263 
    264 
    265   vector::iterator vector::end(void)
    266   {
    267     return vector::iterator(*this, size());
    268   }
    269 
    270 
    271   vector::const_iterator vector::end(void) const
    272   {
    273     return vector::const_iterator(*this, size());
    274   }
    275 
    276 
    277   bool vector::equal(const vector& other, const double d) const
    278   {
    279     if (this==&other)
    280       return true;
    281     if (size()!=other.size())
    282       return false;
    283     // if gsl error handler disabled, out of bounds index will not
    284     // abort the program.
    285     for (size_t i=0; i<size(); ++i)
    286       // The two last condition checks are needed for NaN detection
    287       if (fabs( (*this)(i)-other(i) ) > d ||
    288           (*this)(i)!=(*this)(i) || other(i)!=other(i))
    289         return false;
    290     return true;
    291   }
    292 
    293 
    294   const gsl_vector* vector::gsl_vector_p(void) const
    295   {
    296     return proxy_v_;
    297   }
    298 
    299 
    300   gsl_vector* vector::gsl_vector_p(void)
    301   {
    302     return v_;
    303   }
    304 
    305 
    306   bool vector::isview(void) const
    307   {
    308     return view_ || view_const_;
    309   }
    310 
    311 
    312   void vector::mul(const vector& other)
    313   {
    314     assert(v_);
    315     int status=gsl_vector_mul(v_,other.gsl_vector_p());
    316     if (status)
    317       throw utility::GSL_error(std::string("vector::mul",status));
    318   }
    319 
    320 
    321   void vector::resize(size_t n, double init_value)
    322   {
    323     delete_allocated_memory();
    324     proxy_v_ = v_ = gsl_vector_alloc(n);
    325     if (!v_)
    326       throw utility::GSL_error("vector::vector failed to allocate memory");
    327     all(init_value);
    328   }
    329 
    330 
    331   void vector::reverse(void)
    332   {
    333     assert(v_);
    334     gsl_vector_reverse(v_);
    335   }
    336 
    337 
    338   void vector::all(const double& value)
    339   {
    340     assert(v_);
    341     gsl_vector_set_all(v_,value);
    342   }
    343 
    344 
    345   size_t vector::size(void) const
    346   {
    347     if (!proxy_v_)
    348       return 0;
    349     return proxy_v_->size;
    350   }
    351 
    352 
    353   void vector::swap(size_t i, size_t j)
    354   {
    355     assert(v_);
    356     int status=gsl_vector_swap_elements(v_, i, j);
    357     if (status)
    358       throw utility::GSL_error(std::string("vector::swap_elements",status));
    359   }
    360 
    361 
    362   double& vector::operator()(size_t i)
    363   {
    364     assert(v_);
    365     double* d=gsl_vector_ptr(v_, i);
    366     if (!d)
    367       throw utility::GSL_error("vector::operator()",GSL_EINVAL);
    368     return *d;
    369   }
    370 
    371 
    372   const double& vector::operator()(size_t i) const
    373   {
    374     const double* d=gsl_vector_const_ptr(proxy_v_, i);
    375     if (!d)
    376       throw utility::GSL_error("vector::operator()",GSL_EINVAL);
    377     return *d;
    378   }
    379 
    380 
    381   double& vector::operator[](size_t i)
    382   {
    383     return this->operator()(i);
    384   }
    385 
    386 
    387   const double& vector::operator[](size_t i) const
    388   {
    389     return this->operator()(i);
    390   }
    391 
    392 
    393   bool vector::operator==(const vector& other) const
    394   {
    395     return equal(other);
    396   }
    397 
    398 
    399   bool vector::operator!=(const vector& other) const
    400   {
    401     return !equal(other);
    402   }
    403 
    404 
    405   double vector::operator*( const vector &other ) const
    406   {
    407     double res = 0.0;;
    408     for ( size_t i = 0; i < size(); ++i )
    409       res += other(i) * (*this)(i);
    410     return res;
    411   }
    412 
    413 
    414   const vector& vector::operator=( const vector& other )
    415   {
    416     assert(v_);
    417     if (gsl_vector_memcpy(v_,other.gsl_vector_p()))
    418       throw utility::GSL_error("vector::set dimension mis-match");
     110    if (size()!=p.vec_->size)
     111      throw utility::GSL_error("vectorView::dimension mis-match");
     112    if (!size())
     113      return *this;
     114    clean_up();
     115    view_ = new gsl_vector_view(gsl_vector_subvector(p.vec_,0, p.vec_->size));
     116    const_vec_ = vec_ = &view_->vector;
    419117    return *this;
    420118  }
    421119
    422120
    423   const vector& vector::operator+=(const vector& other)
     121  void vectorView::clean_up(void)
    424122  {
    425     assert(v_);
    426     int status=gsl_vector_add(v_, other.gsl_vector_p());
    427     if (status)
    428       throw utility::GSL_error(std::string("vector::add", status));
    429     return *this;
     123    if (view_)
     124      delete view_;
    430125  }
    431126
    432127
    433   const vector& vector::operator+=(double d)
    434   {
    435     assert(v_);
    436     gsl_vector_add_constant(v_, d);
    437     return *this;
    438   }
    439 
    440 
    441   const vector& vector::operator-=(const vector& other)
    442   {
    443     assert(v_);
    444     int status=gsl_vector_sub(v_, other.gsl_vector_p());
    445     if (status)
    446       throw utility::GSL_error(std::string("vector::sub", status));
    447     return *this;
    448   }
    449 
    450 
    451   const vector& vector::operator-=(const double d)
    452   {
    453     assert(v_);
    454     gsl_vector_add_constant(v_, -d);
    455     return *this;
    456   }
    457 
    458 
    459   const vector& vector::operator*=(const double d)
    460   {
    461     assert(v_);
    462     gsl_vector_scale(v_, d);
    463     return *this;
    464   }
    465 
    466 
    467   bool isnull(const vector& v)
    468   {
    469     return gsl_vector_isnull(v.gsl_vector_p());
    470   }
    471 
    472 
    473   double max(const vector& v)
    474   {
    475     return gsl_vector_max(v.gsl_vector_p());
    476   }
    477 
    478 
    479   size_t max_index(const vector& v)
    480   {
    481     return gsl_vector_max_index(v.gsl_vector_p());
    482   }
    483 
    484 
    485   double min(const vector& v)
    486   {
    487     return gsl_vector_min(v.gsl_vector_p());
    488   }
    489 
    490 
    491   size_t min_index(const vector& v)
    492   {
    493     return gsl_vector_min_index(v.gsl_vector_p());
    494   }
    495 
    496 
    497   bool nan(const vector& templat, vector& flag)
    498   {
    499     size_t vsize=templat.size();
    500     if (vsize!=flag.size())
    501       flag.clone(vector(vsize,1.0));
    502     else
    503       flag.all(1.0);
    504     bool nan=false;
    505     for (size_t i=0; i<vsize; i++)
    506       if (std::isnan(templat(i))) {
    507         flag(i)=0;
    508         nan=true;
    509       }
    510     return nan;
    511   }
    512 
    513 
    514   void set_basis(vector& v, size_t i)
    515   {
    516     assert(v.gsl_vector_p());
    517     gsl_vector_set_basis(v.gsl_vector_p(),i);
    518   }
    519 
    520 
    521   void shuffle(vector& invec)
    522   {
    523     random::DiscreteUniform rnd;
    524     std::random_shuffle(invec.begin(), invec.end(), rnd);
    525   }
    526 
    527 
    528   void sort(vector& v)
    529   {
    530     assert(v.gsl_vector_p());
    531     std::sort(v.begin(), v.end());
    532   }
    533 
    534 
    535   void sort_index(std::vector<size_t>& sort_index, const vector& invec)
    536   {
    537     assert(invec.gsl_vector_p());
    538     gsl_permutation* p = gsl_permutation_alloc(invec.size());
    539     int status=gsl_sort_vector_index(p,invec.gsl_vector_p());
    540     if (status) {
    541       gsl_permutation_free(p);
    542       throw utility::GSL_error(std::string("sort_index(vector&,const vector&)",status));     
    543     }
    544     sort_index=std::vector<size_t>(p->data,p->data+p->size);
    545     gsl_permutation_free(p);
    546   }
    547 
    548 
    549   void sort_smallest_index(std::vector<size_t>& sort_index, size_t k,
    550                             const vector& invec)
    551   {
    552     assert(invec.gsl_vector_p());
    553     assert(k<=invec.size());
    554     sort_index.resize(k);
    555     gsl_sort_vector_smallest_index(&sort_index[0],k,invec.gsl_vector_p());
    556   }
    557  
    558   void sort_largest_index(std::vector<size_t>& sort_index, size_t k,
    559                             const vector& invec)
    560   {
    561     assert(invec.gsl_vector_p());
    562     assert(k<=invec.size());
    563     sort_index.resize(k);
    564     gsl_sort_vector_largest_index(&sort_index[0],k,invec.gsl_vector_p());
    565   }
    566 
    567 
    568   double sum(const vector& v)
    569   {
    570     double sum = 0;
    571     size_t vsize=v.size();
    572     for (size_t i=0; i<vsize; ++i)
    573       sum += v[i];
    574     return sum;
    575   }
    576 
    577 
    578   void swap(vector& v, vector& w)
    579   {
    580     assert(v.gsl_vector_p()); assert(w.gsl_vector_p());
    581     int status=gsl_vector_swap(v.gsl_vector_p(),w.gsl_vector_p());
    582     if (status)
    583       throw utility::GSL_error(std::string("swap(vector&,vector&)",status));
    584   }
    585 
    586 
    587   std::ostream& operator<<(std::ostream& s, const vector& a)
    588   {
    589     s.setf(std::ios::dec);
    590     s.precision(12);
    591     for (size_t j = 0; j < a.size(); ++j) {
    592       s << a[j];
    593       if ( (j+1)<a.size() )
    594         s << s.fill();
    595     }
    596     return s;
    597   }
    598 
    599128}}} // of namespace utility, yat, and thep
  • branches/peter-dev/yat/utility/vectorView.h

    r994 r995  
    1 #ifndef _theplu_yat_utility_vector_
    2 #define _theplu_yat_utility_vector_
     1#ifndef _theplu_yat_utility_vector_view_
     2#define _theplu_yat_utility_vector_view_
    33
    44// $Id$
     
    2929*/
    3030
     31#include "vectorMutable.h"
     32
    3133#include "Exception.h"
    3234#include "Iterator.h"
     
    4446
    4547  class matrix;
     48  class vector;
    4649
    4750  /**
     
    8487  */
    8588
    86   class vector
     89  class vectorView : public vectorMutable
    8790  {
    8891  public:
    89 
    90     /// \brief vector::iterator
    91     typedef Iterator<double&, vector> iterator;
    92     /// \brief vector::const_iterator
    93     typedef Iterator<const double, const vector> const_iterator;
     92    /// \brief vectorView::iterator
     93    typedef Iterator<double&, vectorView> iterator;
    9494
    9595    /**
    96        \brief The default constructor.
     96       \brief Default constructor.
    9797    */
    98     vector(void);
    99 
    100     /**
    101        \brief Allocates memory space for \a n elements, and sets all
    102        elements to \a init_value.
    103        
    104        \throw GSL_error if memory allocation fails.
    105     */
    106     vector(size_t n, double init_value=0);
     98    vectorView(void);
    10799
    108100    /**
    109101       \brief The copy constructor.
    110 
    111        \note If the object to be copied is a vector view, the values
    112        of the view will be copied, i.e. the view is not copied.
    113 
    114        \throw A GSL_error is indirectly thrown if memory allocation
    115        fails.
    116102    */
    117     vector(const vector& other);
     103    vectorView(vectorMutable& other);
    118104
    119105    /**
    120        \brief Vector view constructor.
     106       \brief VectorView constructor.
    121107
    122        Create a view of vector \a v, with starting index \a offset,
     108       Create a view of vectorView \a v, with starting index \a offset,
    123109       size \a n, and an optional \a stride.
    124110
    125        A vector view can be used as any vector with the difference
     111       A vectorView view can be used as any vectorView with the difference
    126112       that changes made to the view will also change the object that
    127113       is viewed. Also, using the copy constructor will create a new
    128        vector object that is a copy of whatever is viewed. If a copy
     114       vectorView object that is a copy of whatever is viewed. If a copy
    129115       of the view is needed then you should use this constructor to
    130116       obtain a copy.
     
    136122       \throw GSL_error if a view cannot be set up.
    137123    */
    138     vector(vector& v, size_t offset, size_t n, size_t stride=1);
    139 
    140     /**
    141        \brief Vector const view constructor.
    142 
    143        Create a view of vector \a v, with starting index \a offset,
    144        size \a n, and an optional \a stride.
    145 
    146        A vector view can be used as any const vector. Using the copy
    147        constructor will create a new vector object that is a copy of
    148        whatever is viewed. If a copy of the view is needed then you
    149        should use this constructor to obtain a copy.
    150 
    151        \note If the object viewed by the view goes out of scope or is
    152        deleted, the view becomes invalid and the result of further use
    153        is undefined.
    154 
    155        \throw GSL_error if a view cannot be set up.
    156     */
    157     vector(const vector& v, size_t offset, size_t n, size_t stride=1);
     124    vectorView(vectorMutable& v, size_t offset, size_t n, size_t stride=1);
    158125
    159126    ///
    160127    /// Matrix row/column view constructor.
    161128    ///
    162     /// Create a row/column vector view of matrix \a m, pointing at
     129    /// Create a row/column vectorView view of matrix \a m, pointing at
    163130    /// row/column \a i. The parameter \a row is used to set whether
    164131    /// the view should be a row or column view. If \a row is set to
     
    166133    /// naturally, a column view otherwise.
    167134    ///
    168     /// A vector view can be used as any vector with the difference
     135    /// A vectorView view can be used as any vectorView with the difference
    169136    /// that changes made to the view will also change the object that
    170137    /// is viewed. Also, using the copy constructor will create a new
    171     /// vector object that is a copy of whatever is viewed. If a copy
    172     /// of the view is needed then you should use the vector view
     138    /// vectorView object that is a copy of whatever is viewed. If a copy
     139    /// of the view is needed then you should use the vectorView view
    173140    /// constructor to obtain a copy.
    174141    ///
     
    177144    /// use is undefined.
    178145    ///
    179     vector(matrix& m, size_t i, bool row=true);
    180 
    181     ///
    182     /// Matrix row/column const view constructor.
    183     ///
    184     /// Create a row/column vector view of matrix \a m, pointing at
    185     /// row/column \a i. The parameter \a row is used to set whether
    186     /// the view should be a row or column view. If \a row is set to
    187     /// true, the view will be a row view (default behaviour), and,
    188     /// naturally, a column view otherwise.
    189     ///
    190     /// A const vector view can be used as any const vector. Using the
    191     /// copy constructor will create a new vector object that is a
    192     /// copy of whatever is viewed. If a copy of the view is needed
    193     /// then you should use the vector view constructor to obtain a
    194     /// copy.
    195     ///
    196     /// @note If the object viewed by the view goes out of scope or is
    197     /// deleted, the view becomes invalid and the result of further
    198     /// use is undefined.
    199     ///
    200     vector(const matrix& m, size_t i, bool row=true);
     146    vectorView(matrix& m, size_t i, bool row=true);
    201147
    202148    /**
    203        \brief The istream constructor.
    204 
    205        Either elements should be separated with white space characters
    206        (default), or elements should be separated by the delimiter \a
    207        sep. When delimiter \a sep is used empty elements are stored as
    208        NaN's (except that empty lines are ignored). The end of input
    209        to the vector is at end of file marker.
    210 
    211        \throw GSL_error if memory allocation fails, IO_error if
    212        unexpected input is found in the input stream.
    213     */
    214     explicit vector(std::istream &, char sep='\0')
    215       throw(utility::IO_error, std::exception);
     149     */
     150    vectorView(proxy p);
    216151
    217152    ///
    218153    /// The destructor.
    219154    ///
    220     ~vector(void);
    221 
    222     /**
    223        \return mutable iterator to start of vector
    224      */
    225     iterator begin(void);
    226 
    227     /**
    228        \return read-only iterator to start of vector
    229      */
    230     const_iterator begin(void) const;
    231 
    232     /**
    233        \brief Make a copy of \a other.
    234 
    235        This function will make a deep copy of \a other. Memory is
    236        resized and view state is changed if needed.
    237     */
    238     const vector& clone(const vector& other);
    239 
    240     /**
    241        \brief This function performs element-wise division, \f$ this_i =
    242        this_i/other_i \; \forall i \f$.
    243 
    244        \throw GSL_error if dimensions mis-match.
    245     */
    246     void div(const vector& other);
    247 
    248     /**
    249        \return mutable iterator to end of vector
    250      */
    251     iterator end(void);
    252 
    253     /**
    254        \return read-only iterator to end of vector
    255      */
    256     const_iterator end(void) const;
    257 
    258     /**
    259        \brief Check whether vectors are equal within a user defined
    260        precision, set by \a precision.
    261 
    262        \return True if each element deviates less or equal than \a
    263        d. If any vector contain a NaN, false is always returned.
    264 
    265        \see operator== and operator!=
    266     */
    267     bool equal(const vector&, const double precision=0) const;
    268 
    269     ///
    270     /// @return A const pointer to the internal GSL vector,
    271     ///
    272     const gsl_vector* gsl_vector_p(void) const;
    273 
    274     ///
    275     /// @return A pointer to the internal GSL vector,
    276     ///
    277     gsl_vector* gsl_vector_p(void);
    278 
    279     ///
    280     /// Check if the vector object is a view (sub-vector) to another
    281     /// vector.
    282     ///
    283     /// @return True if the object is a view, false othwerwise.
    284     ///
    285     bool isview(void) const;
    286 
    287     /**
    288        \brief This function performs element-wise multiplication, \f$
    289        this_i = this_i * other_i \; \forall i \f$.
    290 
    291        \throw GSL_error if dimensions mis-match.
    292     */
    293     void mul(const vector& other);
    294 
    295     /**
    296        @brief Resize vector
    297        
    298        All elements are set to @a init_value.
    299 
    300        \note Underlying GSL vector is destroyed and a view into this
    301        vector becomes invalid.
    302     */
    303     void resize(size_t, double init_value=0);
    304 
    305     /**
    306        \brief Reverse the order of elements in the vector.
    307     */
    308     void reverse(void);
    309 
    310     ///
    311     /// Set all elements to \a value.
    312     ///
    313     void all(const double& value);
    314 
    315     ///
    316     /// @return the number of elements in the vector.
    317     ///
    318     size_t size(void) const;
    319 
    320     /**
    321        \brief Exchange elements \a i and \a j.
    322 
    323        \throw GSL_error if vector lengths differs.
    324     */
    325     void swap(size_t i, size_t j);
    326 
    327     /**
    328        \brief Element access operator.
    329 
    330        \return Reference to element \a i.
    331 
    332        \throw If GSL range checks are enabled in the underlying GSL
    333        library a GSL_error exception is thrown if either index is out
    334        of range.
    335     */
    336     double& operator()(size_t i);
    337 
    338     /**
    339        \brief Element access operator.
    340 
    341        \return Const reference to element \a i.
    342 
    343        \throw If GSL range checks are enabled in the underlying GSL
    344        library a GSL_error exception is thrown if either index is out
    345        of range.
    346     */
    347     const double& operator()(size_t i) const;
    348 
    349     ///
    350     /// Element access operator.
    351     ///
    352     /// @return Reference to element \a i.
    353     ///
    354     double& operator[](size_t i);
    355 
    356     ///
    357     /// Const element access operator.
    358     ///
    359     /// @return The value of element \a i.
    360     ///
    361     const double& operator[](size_t i) const;
    362 
    363     /**
    364        \brief Comparison operator. Takes linear time.
    365 
    366        Checks are performed with exact matching, i.e., rounding off
    367        effects may destroy comparison. Use the equal function for
    368        comparing elements within a user defined precision.
    369 
    370        \return True if all elements are equal otherwise false.
    371 
    372        \see equal(const vector&, const double precision=0)
    373     */
    374     bool operator==(const vector&) const;
    375 
    376     /**
    377        \brief Comparison operator. Takes linear time.
    378 
    379        Checks are performed with exact matching, i.e., rounding off
    380        effects may destroy comparison. Use the equal function for
    381        comparing elements within a user defined precision.
    382 
    383        \return False if all elements are equal otherwise true.
    384 
    385        \see equal(const vector&, const double precision=0)
    386     */
    387     bool operator!=(const vector&) const;
    388 
    389     ///
    390     /// @return The dot product.
    391     ///
    392     double operator*(const vector&) const;
     155    ~vectorView(void);
    393156
    394157    /**
    395158       \brief The assignment operator.
    396 
    397        Dimensions of the vectors must match. If the LHS vector is a
    398        view, the underlying data will be changed.
    399 
    400        \return A const reference to the resulting vector.
    401 
    402        \see void set(const vector&).
    403 
    404        \throw GSL_error if dimensions mis-match.
    405     */
    406     const vector& operator=(const vector&);
    407 
    408     /**
    409        \brief Addition and assign operator. Vector addition, \f$
    410        this_i = this_i + other_i \; \forall i \f$.
    411159
    412160       \return A const reference to the resulting vector.
     
    414162       \throw GSL_error if dimensions mis-match.
    415163    */
    416     const vector& operator+=(const vector&);
     164    const vectorMutable& operator=(vectorMutable&);
    417165
    418     /**
    419        \brief Add a constant to a vector, \f$ this_i = this_i + d \;
    420        \forall i \f$.
    421 
    422        \return A const reference to the resulting vector.
    423     */
    424     const vector& operator+=(double d);
    425 
    426     /**
    427        \brief Subtract and assign operator. Vector subtraction, \f$
    428        this_i = this_i - other_i \; \forall i \f$.
    429 
    430        \return A const reference to the resulting vector.
    431 
    432        \throw GSL_error if dimensions mis-match.
    433     */
    434     const vector& operator-=(const vector&);
    435 
    436     /**
    437        \brief Subtract a constant to a vector, \f$ this_i = this_i - d
    438        \; \forall i \f$.
    439 
    440        \return A const reference to the resulting vector.
    441     */
    442     const vector& operator-=(double d);
    443 
    444     /**
    445        \brief Multiply with scalar and assign operator, \f$ this_i =
    446        this_i * d \; \forall i \f$.
    447 
    448        \return A const reference to the resulting vector.
    449     */
    450     const vector& operator*=(double d);
    451 
     166    const vectorMutable& operator=(proxy);
    452167
    453168  private:
     169    void clean_up(void);
    454170
    455     /**
    456        \brief Create a new copy of the internal GSL vector.
    457 
    458        Necessary memory for the new GSL vector is allocated and the
    459        caller is responsible for freeing the allocated memory.
    460 
    461        \return A pointer to a copy of the internal GSL vector.
    462 
    463        \throw GSL_error if memory cannot be allocated for the new
    464        copy, or if dimensions mis-match.
    465     */
    466     gsl_vector* create_gsl_vector_copy(void) const;
    467 
    468     /**
    469        \brief Clear all dynamically allocated memory.
    470 
    471        Internal utility function.
    472     */
    473     void delete_allocated_memory(void);
    474 
    475     gsl_vector* v_;
    476171    gsl_vector_view* view_;
    477     const gsl_vector_const_view* view_const_;
    478     // proxy_v_ is used to access the proper underlying gsl_vector
    479     // in all const member functions. It is not used by design for
    480     // non-const vector functions and operators. This is to make sure
    481     // that runtime errors occur if a const vector is used in an
    482     // inappropriate manner such as on left hand side in assignment
    483     // (remember, v_ is null for const vector views).
    484     const gsl_vector* proxy_v_;
     172   
    485173  };
    486 
    487   /**
    488      \brief Check if all elements of the vector are zero.
    489 
    490      \return True if all elements in the vector is zero, false
    491      othwerwise.
    492   */
    493   bool isnull(const vector&);
    494 
    495   /**
    496      \brief Get the maximum value of the vector.
    497 
    498      \return The maximum value of the vector.
    499   */
    500   double max(const vector&);
    501 
    502   /**
    503      \brief Locate the maximum value in the vector.
    504 
    505      \return The index to the maximum value of the vector.
    506 
    507      \note Lower index has precedence.
    508   */
    509   size_t max_index(const vector&);
    510 
    511   /**
    512      \brief Get the minimum value of the vector.
    513 
    514      \return The minimum value of the vector.
    515   */
    516   double min(const vector&);
    517 
    518   /**
    519      \brief Locate the minimum value in the vector.
    520 
    521      \return The index to the minimum value of the vector.
    522 
    523      \note Lower index has precedence.
    524   */
    525   size_t min_index(const vector&);
    526 
    527   /**
    528      \brief Create a vector \a flag indicating NaN's in another vector
    529      \a templat.
    530 
    531      The \a flag vector is changed to contain 1's and 0's only. A 1
    532      means that the corresponding element in the \a templat vector is
    533      valid and a zero means that the corresponding element is a NaN.
    534 
    535      \note Space for vector \a flag is reallocated to fit the size of
    536      vector \a templat if sizes mismatch.
    537 
    538      \return True if the \a templat vector contains at least one NaN.
    539   */
    540   bool nan(const vector& templat, vector& flag);
    541 
    542   /**
    543      \brief Transforms a vector to a basis vector.
    544 
    545      All elements are set to zero except the \a i-th element which is
    546      set to one.
    547   */
    548   void set_basis(vector&, size_t i);
    549 
    550   /**
    551      Randomly shuffles the elements in vector \a invec
    552   */
    553   void shuffle(vector& invec);
    554 
    555   /**
    556      Sort the elements in the vector.
    557   */
    558   void sort(vector&);
    559 
    560   /**
    561      Create a vector \a sort_index containing the indeces of
    562      elements in a another vector \a invec.  The elements of \a
    563      sort_index give the index of the vector element which would
    564      have been stored in that position if the vector had been sorted
    565      in place. The first element of \a sort_index gives the index of the least
    566      element in \a invec, and the last element of \a sort_index gives the index of the
    567      greatest element in \a invec . The vector \a invec is not changed.
    568 
    569   */
    570   void sort_index(std::vector<size_t>& sort_index, const vector& invec);
    571 
    572   /** Similar to sort_index but creates a vector with indices to the \a k
    573   smallest elements in \a invec. 
    574   */
    575   void sort_smallest_index(std::vector<size_t>& sort_index, size_t k, const
    576   vector& invec);
    577 
    578   /** Similar to sort_index but creates a vector with indices to the \a k
    579   largest elements in \a invec. 
    580   */
    581   void sort_largest_index(std::vector<size_t>& sort_index, size_t k, const
    582   vector& invec);
    583 
    584  
    585 
    586   /**
    587      \brief Calculate the sum of all vector elements.
    588 
    589      \return The sum.
    590   */
    591   double sum(const vector&);
    592 
    593   /**
    594      \brief Swap vector elements by copying.
    595 
    596      The two vectors must have the same length.
    597 
    598      \throw GSL_error if vector lengths differs.
    599   */
    600   void swap(vector&, vector&);
    601 
    602   /**
    603      \brief The output operator for the vector class.
    604   */
    605   std::ostream& operator<<(std::ostream&, const vector&);
    606174
    607175}}} // of namespace utility, yat, and theplu
Note: See TracChangeset for help on using the changeset viewer.