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

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

Fixes #171 and addresses #2. A few more GSL_error exceptions. Removed Jari comments.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.5 KB
Line 
1// $Id: vector.cc 759 2007-02-19 19:41:25Z 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    view_ = new gsl_vector_view(gsl_vector_subvector_with_stride(v.v_,offset,
70                                                                 stride,n));
71    if (!view_)
72      throw utility::GSL_error("vector::vector failed to setup view");
73    proxy_v_ = v_ = &(view_->vector);
74  }
75
76
77  vector::vector(const vector& v, size_t offset, size_t n, size_t stride)
78    : v_(NULL), view_(NULL)
79  {
80    view_const_ = new gsl_vector_const_view(
81                   gsl_vector_const_subvector_with_stride(v.v_,offset,stride,n));
82    if (!view_const_)
83      throw utility::GSL_error("vector::vector failed to setup view");
84    proxy_v_ = v_const_ = &(view_const_->vector);
85  }
86
87
88  vector::vector(matrix& m, size_t i, bool row)
89    : v_const_(NULL), view_const_(NULL)
90  {
91    view_=new gsl_vector_view(row ?
92                              gsl_matrix_row   (m.gsl_matrix_p(),i) :
93                              gsl_matrix_column(m.gsl_matrix_p(),i)  );
94    if (!view_)
95      throw utility::GSL_error("vector::vector failed to setup view");
96    proxy_v_ = v_ = &(view_->vector);
97  }
98
99
100  vector::vector(const matrix& m, size_t i, bool row)
101    : v_(NULL), view_(NULL)
102  {
103    view_const_ = new gsl_vector_const_view(row ?
104                                   gsl_matrix_const_row   (m.gsl_matrix_p(),i) :
105                                   gsl_matrix_const_column(m.gsl_matrix_p(),i) );
106    if (!view_const_)
107      throw utility::GSL_error("vector::vector failed to setup view");
108    proxy_v_ = v_const_ = &(view_const_->vector);
109  }
110
111
112  vector::vector(std::istream& is, char sep)
113    throw (utility::IO_error, std::exception)
114    : v_const_(NULL), view_(NULL), view_const_(NULL)
115  {
116    // read the data file and store in stl vectors (dynamically
117    // expandable)
118    std::vector<std::vector<double> > data_matrix;
119    u_int nof_columns=0;
120    u_int nof_rows=0;
121    std::string line;
122    while(getline(is, line, '\n')) {
123      // Empty lines
124      if (!line.size())
125          continue;
126      nof_rows++;
127
128      std::vector<double> v;
129      std::string element;
130      std::stringstream ss(line);
131      bool ok=true;
132      while(ok) {
133        if(sep=='\0')
134          ok=(ss>>element);
135        else 
136          ok=getline(ss, element, sep);
137        if(!ok)
138          break;       
139
140        if(utility::is_double(element)) {
141          v.push_back(atof(element.c_str()));
142        }
143        else if (!element.size() || utility::is_nan(element)) {
144          v.push_back(std::numeric_limits<double>::quiet_NaN());
145        }
146        else {
147          std::stringstream ss("Warning: '");
148          ss << element << "' is not an integer.";
149          throw IO_error(ss.str());
150        }
151      } 
152      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
153          v.push_back(std::numeric_limits<double>::quiet_NaN());
154      if (!nof_columns)
155        nof_columns=v.size();
156      else if (nof_rows && (nof_columns>1)) {
157        std::ostringstream s;
158        s << "vector::vector(std::istream&) data file error:\n"
159          << "    File has inconsistent number of rows (" << nof_rows
160          << ") and columns (" << nof_columns
161          << ").\n    Expected a row or a column vector.";
162        throw utility::IO_error(s.str());
163      }
164      else if (v.size()!=nof_columns) {
165        std::ostringstream s;
166        s << "vector::vector(std::istream&) data file error:\n"
167          << "    Line " << nof_rows << " has " << v.size()
168          << " columns; expected " << nof_columns << " column.";
169        throw utility::IO_error(s.str());
170      }
171      data_matrix.push_back(v);
172    }
173
174    // manipulate the state of the stream to be good
175    is.clear(std::ios::goodbit);
176    // convert the data to a gsl vector
177    proxy_v_ = v_ = gsl_vector_alloc(nof_rows*nof_columns);
178    if (!v_)
179      throw utility::GSL_error("vector::vector failed to allocate memory");
180    size_t n=0;
181    // if gsl error handler disabled, out of bounds index will not
182    // abort the program.
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  void vector::add(const vector& other)
201  {
202    int status=gsl_vector_add(v_,other.v_);
203    if (status)
204      throw utility::GSL_error(std::string("vector::add",status));
205  }
206
207
208  void vector::add(double term)
209  {
210    gsl_vector_add_constant(v_,term);
211  }
212
213
214  gsl_vector* vector::create_gsl_vector_copy(void) const
215  {
216    gsl_vector* vec = gsl_vector_alloc(size());
217    if (!vec)
218      throw utility::GSL_error("vector::create_gsl_vector_copy failed to allocate memory");
219    if (gsl_vector_memcpy(vec, proxy_v_))
220      throw utility::GSL_error("vector::create_gsl_matrix_copy dimension mis-match");
221    return vec;
222  }
223
224
225  void vector::div(const vector& other)
226  {
227    int status=gsl_vector_div(v_,other.v_);
228    if (status)
229      throw utility::GSL_error(std::string("vector::div",status));
230  }
231
232
233  const gsl_vector* vector::gsl_vector_p(void) const
234  {
235    return proxy_v_;
236  }
237
238
239  gsl_vector* vector::gsl_vector_p(void)
240  {
241    return v_;
242  }
243
244
245  bool vector::isnull(void) const
246  {
247    return gsl_vector_isnull(proxy_v_);
248  }
249
250
251  bool vector::isview(void) const
252  {
253    return view_ || view_const_;
254  }
255
256
257  double vector::max(void) const
258  {
259    return gsl_vector_max(proxy_v_);
260  }
261
262
263  size_t vector::max_index(void) const
264  {
265    return gsl_vector_max_index(proxy_v_);
266  }
267
268
269  double vector::min(void) const
270  {
271    return gsl_vector_min(proxy_v_);
272  }
273
274
275  size_t vector::min_index(void) const
276  {
277    return gsl_vector_min_index(proxy_v_);
278  }
279
280
281  std::pair<double,double> vector::minmax(void) const
282  {
283    double min, max;
284    gsl_vector_minmax(proxy_v_, &min, &max);
285    return std::pair<double,double>(min,max);
286  }
287
288
289  std::pair<size_t,size_t> vector::minmax_index(void) const
290  {
291    size_t min_index, max_index;
292    gsl_vector_minmax_index(proxy_v_, &min_index, &max_index);
293    return std::pair<size_t,size_t>(min_index,max_index);
294  }
295
296
297  void vector::mul(const vector& other)
298  {
299    int status=gsl_vector_mul(v_,other.v_);
300    if (status)
301      throw utility::GSL_error(std::string("vector::mul",status));
302  }
303
304
305  void vector::reverse(void)
306  {
307    gsl_vector_reverse(v_);
308  }
309
310
311  void vector::scale(double factor)
312  {
313    gsl_vector_scale(v_,factor);
314  }
315
316
317  void vector::set(const vector& vec)
318  {
319    if (gsl_vector_memcpy(v_,vec.v_))
320      throw utility::GSL_error("vector::set dimension mis-match");
321  }
322
323
324  void vector::set_all(const double& value)
325  {
326    gsl_vector_set_all(v_,value);
327  }
328
329
330  void vector::set_basis(const size_t i)
331  {
332    gsl_vector_set_basis(v_,i);
333  }
334
335
336  void vector::set_zero(void)
337  {
338    gsl_vector_set_zero(v_);
339  }
340
341
342  size_t vector::size(void) const
343  {
344    if (!proxy_v_)
345      return 0;
346    return proxy_v_->size;
347  }
348
349
350  void vector::sort(void)
351  {
352    gsl_sort_vector(v_);
353  }
354
355
356  void vector::sub(const vector& other)
357  {
358    int status=gsl_vector_sub(v_,other.v_);
359    if (status)
360      throw utility::GSL_error(std::string("vector::sub",status));
361  }
362
363
364  double vector::sum(void) const
365  {
366    double sum = 0;
367    for (size_t i=0; i<size(); i++)
368      sum += (*this)(i);
369    return( sum );
370  }
371
372
373  void vector::swap(vector& other)
374  {
375    int status=gsl_vector_swap(v_,other.v_);
376    if (status)
377      throw utility::GSL_error(std::string("vector::swap",status));
378  }
379
380
381  void vector::swap_elements(size_t i, size_t j)
382  {
383    int status=gsl_vector_swap_elements(v_, i, j);
384    if (status)
385      throw utility::GSL_error(std::string("vector::swap_elements",status));
386  }
387
388
389  double& vector::operator()(size_t i)
390  {
391    double* d=gsl_vector_ptr(v_, i);
392    if (!d)
393      throw utility::GSL_error("vector::operator()",GSL_EINVAL);
394    return *d;
395  }
396
397
398  const double& vector::operator()(size_t i) const
399  {
400    const double* d=gsl_vector_const_ptr(proxy_v_, i);
401    if (!d)
402      throw utility::GSL_error("vector::operator()",GSL_EINVAL);
403    return *d;
404  }
405
406
407  double& vector::operator[](size_t i)
408  {
409    return this->operator()(i);
410  }
411
412
413  const double& vector::operator[](size_t i) const
414  {
415    return this->operator()(i);
416  }
417
418
419  bool vector::operator==(const vector& a) const
420  {
421    if (size()!=a.size())
422      return false;
423    // if gsl error handler disabled, out of bounds index will not
424    // abort the program.
425    for (size_t i=0; i<size(); ++i)
426      if (gsl_vector_get(proxy_v_,i)!=a(i))
427        return false;
428    return true;
429  }
430
431
432  double vector::operator*( const vector &other ) const
433  {
434    double res = 0.0;;
435    for ( size_t i = 0; i < size(); ++i ) 
436      res += other(i) * (*this)(i);
437    return res;
438  }
439
440
441  const vector& vector::operator=( const vector& other )
442  {
443    if( this != &other ) {
444      if (view_)
445        delete view_;
446      else if (view_const_)
447        delete view_const_;
448      else if ( v_ )
449        gsl_vector_free( v_ );
450      v_const_=NULL;
451      proxy_v_ = v_ = other.create_gsl_vector_copy();
452    }
453    return *this;
454  } 
455
456
457  const vector& vector::operator+=(const vector& other)
458  {
459    add(other);
460    return *this;
461  }
462
463
464  const vector& vector::operator-=(const vector& other)
465  {
466    sub(other);
467    return *this;
468  }
469
470
471  const vector& vector::operator*=(const double d)
472  {
473    scale(d);
474    return *this;
475  }
476
477
478  std::ostream& operator<<(std::ostream& s, const vector& a)
479  {
480    s.setf(std::ios::dec);
481    s.precision(12);
482    for (size_t j = 0; j < a.size(); ++j) {
483      s << a[j];
484      if ( (j+1)<a.size() )
485        s << s.fill();
486    }
487    return s;
488  }
489
490}}} // of namespace utility, yat, and thep
Note: See TracBrowser for help on using the repository browser.