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

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

removing includes

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