source: trunk/yat/utility/vectorBase.cc @ 1009

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

merging branch peter-dev into trunk delta 1008:994

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.8 KB
Line 
1// $Id: vectorBase.cc 1009 2008-02-01 04:15:44Z 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(void)
46    : vec_(NULL), const_vec_(NULL)
47  {
48  }
49
50
51  vectorBase::vectorBase(gsl_vector* v)
52    : vec_(v), const_vec_(v)
53  {
54  }
55
56
57  vectorBase::vectorBase(const gsl_vector* v)
58    : vec_(NULL), const_vec_(v)
59  {
60  }
61
62
63  vectorBase::~vectorBase(void)
64  {
65  }
66
67
68  void vectorBase::all(double value)
69  {
70    assert(vec_);
71    gsl_vector_set_all(vec_,value);
72  }
73
74
75  vectorBase::const_iterator vectorBase::begin(void) const
76  {
77    return const_iterator(*this, 0);
78  }
79
80
81  vectorBase::iterator vectorBase::begin(void)
82  {
83    return iterator(*this, 0);
84  }
85
86
87  void vectorBase::div(const vectorBase& other)
88  {
89    assert(vec_);
90    int status=gsl_vector_div(vec_,other.gsl_vector_p());
91    if (status)
92      throw utility::GSL_error(std::string("vectorBase::div",status));
93  }
94
95
96  vectorBase::const_iterator vectorBase::end(void) const
97  {
98    return const_iterator(*this, size());
99  }
100
101
102  vectorBase::iterator vectorBase::end(void)
103  {
104    return iterator(*this, size());
105  }
106
107
108  bool vectorBase::equal(const vectorBase& other, const double d) const
109  {
110    if (this==&other)
111      return true;
112    if (size()!=other.size())
113      return false;
114    // if gsl error handler disabled, out of bounds index will not
115    // abort the program.
116    for (size_t i=0; i<size(); ++i)
117      if (fabs( (*this)(i)-other(i) ) > d ||
118          std::isnan((*this)(i)) || std::isnan(other(i)) )
119        return false;
120    return true;
121  }
122
123
124  const gsl_vector* vectorBase::gsl_vector_p(void) const
125  {
126    return const_vec_;
127  }
128
129
130  gsl_vector* vectorBase::gsl_vector_p(void)
131  {
132    return vec_;
133  }
134
135
136  void vectorBase::mul(const vectorBase& other)
137  {
138    assert(vec_);
139    int status=gsl_vector_mul(vec_,other.gsl_vector_p());
140    if (status)
141      throw utility::GSL_error(std::string("vectorBase::div",status));
142  }
143
144
145  size_t vectorBase::size(void) const
146  {
147    if (!const_vec_)
148      return 0;
149    return const_vec_->size;
150  }
151
152
153  const double& vectorBase::operator()(size_t i) const
154  {
155    const double* d=gsl_vector_const_ptr(const_vec_, i);
156    if (!d)
157      throw utility::GSL_error("vectorBase::operator()",GSL_EINVAL);
158    return *d;
159  }
160
161
162  double& vectorBase::operator()(size_t i)
163  {
164    double* d=gsl_vector_ptr(vec_, i);
165    if (!d)
166      throw utility::GSL_error("vectorBase::operator()",GSL_EINVAL);
167    return *d;
168  }
169
170
171  double& vectorBase::operator[](size_t i)
172  {
173    return this->operator()(i);
174  }
175
176
177  const double& vectorBase::operator[](size_t i) const
178  {
179    return this->operator()(i);
180  }
181
182
183  bool vectorBase::operator==(const vectorBase& other) const
184  {
185    return equal(other);
186  }
187
188
189  bool vectorBase::operator!=(const vectorBase& other) const
190  {
191    return !equal(other);
192  }
193
194
195  double vectorBase::operator*( const vectorBase &other ) const
196  {
197    double res = 0.0;;
198    for ( size_t i = 0; i < size(); ++i ) 
199      res += other(i) * (*this)(i);
200    return res;
201  }
202
203
204  const vectorBase& vectorBase::operator+=(double d)
205  {
206    assert(vec_);
207    gsl_vector_add_constant(vec_, d);
208    return *this;
209  }
210
211
212  const vectorBase& vectorBase::operator-=(const vectorBase& other)
213  {
214    assert(vec_);
215    int status=gsl_vector_sub(vec_, other.gsl_vector_p());
216    if (status)
217      throw utility::GSL_error(std::string("vectorBase::sub", status));
218    return *this;
219  }
220
221
222  const vectorBase& vectorBase::operator-=(const double d)
223  {
224    assert(vec_);
225    gsl_vector_add_constant(vec_, -d);
226    return *this;
227  }
228
229
230  const vectorBase& vectorBase::operator*=(double d)
231  {
232    assert(vec_);
233    gsl_vector_scale(vec_, d);
234    return *this;
235  }
236
237
238  bool isnull(const vectorBase& v)
239  {
240    return gsl_vector_isnull(v.gsl_vector_p());
241  }
242
243
244  double max(const vectorBase& v)
245  {
246    return gsl_vector_max(v.gsl_vector_p());
247  }
248
249
250  size_t max_index(const vectorBase& v)
251  {
252    return gsl_vector_max_index(v.gsl_vector_p());
253  }
254
255
256  double min(const vectorBase& v)
257  {
258    return gsl_vector_min(v.gsl_vector_p());
259  }
260
261
262  size_t min_index(const vectorBase& v)
263  {
264    return gsl_vector_min_index(v.gsl_vector_p());
265  }
266
267
268  bool nan(const vectorBase& templat, vector& flag)
269  {
270    size_t vsize(templat.size());
271    flag = vector(vsize, 1.0);
272    bool nan=false;
273    for (size_t i=0; i<vsize; i++)
274      if (std::isnan(templat(i))) {
275        flag(i)=0;
276        nan=true;
277      }
278    return nan;
279  }
280
281
282  void shuffle(vectorBase& invec)
283  {
284    random::DiscreteUniform rnd;
285    std::random_shuffle(invec.begin(), invec.end(), rnd);
286  }
287
288
289  void sort_index(std::vector<size_t>& sort_index, const vectorBase& invec)
290  {
291    assert(invec.gsl_vector_p());
292    gsl_permutation* p = gsl_permutation_alloc(invec.size());
293    int status=gsl_sort_vector_index(p,invec.gsl_vector_p());
294    if (status) {
295      gsl_permutation_free(p);
296      throw utility::GSL_error(std::string("sort_index(vector&,const vectorBase&)",status));     
297    }
298    sort_index=std::vector<size_t>(p->data,p->data+p->size);
299    gsl_permutation_free(p);
300  }
301
302
303  void sort_smallest_index(std::vector<size_t>& sort_index, size_t k, 
304                            const vectorBase& invec)
305  {
306    assert(invec.gsl_vector_p());
307    assert(k<=invec.size());
308    sort_index.resize(k);
309    gsl_sort_vector_smallest_index(&sort_index[0],k,invec.gsl_vector_p());
310  }
311 
312  void sort_largest_index(std::vector<size_t>& sort_index, size_t k, 
313                            const vectorBase& invec)
314  {
315    assert(invec.gsl_vector_p());
316    assert(k<=invec.size());
317    sort_index.resize(k);
318    gsl_sort_vector_largest_index(&sort_index[0],k,invec.gsl_vector_p());
319  }
320
321
322  double sum(const vectorBase& v)
323  {
324    double sum = 0;
325    size_t vsize=v.size();
326    for (size_t i=0; i<vsize; ++i)
327      sum += v(i);
328    return sum;
329  }
330
331
332  std::ostream& operator<<(std::ostream& s, const vectorBase& a)
333  {
334    s.setf(std::ios::dec);
335    s.precision(12);
336    for (size_t j = 0; j < a.size(); ++j) {
337      s << a(j);
338      if ( (j+1)<a.size() )
339        s << s.fill();
340    }
341    return s;
342  }
343
344}}} // of namespace utility, yat, and thep
Note: See TracBrowser for help on using the repository browser.