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

Last change on this file since 880 was 880, checked in by Peter, 14 years ago

refs #244

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.1 KB
Line 
1// $Id: vector.cc 880 2007-09-21 23:43:22Z peter $
2
3/*
4  Copyright (C) 2003 Daniel Dalevi, Peter Johansson
5  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
6  Copyright (C) 2005, 2006 Jari Häkkinen, Markus Ringnér, Peter Johansson
7  Copyright (C) 2007 Jari Häkkinen, Peter Johansson
8
9  This file is part of the yat library, http://trac.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#include "yat/random/random.h"
31
32#include <cassert>
33#include <cmath>
34#include <iostream>
35#include <sstream>
36#include <utility>
37#include <vector>
38
39namespace theplu {
40namespace yat {
41namespace utility {
42
43
44  vector::vector(void)
45    : v_(NULL), view_(NULL), view_const_(NULL), proxy_v_(NULL)
46  {
47  }
48
49
50  vector::vector(size_t n, double init_value)
51    : v_(gsl_vector_alloc(n)), view_(NULL), view_const_(NULL),
52      proxy_v_(v_)
53  {
54    if (!v_)
55      throw utility::GSL_error("vector::vector failed to allocate memory");
56
57    all(init_value);
58  }
59
60
61  vector::vector(const vector& other)
62    : v_(other.create_gsl_vector_copy()), view_(NULL),
63      view_const_(NULL), proxy_v_(v_)
64  {
65  }
66
67
68  vector::vector(vector& v, size_t offset, size_t n, size_t stride)
69    : view_const_(NULL)
70  {
71    view_ = new gsl_vector_view(gsl_vector_subvector_with_stride(v.v_,offset,
72                                                                 stride,n));
73    if (!view_)
74      throw utility::GSL_error("vector::vector failed to setup view");
75    proxy_v_ = v_ = &(view_->vector);
76  }
77
78
79  vector::vector(const vector& v, size_t offset, size_t n, size_t stride)
80    : v_(NULL), view_(NULL)
81  {
82    view_const_ = new gsl_vector_const_view(
83                   gsl_vector_const_subvector_with_stride(v.v_,offset,stride,n));
84    if (!view_const_)
85      throw utility::GSL_error("vector::vector failed to setup view");
86    proxy_v_ = &(view_const_->vector);
87  }
88
89
90  vector::vector(matrix& m, size_t i, bool row)
91    : view_const_(NULL)
92  {
93    view_=new gsl_vector_view(row ?
94                              gsl_matrix_row   (m.gsl_matrix_p(),i) :
95                              gsl_matrix_column(m.gsl_matrix_p(),i)  );
96    if (!view_)
97      throw utility::GSL_error("vector::vector failed to setup view");
98    proxy_v_ = v_ = &(view_->vector);
99  }
100
101
102  vector::vector(const matrix& m, size_t i, bool row)
103    : v_(NULL), view_(NULL)
104  {
105    view_const_ = new gsl_vector_const_view(row ?
106                                   gsl_matrix_const_row   (m.gsl_matrix_p(),i) :
107                                   gsl_matrix_const_column(m.gsl_matrix_p(),i) );
108    if (!view_const_)
109      throw utility::GSL_error("vector::vector failed to setup view");
110    proxy_v_ = &(view_const_->vector);
111  }
112
113
114  vector::vector(std::istream& is, char sep)
115    throw (utility::IO_error, std::exception)
116    : 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          std::stringstream ss("Warning: '");
150          ss << element << "' is not a double.";
151          throw IO_error(ss.str());
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    // if gsl error handler disabled, out of bounds index will not
184    // abort the program.
185    for (size_t i=0; i<nof_rows; i++)
186      for (size_t j=0; j<nof_columns; j++) 
187        gsl_vector_set( v_, n++, data_matrix[i][j] );
188  }
189
190
191  vector::~vector(void)
192  {
193    delete_allocated_memory();
194  }
195
196
197  vector::iterator vector::begin(void)
198  {
199    return vector::iterator(*this, 0);
200  }
201
202
203  vector::const_iterator vector::begin(void) const
204  {
205    return vector::const_iterator(*this, 0);
206  }
207
208
209  const vector& vector::clone(const vector& other)
210  {
211    if (this!=&other) {
212
213      delete_allocated_memory();
214
215      if (other.view_) {
216        view_ = new gsl_vector_view(*other.view_);
217        proxy_v_ = v_ = &(view_->vector);
218      }
219      else if (other.view_const_) {
220        view_const_ = new gsl_vector_const_view(*other.view_const_);
221        proxy_v_ = &(view_const_->vector);
222        v_=NULL;
223      }
224      else if (other.v_)
225        proxy_v_ = v_ = other.create_gsl_vector_copy();
226
227    }
228    return *this;
229  } 
230
231
232  gsl_vector* vector::create_gsl_vector_copy(void) const
233  {
234    gsl_vector* vec = gsl_vector_alloc(size());
235    if (!vec)
236      throw utility::GSL_error("vector::create_gsl_vector_copy failed to allocate memory");
237    if (gsl_vector_memcpy(vec, proxy_v_))
238      throw utility::GSL_error("vector::create_gsl_matrix_copy dimension mis-match");
239    return vec;
240  }
241
242
243  void vector::delete_allocated_memory(void)
244  {
245    if (view_)
246      delete view_;
247    else if (view_const_)
248      delete view_const_;
249    else if (v_)
250      gsl_vector_free(v_);
251    proxy_v_=v_=NULL;
252  }
253
254
255  void vector::div(const vector& other)
256  {
257    assert(v_);
258    int status=gsl_vector_div(v_,other.gsl_vector_p());
259    if (status)
260      throw utility::GSL_error(std::string("vector::div",status));
261  }
262
263
264  vector::iterator vector::end(void)
265  {
266    return vector::iterator(*this, size());
267  }
268
269
270  vector::const_iterator vector::end(void) const
271  {
272    return vector::const_iterator(*this, size());
273  }
274
275
276  bool vector::equal(const vector& other, const double d) const
277  {
278    if (this==&other)
279      return true;
280    if (size()!=other.size())
281      return false;
282    // if gsl error handler disabled, out of bounds index will not
283    // abort the program.
284    for (size_t i=0; i<size(); ++i)
285      // The two last condition checks are needed for NaN detection
286      if (fabs( (*this)(i)-other(i) ) > d ||
287          (*this)(i)!=(*this)(i) || other(i)!=other(i))
288        return false;
289    return true;
290  }
291
292
293  const gsl_vector* vector::gsl_vector_p(void) const
294  {
295    return proxy_v_;
296  }
297
298
299  gsl_vector* vector::gsl_vector_p(void)
300  {
301    return v_;
302  }
303
304
305  bool vector::isview(void) const
306  {
307    return view_ || view_const_;
308  }
309
310
311  void vector::mul(const vector& other)
312  {
313    assert(v_);
314    int status=gsl_vector_mul(v_,other.gsl_vector_p());
315    if (status)
316      throw utility::GSL_error(std::string("vector::mul",status));
317  }
318
319
320  void vector::resize(size_t n, double init_value)
321  {
322    delete_allocated_memory();
323    proxy_v_ = v_ = gsl_vector_alloc(n);
324    if (!v_)
325      throw utility::GSL_error("vector::vector failed to allocate memory");
326    all(init_value);
327  } 
328
329
330  void vector::reverse(void)
331  {
332    assert(v_);
333    gsl_vector_reverse(v_);
334  }
335
336
337  void vector::all(const double& value)
338  {
339    assert(v_);
340    gsl_vector_set_all(v_,value);
341  }
342
343
344  size_t vector::size(void) const
345  {
346    if (!proxy_v_)
347      return 0;
348    return proxy_v_->size;
349  }
350
351
352  void vector::swap(size_t i, size_t j)
353  {
354    assert(v_);
355    int status=gsl_vector_swap_elements(v_, i, j);
356    if (status)
357      throw utility::GSL_error(std::string("vector::swap_elements",status));
358  }
359
360
361  double& vector::operator()(size_t i)
362  {
363    assert(v_);
364    double* d=gsl_vector_ptr(v_, i);
365    if (!d)
366      throw utility::GSL_error("vector::operator()",GSL_EINVAL);
367    return *d;
368  }
369
370
371  const double& vector::operator()(size_t i) const
372  {
373    const double* d=gsl_vector_const_ptr(proxy_v_, i);
374    if (!d)
375      throw utility::GSL_error("vector::operator()",GSL_EINVAL);
376    return *d;
377  }
378
379
380  double& vector::operator[](size_t i)
381  {
382    return this->operator()(i);
383  }
384
385
386  const double& vector::operator[](size_t i) const
387  {
388    return this->operator()(i);
389  }
390
391
392  bool vector::operator==(const vector& other) const
393  {
394    return equal(other);
395  }
396
397
398  bool vector::operator!=(const vector& other) const
399  {
400    return !equal(other);
401  }
402
403
404  double vector::operator*( const vector &other ) const
405  {
406    double res = 0.0;;
407    for ( size_t i = 0; i < size(); ++i ) 
408      res += other(i) * (*this)(i);
409    return res;
410  }
411
412
413  const vector& vector::operator=( const vector& other )
414  {
415    assert(v_);
416    if (gsl_vector_memcpy(v_,other.gsl_vector_p()))
417      throw utility::GSL_error("vector::set dimension mis-match");
418    return *this;
419  } 
420
421
422  const vector& vector::operator+=(const vector& other)
423  {
424    assert(v_);
425    int status=gsl_vector_add(v_, other.gsl_vector_p());
426    if (status)
427      throw utility::GSL_error(std::string("vector::add", status));
428    return *this;
429  }
430
431
432  const vector& vector::operator+=(double d)
433  {
434    assert(v_);
435    gsl_vector_add_constant(v_, d);
436    return *this;
437  }
438
439
440  const vector& vector::operator-=(const vector& other)
441  {
442    assert(v_);
443    int status=gsl_vector_sub(v_, other.gsl_vector_p());
444    if (status)
445      throw utility::GSL_error(std::string("vector::sub", status));
446    return *this;
447  }
448
449
450  const vector& vector::operator-=(const double d)
451  {
452    assert(v_);
453    gsl_vector_add_constant(v_, -d);
454    return *this;
455  }
456
457
458  const vector& vector::operator*=(const double d)
459  {
460    assert(v_);
461    gsl_vector_scale(v_, d);
462    return *this;
463  }
464
465
466  bool isnull(const vector& v)
467  {
468    return gsl_vector_isnull(v.gsl_vector_p());
469  }
470
471
472  double max(const vector& v)
473  {
474    return gsl_vector_max(v.gsl_vector_p());
475  }
476
477
478  size_t max_index(const vector& v)
479  {
480    return gsl_vector_max_index(v.gsl_vector_p());
481  }
482
483
484  double min(const vector& v)
485  {
486    return gsl_vector_min(v.gsl_vector_p());
487  }
488
489
490  size_t min_index(const vector& v)
491  {
492    return gsl_vector_min_index(v.gsl_vector_p());
493  }
494
495
496  bool nan(const vector& templat, vector& flag)
497  {
498    size_t vsize=templat.size();
499    if (vsize!=flag.size())
500      flag.clone(vector(vsize,1.0));
501    else
502      flag.all(1.0);
503    bool nan=false;
504    for (size_t i=0; i<vsize; i++)
505      if (std::isnan(templat(i))) {
506        flag(i)=0;
507        nan=true;
508      }
509    return nan;
510  }
511
512
513  void set_basis(vector& v, size_t i)
514  {
515    assert(v.gsl_vector_p());
516    gsl_vector_set_basis(v.gsl_vector_p(),i);
517  }
518
519
520  void shuffle(vector& invec)
521  {
522    vector vec(invec);
523    random::DiscreteUniform rnd;
524    for (size_t i=0; i<vec.size(); i++){
525      size_t index = rnd(vec.size()-i);
526      invec[i]=vec(index);
527      vec(index)=vec(vec.size()-i-1);
528    }
529  }
530
531
532  void sort(vector& v)
533  {
534    assert(v.gsl_vector_p());
535    gsl_sort_vector(v.gsl_vector_p());
536  }
537
538
539  void sort_index(std::vector<size_t>& sort_index, const vector& invec)
540  {
541    assert(invec.gsl_vector_p());
542    gsl_permutation* p = gsl_permutation_alloc(invec.size());
543    int status=gsl_sort_vector_index(p,invec.gsl_vector_p());
544    if (status) {
545      gsl_permutation_free(p);
546      throw utility::GSL_error(std::string("sort_index(vector&,const vector&)",status));     
547    }
548    sort_index=std::vector<size_t>(p->data,p->data+p->size);
549    gsl_permutation_free(p);
550  }
551
552
553  void sort_smallest_index(std::vector<size_t>& sort_index, size_t k, 
554                            const vector& invec)
555  {
556    assert(invec.gsl_vector_p());
557    assert(k<=invec.size());
558    sort_index.resize(k);
559    gsl_sort_vector_smallest_index(&sort_index[0],k,invec.gsl_vector_p());
560  }
561 
562  void sort_largest_index(std::vector<size_t>& sort_index, size_t k, 
563                            const vector& invec)
564  {
565    assert(invec.gsl_vector_p());
566    assert(k<=invec.size());
567    sort_index.resize(k);
568    gsl_sort_vector_largest_index(&sort_index[0],k,invec.gsl_vector_p());
569  }
570
571
572  double sum(const vector& v)
573  {
574    double sum = 0;
575    size_t vsize=v.size();
576    for (size_t i=0; i<vsize; ++i)
577      sum += v[i];
578    return sum;
579  }
580
581
582  void swap(vector& v, vector& w)
583  {
584    assert(v.gsl_vector_p()); assert(w.gsl_vector_p());
585    int status=gsl_vector_swap(v.gsl_vector_p(),w.gsl_vector_p());
586    if (status)
587      throw utility::GSL_error(std::string("swap(vector&,vector&)",status));
588  }
589
590
591  std::ostream& operator<<(std::ostream& s, const vector& a)
592  {
593    s.setf(std::ios::dec);
594    s.precision(12);
595    for (size_t j = 0; j < a.size(); ++j) {
596      s << a[j];
597      if ( (j+1)<a.size() )
598        s << s.fill();
599    }
600    return s;
601  }
602
603}}} // of namespace utility, yat, and thep
Note: See TracBrowser for help on using the repository browser.