source: trunk/yat/utility/VectorBase.cc @ 1121

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

fixes #308

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