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

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

Addresses #2. Continued adding GSL_error exceptions.

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