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

Last change on this file since 808 was 808, checked in by Peter, 15 years ago

previous argument was invalid, but here is an implementation. Fixes #205

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.0 KB
Line 
1// $Id: vector.cc 808 2007-03-15 20:07:01Z peter $
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, 2007 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 "utility.h"
30
31#include <cassert>
32#include <cmath>
33#include <iostream>
34#include <sstream>
35#include <utility>
36#include <vector>
37
38namespace theplu {
39namespace yat {
40namespace utility {
41
42
43  vector::vector(void)
44    : v_(NULL), view_(NULL), view_const_(NULL), proxy_v_(NULL)
45  {
46  }
47
48
49  vector::vector(size_t n, double init_value)
50    : v_(gsl_vector_alloc(n)), view_(NULL), view_const_(NULL),
51      proxy_v_(v_)
52  {
53    if (!v_)
54      throw utility::GSL_error("vector::vector failed to allocate memory");
55
56    set_all(init_value);
57  }
58
59
60  vector::vector(const vector& other)
61    : v_(other.create_gsl_vector_copy()), view_(NULL),
62      view_const_(NULL), proxy_v_(v_)
63  {
64  }
65
66
67  vector::vector(vector& v, size_t offset, size_t n, size_t stride)
68    : view_const_(NULL)
69  {
70    view_ = new gsl_vector_view(gsl_vector_subvector_with_stride(v.v_,offset,
71                                                                 stride,n));
72    if (!view_)
73      throw utility::GSL_error("vector::vector failed to setup view");
74    proxy_v_ = v_ = &(view_->vector);
75  }
76
77
78  vector::vector(const vector& v, size_t offset, size_t n, size_t stride)
79    : v_(NULL), view_(NULL)
80  {
81    view_const_ = new gsl_vector_const_view(
82                   gsl_vector_const_subvector_with_stride(v.v_,offset,stride,n));
83    if (!view_const_)
84      throw utility::GSL_error("vector::vector failed to setup view");
85    proxy_v_ = &(view_const_->vector);
86  }
87
88
89  vector::vector(matrix& m, size_t i, bool row)
90    : view_const_(NULL)
91  {
92    view_=new gsl_vector_view(row ?
93                              gsl_matrix_row   (m.gsl_matrix_p(),i) :
94                              gsl_matrix_column(m.gsl_matrix_p(),i)  );
95    if (!view_)
96      throw utility::GSL_error("vector::vector failed to setup view");
97    proxy_v_ = v_ = &(view_->vector);
98  }
99
100
101  vector::vector(const matrix& m, size_t i, bool row)
102    : v_(NULL), view_(NULL)
103  {
104    view_const_ = new gsl_vector_const_view(row ?
105                                   gsl_matrix_const_row   (m.gsl_matrix_p(),i) :
106                                   gsl_matrix_const_column(m.gsl_matrix_p(),i) );
107    if (!view_const_)
108      throw utility::GSL_error("vector::vector failed to setup view");
109    proxy_v_ = &(view_const_->vector);
110  }
111
112
113  vector::vector(std::istream& is, char sep)
114    throw (utility::IO_error, std::exception)
115    : view_(NULL), view_const_(NULL)
116  {
117    // read the data file and store in stl vectors (dynamically
118    // expandable)
119    std::vector<std::vector<double> > data_matrix;
120    u_int nof_columns=0;
121    u_int nof_rows=0;
122    std::string line;
123    while(getline(is, line, '\n')) {
124      // Empty lines
125      if (!line.size())
126          continue;
127      nof_rows++;
128
129      std::vector<double> v;
130      std::string element;
131      std::stringstream ss(line);
132      bool ok=true;
133      while(ok) {
134        if(sep=='\0')
135          ok=(ss>>element);
136        else 
137          ok=getline(ss, element, sep);
138        if(!ok)
139          break;       
140
141        if(utility::is_double(element)) {
142          v.push_back(atof(element.c_str()));
143        }
144        else if (!element.size() || utility::is_nan(element)) {
145          v.push_back(std::numeric_limits<double>::quiet_NaN());
146        }
147        else {
148          std::stringstream ss("Warning: '");
149          ss << element << "' is not an integer.";
150          throw IO_error(ss.str());
151        }
152      } 
153      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
154          v.push_back(std::numeric_limits<double>::quiet_NaN());
155      if (!nof_columns)
156        nof_columns=v.size();
157      else if (nof_rows && (nof_columns>1)) {
158        std::ostringstream s;
159        s << "vector::vector(std::istream&) data file error:\n"
160          << "    File has inconsistent number of rows (" << nof_rows
161          << ") and columns (" << nof_columns
162          << ").\n    Expected a row or a column vector.";
163        throw utility::IO_error(s.str());
164      }
165      else if (v.size()!=nof_columns) {
166        std::ostringstream s;
167        s << "vector::vector(std::istream&) data file error:\n"
168          << "    Line " << nof_rows << " has " << v.size()
169          << " columns; expected " << nof_columns << " column.";
170        throw utility::IO_error(s.str());
171      }
172      data_matrix.push_back(v);
173    }
174
175    // manipulate the state of the stream to be good
176    is.clear(std::ios::goodbit);
177    // convert the data to a gsl vector
178    proxy_v_ = v_ = gsl_vector_alloc(nof_rows*nof_columns);
179    if (!v_)
180      throw utility::GSL_error("vector::vector failed to allocate memory");
181    size_t n=0;
182    // if gsl error handler disabled, out of bounds index will not
183    // abort the program.
184    for (size_t i=0; i<nof_rows; i++)
185      for (size_t j=0; j<nof_columns; j++) 
186        gsl_vector_set( v_, n++, data_matrix[i][j] );
187  }
188
189
190  vector::~vector(void)
191  {
192    if (view_)
193      delete view_;
194    else if (view_const_)
195      delete view_const_;
196    else if (v_)
197      gsl_vector_free(v_);
198  }
199
200
201  const vector& vector::clone(const vector& other)
202  {
203    if (this!=&other) {
204
205      if (view_)
206        delete view_;
207      else if (view_const_)
208        delete view_const_;
209      else if (v_)
210        gsl_vector_free(v_);
211
212      if (other.view_) {
213        view_ = new gsl_vector_view(*other.view_);
214        proxy_v_ = v_ = &(view_->vector);
215      }
216      else if (other.view_const_) {
217        view_const_ = new gsl_vector_const_view(*other.view_const_);
218        proxy_v_ = &(view_const_->vector);
219        v_=NULL;
220      }
221      else if (other.v_)
222        proxy_v_ = v_ = other.create_gsl_vector_copy();
223
224    }
225    return *this;
226  } 
227
228
229  gsl_vector* vector::create_gsl_vector_copy(void) const
230  {
231    gsl_vector* vec = gsl_vector_alloc(size());
232    if (!vec)
233      throw utility::GSL_error("vector::create_gsl_vector_copy failed to allocate memory");
234    if (gsl_vector_memcpy(vec, proxy_v_))
235      throw utility::GSL_error("vector::create_gsl_matrix_copy dimension mis-match");
236    return vec;
237  }
238
239
240  void vector::div(const vector& other)
241  {
242    assert(v_);
243    int status=gsl_vector_div(v_,other.gsl_vector_p());
244    if (status)
245      throw utility::GSL_error(std::string("vector::div",status));
246  }
247
248
249  bool vector::equal(const vector& other, const double d) const
250  {
251    if (this==&other)
252      return true;
253    if (size()!=other.size())
254      return false;
255    // if gsl error handler disabled, out of bounds index will not
256    // abort the program.
257    for (size_t i=0; i<size(); ++i)
258      // The two last condition checks are needed for NaN detection
259      if (fabs( (*this)(i)-other(i) ) > d ||
260          (*this)(i)!=(*this)(i) || other(i)!=other(i))
261        return false;
262    return true;
263  }
264
265
266  const gsl_vector* vector::gsl_vector_p(void) const
267  {
268    return proxy_v_;
269  }
270
271
272  gsl_vector* vector::gsl_vector_p(void)
273  {
274    return v_;
275  }
276
277
278  bool vector::isview(void) const
279  {
280    return view_ || view_const_;
281  }
282
283
284  void vector::mul(const vector& other)
285  {
286    assert(v_);
287    int status=gsl_vector_mul(v_,other.gsl_vector_p());
288    if (status)
289      throw utility::GSL_error(std::string("vector::mul",status));
290  }
291
292
293  void vector::resize(size_t n, double init_value)
294  {
295    if (view_)
296      delete view_;
297    else if (view_const_)
298      delete view_const_;
299    else if (v_)
300      gsl_vector_free(v_);
301    proxy_v_ = v_ = gsl_vector_alloc(n);
302
303    if (!v_)
304      throw utility::GSL_error("vector::vector failed to allocate memory");
305    set_all(init_value);
306  } 
307
308
309  void vector::reverse(void)
310  {
311    assert(v_);
312    gsl_vector_reverse(v_);
313  }
314
315
316  void vector::set(const vector& vec)
317  {
318    assert(v_);
319    if (gsl_vector_memcpy(v_,vec.gsl_vector_p()))
320      throw utility::GSL_error("vector::set dimension mis-match");
321  }
322
323
324  void vector::set_all(const double& value)
325  {
326    assert(v_);
327    gsl_vector_set_all(v_,value);
328  }
329
330
331  size_t vector::size(void) const
332  {
333    if (!proxy_v_)
334      return 0;
335    return proxy_v_->size;
336  }
337
338
339  void vector::swap(size_t i, size_t j)
340  {
341    assert(v_);
342    int status=gsl_vector_swap_elements(v_, i, j);
343    if (status)
344      throw utility::GSL_error(std::string("vector::swap_elements",status));
345  }
346
347
348  double& vector::operator()(size_t i)
349  {
350    assert(v_);
351    double* d=gsl_vector_ptr(v_, i);
352    if (!d)
353      throw utility::GSL_error("vector::operator()",GSL_EINVAL);
354    return *d;
355  }
356
357
358  const double& vector::operator()(size_t i) const
359  {
360    const double* d=gsl_vector_const_ptr(proxy_v_, i);
361    if (!d)
362      throw utility::GSL_error("vector::operator()",GSL_EINVAL);
363    return *d;
364  }
365
366
367  double& vector::operator[](size_t i)
368  {
369    return this->operator()(i);
370  }
371
372
373  const double& vector::operator[](size_t i) const
374  {
375    return this->operator()(i);
376  }
377
378
379  bool vector::operator==(const vector& other) const
380  {
381    return equal(other);
382  }
383
384
385  bool vector::operator!=(const vector& other) const
386  {
387    return !equal(other);
388  }
389
390
391  double vector::operator*( const vector &other ) const
392  {
393    double res = 0.0;;
394    for ( size_t i = 0; i < size(); ++i ) 
395      res += other(i) * (*this)(i);
396    return res;
397  }
398
399
400  const vector& vector::operator=( const vector& other )
401  {
402    set(other);
403    return *this;
404  } 
405
406
407  const vector& vector::operator+=(const vector& other)
408  {
409    assert(v_);
410    int status=gsl_vector_add(v_, other.gsl_vector_p());
411    if (status)
412      throw utility::GSL_error(std::string("vector::add", status));
413    return *this;
414  }
415
416
417  const vector& vector::operator+=(double d)
418  {
419    assert(v_);
420    gsl_vector_add_constant(v_, d);
421    return *this;
422  }
423
424
425  const vector& vector::operator-=(const vector& other)
426  {
427    assert(v_);
428    int status=gsl_vector_sub(v_, other.gsl_vector_p());
429    if (status)
430      throw utility::GSL_error(std::string("vector::sub", status));
431    return *this;
432  }
433
434
435  const vector& vector::operator-=(const double d)
436  {
437    assert(v_);
438    gsl_vector_add_constant(v_, -d);
439    return *this;
440  }
441
442
443  const vector& vector::operator*=(const double d)
444  {
445    assert(v_);
446    gsl_vector_scale(v_, d);
447    return *this;
448  }
449
450
451  bool isnull(const vector& v)
452  {
453    return gsl_vector_isnull(v.gsl_vector_p());
454  }
455
456
457  double max(const vector& v)
458  {
459    return gsl_vector_max(v.gsl_vector_p());
460  }
461
462
463  size_t max_index(const vector& v)
464  {
465    return gsl_vector_max_index(v.gsl_vector_p());
466  }
467
468
469  double min(const vector& v)
470  {
471    return gsl_vector_min(v.gsl_vector_p());
472  }
473
474
475  size_t min_index(const vector& v)
476  {
477    return gsl_vector_min_index(v.gsl_vector_p());
478  }
479
480
481  std::pair<double,double> minmax(const vector& v)
482  {
483    std::pair<double,double> minmax;
484    gsl_vector_minmax(v.gsl_vector_p(), &minmax.first, &minmax.second);
485    return minmax;
486  }
487
488
489  std::pair<size_t,size_t> minmax_index(const vector& v)
490  {
491    std::pair<size_t,size_t> minmax;
492    gsl_vector_minmax_index(v.gsl_vector_p(), &minmax.first, &minmax.second);
493    return minmax;
494  }
495
496
497  bool nan(const vector& templat, vector& flag)
498  {
499    size_t rows=templat.size();
500    if (rows!=flag.size())
501      flag.clone(vector(rows,1.0));
502    else
503      flag.set_all(1.0);
504    bool nan=false;
505    for (size_t i=0; i<rows; i++)
506      if (std::isnan(templat(i))) {
507        flag(i)=0;
508        nan=true;
509      }
510    return nan;
511  }
512
513
514  void set_basis(vector& v, size_t i)
515  {
516    assert(v.gsl_vector_p());
517    gsl_vector_set_basis(v.gsl_vector_p(),i);
518  }
519
520
521  void sort(vector& v)
522  {
523    assert(v.gsl_vector_p());
524    gsl_sort_vector(v.gsl_vector_p());
525  }
526
527
528  double sum(const vector& v)
529  {
530    double sum = 0;
531    size_t vsize=v.size();
532    for (size_t i=0; i<vsize; ++i)
533      sum += v[i];
534    return sum;
535  }
536
537
538  void swap(vector& v, vector& w)
539  {
540    assert(v.gsl_vector_p()); assert(w.gsl_vector_p());
541    int status=gsl_vector_swap(v.gsl_vector_p(),w.gsl_vector_p());
542    if (status)
543      throw utility::GSL_error(std::string("swap(vector&,vector&)",status));
544  }
545
546
547  std::ostream& operator<<(std::ostream& s, const vector& a)
548  {
549    s.setf(std::ios::dec);
550    s.precision(12);
551    for (size_t j = 0; j < a.size(); ++j) {
552      s << a[j];
553      if ( (j+1)<a.size() )
554        s << s.fill();
555    }
556    return s;
557  }
558
559}}} // of namespace utility, yat, and thep
Note: See TracBrowser for help on using the repository browser.