source: trunk/yat/utility/vector.cc @ 714

Last change on this file since 714 was 714, checked in by Jari Häkkinen, 15 years ago

Addresses #170..

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.4 KB
Line 
1// $Id: vector.cc 714 2006-12-22 00:10:54Z jari $
2
3/*
4  Copyright (C) 2003 Daniel Dalevi, Peter Johansson
5  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
6  Copyright (C) 2005 Jari Häkkinen, Peter Johansson, Markus Ringnér
7  Copyright (C) 2006 Jari Häkkinen, Peter Johansson
8
9  This file is part of the yat library, http://lev.thep.lu.se/trac/yat
10
11  The yat library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU General Public License as
13  published by the Free Software Foundation; either version 2 of the
14  License, or (at your option) any later version.
15
16  The yat library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  General Public License for more details.
20
21  You should have received a copy of the GNU General Public License
22  along with this program; if not, write to the Free Software
23  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
24  02111-1307, USA.
25*/
26
27#include "vector.h"
28#include "matrix.h"
29#include "stl_utility.h"
30#include "utility.h"
31
32#include <iostream>
33#include <sstream>
34#include <vector>
35#include <utility>
36
37namespace theplu {
38namespace yat {
39namespace utility {
40
41
42  vector::vector(void)
43    : v_(NULL), v_const_(NULL), view_(NULL), view_const_(NULL), proxy_v_(NULL)
44  {
45  }
46
47
48  vector::vector(size_t n, double init_value)
49    : v_(gsl_vector_alloc(n)), v_const_(NULL), view_(NULL), view_const_(NULL),
50      proxy_v_(v_)
51  {
52    set_all(init_value);
53  }
54
55
56  vector::vector(const vector& other)
57    : v_(other.create_gsl_vector_copy()), v_const_(NULL), view_(NULL),
58      view_const_(NULL), proxy_v_(v_)
59  {
60  }
61
62
63  vector::vector(vector& v, size_t offset, size_t n, size_t stride)
64    : v_const_(NULL), view_const_(NULL)
65  {
66    // Jari, exception handling needed here. Failure in setting up a
67    // proper gsl_vector_view is communicated by NULL pointer in the
68    // view structure (cf. GSL manual). How about GSL error state?
69    view_ = new gsl_vector_view(gsl_vector_subvector_with_stride(v.v_,offset,
70                                                                 stride,n));
71    proxy_v_ = v_ = &(view_->vector);
72  }
73
74
75  vector::vector(const vector& v, size_t offset, size_t n, size_t stride)
76    : v_(NULL), view_(NULL)
77  {
78    // Jari, exception handling needed here. Failure in setting up a
79    // proper gsl_vector_view is communicated by NULL pointer in the
80    // view structure (cf. GSL manual). How about GSL error state?
81    view_const_ = new gsl_vector_const_view(
82                   gsl_vector_const_subvector_with_stride(v.v_,offset,stride,n));
83    proxy_v_ = v_const_ = &(view_const_->vector);
84  }
85
86
87  vector::vector(matrix& m, size_t i, bool row)
88    : v_const_(NULL), view_const_(NULL)
89  {
90    view_=new gsl_vector_view(row ?
91                              gsl_matrix_row   (m.gsl_matrix_p(),i) :
92                              gsl_matrix_column(m.gsl_matrix_p(),i)  );
93    proxy_v_ = v_ = &(view_->vector);
94  }
95
96
97  vector::vector(const matrix& m, size_t i, bool row)
98    : v_(NULL), view_(NULL)
99  {
100    view_const_ = new gsl_vector_const_view(row ?
101                                   gsl_matrix_const_row   (m.gsl_matrix_p(),i) :
102                                   gsl_matrix_const_column(m.gsl_matrix_p(),i) );
103    proxy_v_ = v_const_ = &(view_const_->vector);
104  }
105
106
107  vector::vector(std::istream& is, char sep)
108    throw (utility::IO_error, std::exception)
109    : v_const_(NULL), view_(NULL), view_const_(NULL)
110  {
111    // read the data file and store in stl vectors (dynamically
112    // expandable)
113    std::vector<std::vector<double> > data_matrix;
114    u_int nof_columns=0;
115    u_int nof_rows=0;
116    std::string line;
117    while(getline(is, line, '\n')) {
118      // Empty lines
119      if (!line.size())
120          continue;
121      nof_rows++;
122
123      std::vector<double> v;
124      std::string element;
125      std::stringstream ss(line);
126      bool ok=true;
127      while(ok) {
128        if(sep=='\0')
129          ok=(ss>>element);
130        else 
131          ok=getline(ss, element, sep);
132        if(!ok)
133          break;       
134
135        if(utility::is_double(element)) {
136          v.push_back(atof(element.c_str()));
137        }
138        else if (!element.size() || utility::is_nan(element)) {
139          v.push_back(std::numeric_limits<double>::quiet_NaN());
140        }
141        else {
142          // Jari, this should be communicated with as an exception.
143          // std::cerr << "Warning: '" << element
144          //           << "' is not an integer." << std::endl;
145        }
146      } 
147      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
148          v.push_back(std::numeric_limits<double>::quiet_NaN());
149      if (!nof_columns)
150        nof_columns=v.size();
151      else if (nof_rows && (nof_columns>1)) {
152        std::ostringstream s;
153        s << "vector::vector(std::istream&) data file error:\n"
154          << "    File has inconsistent number of rows (" << nof_rows
155          << ") and columns (" << nof_columns
156          << ").\n    Expected a row or a column vector.";
157        throw utility::IO_error(s.str());
158      }
159      else if (v.size()!=nof_columns) {
160        std::ostringstream s;
161        s << "vector::vector(std::istream&) data file error:\n"
162          << "    Line " << nof_rows << " has " << v.size()
163          << " columns; expected " << nof_columns << " column.";
164        throw utility::IO_error(s.str());
165      }
166      data_matrix.push_back(v);
167    }
168
169    // manipulate the state of the stream to be good
170    is.clear(std::ios::goodbit);
171    // convert the data to a gsl vector
172    proxy_v_ = v_ = gsl_vector_alloc(nof_rows*nof_columns);
173    size_t n=0;
174    for (size_t i=0; i<nof_rows; i++)
175      for (size_t j=0; j<nof_columns; j++) 
176        gsl_vector_set( v_, n++, data_matrix[i][j] );
177  }
178
179
180  vector::~vector(void)
181  {
182    if (view_)
183      delete view_;
184    else if (view_const_)
185      delete view_const_;
186    else if (v_)
187      gsl_vector_free(v_);
188  }
189
190
191  int vector::add(const vector& other)
192  {
193    return gsl_vector_add(v_,other.v_);
194  }
195
196
197  int vector::add(double term)
198  {
199    return gsl_vector_add_constant(v_,term);
200  }
201
202
203  gsl_vector* vector::create_gsl_vector_copy(void) const
204  {
205    gsl_vector* vec = gsl_vector_alloc(size());
206    gsl_vector_memcpy(vec, proxy_v_);
207    return vec;
208  }
209
210
211  int vector::div(const vector& other)
212  {
213    return gsl_vector_div(v_,other.v_);
214  }
215
216
217  const gsl_vector* vector::gsl_vector_p(void) const
218  {
219    return proxy_v_;
220  }
221
222
223  gsl_vector* vector::gsl_vector_p(void)
224  {
225    return v_;
226  }
227
228
229  bool vector::isnull(void) const
230  {
231    return gsl_vector_isnull(proxy_v_);
232  }
233
234
235  bool vector::isview(void) const
236  {
237    return view_ || view_const_;
238  }
239
240
241  double vector::max(void) const
242  {
243    return gsl_vector_max(proxy_v_);
244  }
245
246
247  size_t vector::max_index(void) const
248  {
249    return gsl_vector_max_index(proxy_v_);
250  }
251
252
253  double vector::min(void) const
254  {
255    return gsl_vector_min(proxy_v_);
256  }
257
258
259  size_t vector::min_index(void) const
260  {
261    return gsl_vector_min_index(proxy_v_);
262  }
263
264
265  std::pair<double,double> vector::minmax(void) const
266  {
267    double min, max;
268    gsl_vector_minmax(proxy_v_, &min, &max);
269    return std::pair<double,double>(min,max);
270  }
271
272
273  std::pair<size_t,size_t> vector::minmax_index(void) const
274  {
275    size_t min_index, max_index;
276    gsl_vector_minmax_index(proxy_v_, &min_index, &max_index);
277    return std::pair<size_t,size_t>(min_index,max_index);
278  }
279
280
281  int vector::mul(const vector& other)
282  {
283    return gsl_vector_mul(v_,other.v_);
284  }
285
286
287  int vector::reverse(void)
288  {
289    return gsl_vector_reverse(v_);
290  }
291
292
293  int vector::scale(double factor)
294  {
295    return gsl_vector_scale(v_,factor);
296  }
297
298
299  int vector::set(const vector& vec)
300  {
301    return gsl_vector_memcpy(v_,vec.v_);
302  }
303
304
305  void vector::set_all(const double& value)
306  {
307    gsl_vector_set_all(v_,value);
308  }
309
310
311  void vector::set_basis(const size_t i)
312  {
313    gsl_vector_set_basis(v_,i);
314  }
315
316
317  void vector::set_zero(void)
318  {
319    gsl_vector_set_zero(v_);
320  }
321
322
323  size_t vector::size(void) const
324  {
325    return proxy_v_->size;
326  }
327
328
329  void vector::sort(void)
330  {
331    gsl_sort_vector(v_);
332  }
333
334
335  int vector::sub(const vector& other)
336  {
337    return gsl_vector_sub(v_,other.v_);
338  }
339
340
341  double vector::sum(void) const
342  {
343    double sum = 0;
344    for (size_t i=0; i<size(); i++)
345      sum += (*this)(i);
346    return( sum );
347  }
348
349
350  int vector::swap(vector& other)
351  {
352    return gsl_vector_swap(v_,other.v_);
353  }
354
355
356  int vector::swap_elements(size_t i, size_t j)
357  {
358    return gsl_vector_swap_elements(v_,i,j);
359  }
360
361
362  double& vector::operator()(size_t i)
363  {
364    return *gsl_vector_ptr(v_,i);
365  }
366
367
368  const double& vector::operator()(size_t i) const
369  {
370    return *gsl_vector_const_ptr(proxy_v_,i);
371  }
372
373
374  double& vector::operator[](size_t i)
375  {
376    return *gsl_vector_ptr(v_,i);
377  }
378
379
380  const double& vector::operator[](size_t i) const
381  {
382    return *gsl_vector_const_ptr(proxy_v_,i);
383  }
384
385
386  bool vector::operator==(const vector& a) const
387  {
388    if (size()!=a.size())
389      return false;
390    for (size_t i=0; i<size(); ++i)
391      if (gsl_vector_get(proxy_v_,i)!=a(i))
392        return false;
393    return true;
394  }
395
396
397  double vector::operator*( const vector &other ) const
398  {
399    double res = 0.0;;
400    for ( size_t i = 0; i < size(); ++i ) 
401      res += other(i) * (*this)(i);
402    return res;
403  }
404
405
406  const vector& vector::operator=( const vector& other )
407  {
408    if( this != &other ) {
409      if (view_)
410        delete view_;
411      else if (view_const_)
412        delete view_const_;
413      else if ( v_ )
414        gsl_vector_free( v_ );
415      v_const_=NULL;
416      proxy_v_ = v_ = other.create_gsl_vector_copy();
417    }
418    return *this;
419  } 
420
421
422  const vector& vector::operator+=(const vector& other)
423  {
424    gsl_vector_add(v_,other.v_);
425    return *this;
426  }
427
428
429  const vector& vector::operator-=(const vector& other)
430  {
431    gsl_vector_sub(v_,other.v_);
432    return *this;
433  }
434
435
436  const vector& vector::operator*=(const double d)
437  {
438    gsl_vector_scale(v_,d);
439    return *this;
440  }
441
442
443  std::ostream& operator<<(std::ostream& s, const vector& a)
444  {
445    s.setf(std::ios::dec);
446    s.precision(12);
447    for (size_t j = 0; j < a.size(); ++j) {
448      s << a[j];
449      if ( (j+1)<a.size() )
450        s << s.fill();
451    }
452    return s;
453  }
454
455}}} // of namespace utility, yat, and thep
Note: See TracBrowser for help on using the repository browser.