source: trunk/lib/gslapi/vector.cc @ 560

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

fixed problem with ConsensusInputranker? - Now all tests do pass

  • 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 536 2006-03-03 21:19: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#include <c++_tools/classifier/DataLookup1D.h>
8
9#include <iostream>
10#include <sstream>
11#include <vector>
12#include <utility>
13
14
15namespace theplu {
16namespace gslapi {
17
18
19
20  vector::vector(vector& v, size_t offset, size_t n, size_t stride)
21    : const_view_(NULL)
22  {
23    // Jari, exception handling needed here. Failure in setting up a
24    // proper gsl_vector_view is communicated by NULL pointer in the
25    // view structure (cf. GSL manual). How about GSL error state?
26    view_ = new gsl_vector_view(gsl_vector_subvector_with_stride(v.v_,offset,
27                                                                 stride,n));
28    v_ = &(view_->vector);
29  }
30
31
32
33  vector::vector(const vector& v, size_t offset, size_t n, size_t stride)
34    : view_(NULL)
35  {
36    // Jari, exception handling needed here. Failure in setting up a
37    // proper gsl_vector_view is communicated by NULL pointer in the
38    // view structure (cf. GSL manual). How about GSL error state?
39    const_view_ = new gsl_vector_const_view(
40                   gsl_vector_const_subvector_with_stride(v.v_,offset,stride,n));
41    v_ = const_cast<gsl_vector*>(&(const_view_->vector));
42  }
43
44
45
46  vector::vector(matrix& m, size_t i, bool row)
47    : const_view_(NULL)
48  {
49    view_=new gsl_vector_view(row ?
50                              gsl_matrix_row   (m.gsl_matrix_p(),i) :
51                              gsl_matrix_column(m.gsl_matrix_p(),i)  );
52    v_ = &(view_->vector);
53  }
54
55
56
57  vector::vector(const matrix& m, size_t i, bool row)
58    : view_(NULL)
59  {
60    const_view_ = new gsl_vector_const_view(row ?
61                                   gsl_matrix_const_row   (m.gsl_matrix_p(),i) :
62                                   gsl_matrix_const_column(m.gsl_matrix_p(),i) );
63    v_ = const_cast<gsl_vector*>(&(const_view_->vector));
64  }
65
66
67
68  vector::vector(std::istream& is, char sep) 
69    throw (utility::IO_error,std::exception)
70    : view_(NULL), const_view_(NULL)
71  {
72    // read the data file and store in stl vectors (dynamically
73    // expandable)
74    std::vector<std::vector<double> > data_matrix;
75    u_int nof_columns=0;
76    u_int nof_rows=0;
77    std::string line;
78    while(getline(is, line, '\n')) {
79      // Empty lines
80      if (!line.size())
81          continue;
82      nof_rows++;
83     
84      std::vector<double> v;
85      std::string element;
86      std::stringstream ss(line);
87      bool ok=true;
88      while(ok) {
89        if(sep=='\0')
90          ok=(ss>>element);
91        else 
92          ok=getline(ss, element, sep);
93        if(!ok)
94          break;       
95
96        if(utility::is_double(element)) {
97          v.push_back(atof(element.c_str()));
98        }
99        else if (!element.size() || utility::is_nan(element)) {
100          v.push_back(std::numeric_limits<double>::quiet_NaN());
101        }
102        else {
103          // Jari, this should be communicated with as an exception.
104          // std::cerr << "Warning: '" << element
105          //           << "' is not an integer." << std::endl;
106        }
107      } 
108      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
109          v.push_back(std::numeric_limits<double>::quiet_NaN());
110      if (!nof_columns)
111        nof_columns=v.size();
112      else if (nof_rows && (nof_columns>1)) {
113        std::ostringstream s;
114        s << "vector::vector(std::istream&) data file error:\n"
115          << "    File has inconsistent number of rows (" << nof_rows
116          << ") and columns (" << nof_columns
117          << ").\n    Expected a row or a column vector.";
118        throw utility::IO_error(s.str());
119      }
120      else if (v.size()!=nof_columns) {
121        std::ostringstream s;
122        s << "vector::vector(std::istream&) data file error:\n"
123          << "    Line " << nof_rows << " has " << v.size()
124          << " columns; expected " << nof_columns << " column.";
125        throw utility::IO_error(s.str());
126      }
127      data_matrix.push_back(v);
128    }
129
130    // manipulate the state of the stream to be good
131    is.clear(std::ios::goodbit);
132    // convert the data to a gsl vector
133    v_ = gsl_vector_alloc(nof_rows*nof_columns);
134    size_t n=0;
135    for (size_t i=0; i<nof_rows; i++)
136      for (size_t j=0; j<nof_columns; j++) 
137        gsl_vector_set( v_, n++, data_matrix[i][j] );
138  }
139
140
141
142  vector::vector(const classifier::DataLookup1D& data)
143    : view_(NULL), const_view_(NULL)
144  {
145    v_ = gsl_vector_alloc(data.size());
146    for (size_t i=0; i<data.size(); i++)
147      gsl_vector_set( v_, i, data(i) );
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.