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

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

compiling 3 tests fail

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