source: branches/peter-dev/yat/utility/vectorBase.cc @ 995

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

Splitting vector class into three different classes for vector, view and const_view, and adding abstract classes appropriately. NOTE, the tests are not going through, not even compiling; sorry about that

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.0 KB
Line 
1// $Id: vectorBase.cc 995 2007-11-30 22:55:30Z 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 "vectorBase.h"
28#include "matrix.h"
29#include "utility.h"
30#include "yat/random/random.h"
31
32#include <algorithm>
33#include <cassert>
34#include <cmath>
35#include <iostream>
36#include <sstream>
37#include <utility>
38#include <vector>
39
40namespace theplu {
41namespace yat {
42namespace utility {
43
44
45  vectorBase::vectorBase(const gsl_vector* v)
46    : const_vec_(v)
47  {
48  }
49
50
51  vectorBase::~vectorBase(void)
52  {
53  }
54
55
56  vectorBase::const_iterator vectorBase::begin(void) const
57  {
58    return vectorBase::const_iterator(*this, 0);
59  }
60
61
62  vectorBase::const_iterator vectorBase::end(void) const
63  {
64    return vectorBase::const_iterator(*this, size());
65  }
66
67
68  bool vectorBase::equal(const vectorBase& other, const double d) const
69  {
70    if (this==&other)
71      return true;
72    if (size()!=other.size())
73      return false;
74    // if gsl error handler disabled, out of bounds index will not
75    // abort the program.
76    for (size_t i=0; i<size(); ++i)
77      if (fabs( (*this)(i)-other(i) ) > d ||
78          std::isnan((*this)(i)) || std::isnan(other(i)) )
79        return false;
80    return true;
81  }
82
83
84  const gsl_vector* vectorBase::gsl_vector_p(void) const
85  {
86    return const_vec_;
87  }
88
89
90  size_t vectorBase::size(void) const
91  {
92    if (const_vec_)
93      return 0;
94    return const_vec_->size;
95  }
96
97
98  const double& vectorBase::operator()(size_t i) const
99  {
100    const double* d=gsl_vector_const_ptr(const_vec_, i);
101    if (!d)
102      throw utility::GSL_error("vectorBase::operator()",GSL_EINVAL);
103    return *d;
104  }
105
106
107  const double& vectorBase::operator[](size_t i) const
108  {
109    const double* d=gsl_vector_const_ptr(const_vec_, i);
110    if (!d)
111      throw utility::GSL_error("vectorBase::operator[]",GSL_EINVAL);
112    return *d;
113  }
114
115
116  bool vectorBase::operator==(const vectorBase& other) const
117  {
118    return equal(other);
119  }
120
121
122  bool vectorBase::operator!=(const vectorBase& other) const
123  {
124    return !equal(other);
125  }
126
127
128  double vectorBase::operator*( const vectorBase &other ) const
129  {
130    double res = 0.0;;
131    for ( size_t i = 0; i < size(); ++i ) 
132      res += other(i) * (*this)(i);
133    return res;
134  }
135
136
137  bool isnull(const vectorBase& v)
138  {
139    return gsl_vector_isnull(v.gsl_vector_p());
140  }
141
142
143  double max(const vectorBase& v)
144  {
145    return gsl_vector_max(v.gsl_vector_p());
146  }
147
148
149  size_t max_index(const vectorBase& v)
150  {
151    return gsl_vector_max_index(v.gsl_vector_p());
152  }
153
154
155  double min(const vectorBase& v)
156  {
157    return gsl_vector_min(v.gsl_vector_p());
158  }
159
160
161  size_t min_index(const vectorBase& v)
162  {
163    return gsl_vector_min_index(v.gsl_vector_p());
164  }
165
166
167  bool nan(const vectorBase& templat, vector& flag)
168  {
169    size_t vsize(templat.size());
170    flag = vector(vsize, 1.0);
171    bool nan=false;
172    for (size_t i=0; i<vsize; i++)
173      if (std::isnan(templat(i))) {
174        flag(i)=0;
175        nan=true;
176      }
177    return nan;
178  }
179
180
181  void sort_index(std::vector<size_t>& sort_index, const vectorBase& invec)
182  {
183    assert(invec.gsl_vector_p());
184    gsl_permutation* p = gsl_permutation_alloc(invec.size());
185    int status=gsl_sort_vector_index(p,invec.gsl_vector_p());
186    if (status) {
187      gsl_permutation_free(p);
188      throw utility::GSL_error(std::string("sort_index(vector&,const vectorBase&)",status));     
189    }
190    sort_index=std::vector<size_t>(p->data,p->data+p->size);
191    gsl_permutation_free(p);
192  }
193
194
195  void sort_smallest_index(std::vector<size_t>& sort_index, size_t k, 
196                            const vectorBase& invec)
197  {
198    assert(invec.gsl_vector_p());
199    assert(k<=invec.size());
200    sort_index.resize(k);
201    gsl_sort_vector_smallest_index(&sort_index[0],k,invec.gsl_vector_p());
202  }
203 
204  void sort_largest_index(std::vector<size_t>& sort_index, size_t k, 
205                            const vectorBase& invec)
206  {
207    assert(invec.gsl_vector_p());
208    assert(k<=invec.size());
209    sort_index.resize(k);
210    gsl_sort_vector_largest_index(&sort_index[0],k,invec.gsl_vector_p());
211  }
212
213
214  double sum(const vectorBase& v)
215  {
216    double sum = 0;
217    size_t vsize=v.size();
218    for (size_t i=0; i<vsize; ++i)
219      sum += v(i);
220    return sum;
221  }
222
223
224  std::ostream& operator<<(std::ostream& s, const vectorBase& a)
225  {
226    s.setf(std::ios::dec);
227    s.precision(12);
228    for (size_t j = 0; j < a.size(); ++j) {
229      s << a(j);
230      if ( (j+1)<a.size() )
231        s << s.fill();
232    }
233    return s;
234  }
235
236}}} // of namespace utility, yat, and thep
Note: See TracBrowser for help on using the repository browser.