source: branches/peters_vector/lib/gslapi/vector.cc @ 469

Last change on this file since 469 was 469, checked in by Peter, 16 years ago

non compiling checking before revision after design meeting

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.6 KB
Line 
1// $Id: vector.cc 469 2005-12-19 14:58:29Z peter $
2
3#include <c++_tools/gslapi/vector.h>
4#include <c++_tools/gslapi/matrix.h>
5#include <c++_tools/utility/stl_utility.h>
6#include <c++_tools/utility/utility.h>
7
8#include <sstream>
9#include <vector>
10#include <utility>
11
12
13namespace theplu {
14namespace gslapi {
15
16
17
18  vector::vector(vector& v, size_t offset, size_t n, size_t stride)
19    : const_view_(NULL)
20  {
21    // Jari, exception handling needed here. Failure in setting up a
22    // proper gsl_vector_view is communicated by NULL pointer in the
23    // view structure (cf. GSL manual). How about GSL error state?
24    view_ = new gsl_vector_view(gsl_vector_subvector_with_stride(v.v_,offset,
25                                                                 stride,n));
26    v_ = &(view_->vector);
27  }
28
29
30
31  vector::vector(const vector& v, size_t offset, size_t n, size_t stride)
32    : view_(NULL)
33  {
34    // Jari, exception handling needed here. Failure in setting up a
35    // proper gsl_vector_view is communicated by NULL pointer in the
36    // view structure (cf. GSL manual). How about GSL error state?
37    const_view_ = new gsl_vector_const_view(
38                   gsl_vector_const_subvector_with_stride(v.v_,offset,stride,n));
39    v_ = const_cast<gsl_vector*>(&(const_view_->vector));
40  }
41
42
43
44  vector::vector(matrix& m, size_t i, bool row)
45    : const_view_(NULL)
46  {
47    view_=new gsl_vector_view(row ?
48                              gsl_matrix_row   (m.gsl_matrix_p(),i) :
49                              gsl_matrix_column(m.gsl_matrix_p(),i)  );
50    v_ = &(view_->vector);
51  }
52
53
54
55  vector::vector(const matrix& m, size_t i, bool row)
56    : view_(NULL)
57  {
58    const_view_ = new gsl_vector_const_view(row ?
59                                   gsl_matrix_const_row   (m.gsl_matrix_p(),i) :
60                                   gsl_matrix_const_column(m.gsl_matrix_p(),i) );
61    v_ = const_cast<gsl_vector*>(&(const_view_->vector));
62  }
63
64
65
66  vector::vector(std::istream& is, char sep) 
67    throw (utility::IO_error,std::exception)
68    : view_(NULL), const_view_(NULL)
69  {
70    // read the data file and store in stl vectors (dynamically
71    // expandable)
72    std::vector<std::vector<double> > data_matrix;
73    u_int nof_columns=0;
74    u_int nof_rows=0;
75    std::string line;
76    while(getline(is, line, '\n')) {
77      // Empty lines
78      if (!line.size())
79          continue;
80      nof_rows++;
81     
82      std::vector<double> v;
83      std::string element;
84      std::stringstream ss(line);
85      bool ok=true;
86      while(ok) {
87        if(sep=='\0')
88          ok=(ss>>element);
89        else 
90          ok=getline(ss, element, sep);
91        if(!ok)
92          break;       
93
94        if(utility::is_double(element)) {
95          v.push_back(atof(element.c_str()));
96        }
97        else if (!element.size() || utility::is_nan(element)) {
98          v.push_back(std::numeric_limits<double>::quiet_NaN());
99        }
100        else {
101          // Jari, this should be communicated with as an exception.
102          // std::cerr << "Warning: '" << element
103          //           << "' is not an integer." << std::endl;
104        }
105      } 
106      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
107          v.push_back(std::numeric_limits<double>::quiet_NaN());
108      if (!nof_columns)
109        nof_columns=v.size();
110      else if (nof_rows && (nof_columns>1)) {
111        std::ostringstream s;
112        s << "vector::vector(std::istream&) data file error:\n"
113          << "    File has inconsistent number of rows (" << nof_rows
114          << ") and columns (" << nof_columns
115          << ").\n    Expected a row or a column vector.";
116        throw utility::IO_error(s.str());
117      }
118      else if (v.size()!=nof_columns) {
119        std::ostringstream s;
120        s << "vector::vector(std::istream&) data file error:\n"
121          << "    Line " << nof_rows << " has " << v.size()
122          << " columns; expected " << nof_columns << " column.";
123        throw utility::IO_error(s.str());
124      }
125      data_matrix.push_back(v);
126    }
127
128    // manipulate the state of the stream to be good
129    is.clear(std::ios::goodbit);
130    // convert the data to a gsl vector
131    v_ = gsl_vector_alloc(nof_rows*nof_columns);
132    size_t n=0;
133    for (size_t i=0; i<nof_rows; i++)
134      for (size_t j=0; j<nof_columns; j++) 
135        gsl_vector_set( v_, n++, data_matrix[i][j] );
136  }
137
138
139
140  vector::~vector(void)
141  {
142    if (view_)
143      delete view_;
144    else if (const_view_)
145      delete const_view_;
146    else if (v_)
147      gsl_vector_free(v_);
148    v_=NULL;
149  }
150
151
152
153  gsl_vector* vector::create_gsl_vector_copy(void) const
154  {
155    gsl_vector* vec = gsl_vector_alloc(size());
156    gsl_vector_memcpy(vec,v_);
157    return vec;
158  }
159
160
161
162  std::pair<double,double> vector::minmax(void) const
163  {
164    double min, max;
165    gsl_vector_minmax(v_,&min,&max);
166    return std::pair<double,double>(min,max);
167  }
168
169
170
171  std::pair<size_t,size_t> vector::minmax_index(void) const
172  {
173    size_t min_index, max_index;
174    gsl_vector_minmax_index(v_,&min_index,&max_index);
175    return std::pair<size_t,size_t>(min_index,max_index);
176  }
177
178
179
180  double VectorAbstract::sum(void) const
181  {
182    double sum = 0;
183    for (size_t i=0; i<size(); i++)
184      sum += (*this)(i);
185    return( sum );
186  } 
187
188  bool vector::operator==(const vector& a) const
189  {
190    if (size()!=a.size())
191      return false;
192    for (size_t i=0; i<size(); ++i)
193      if (gsl_vector_get(v_,i)!=a(i))
194        return false;
195    return true;
196  }
197
198
199
200  double VectorAbstract::operator*( const VectorAbstract &other ) const
201  {
202    double res = 0.0;;
203    for ( size_t i = 0; i < size(); ++i ) 
204      res += other(i) * (*this)(i);
205    return res;
206  }
207
208
209
210  const vector& vector::operator=( const vector& other )
211  {
212    if( this != &other ) {
213      if (view_)
214        delete view_;
215      else if (const_view_)
216        delete const_view_;
217      else if ( v_ )
218        gsl_vector_free( v_ );
219      v_ = other.create_gsl_vector_copy();
220    }
221    return *this;
222  } 
223
224
225
226  std::ostream& operator<<(std::ostream& s, const vector& a)
227  {
228    s.setf(std::ios::dec);
229    s.precision(12);
230    for (size_t j = 0; j < a.size(); ++j) {
231      s << a[j];
232      if ( (j+1)<a.size() )
233        s << " ";
234    }
235
236    return s;
237  }
238
239
240}} // of namespace gslapi and namespace thep
Note: See TracBrowser for help on using the repository browser.