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

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

remove operator [] in Vector.

  • 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 1151 2008-02-25 22:32:04Z 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  bool VectorBase::operator==(const VectorBase& other) const
113  {
114    return equal(other);
115  }
116
117
118  bool VectorBase::operator!=(const VectorBase& other) const
119  {
120    return !equal(other);
121  }
122
123
124  double VectorBase::operator*( const VectorBase &other ) const
125  {
126    double res = 0.0;;
127    for ( size_t i = 0; i < size(); ++i ) 
128      res += other(i) * (*this)(i);
129    return res;
130  }
131
132
133  bool isnull(const VectorBase& v)
134  {
135    return gsl_vector_isnull(v.gsl_vector_p());
136  }
137
138
139  double max(const VectorBase& v)
140  {
141    return gsl_vector_max(v.gsl_vector_p());
142  }
143
144
145  size_t max_index(const VectorBase& v)
146  {
147    return gsl_vector_max_index(v.gsl_vector_p());
148  }
149
150
151  double min(const VectorBase& v)
152  {
153    return gsl_vector_min(v.gsl_vector_p());
154  }
155
156
157  size_t min_index(const VectorBase& v)
158  {
159    return gsl_vector_min_index(v.gsl_vector_p());
160  }
161
162
163  bool nan(const VectorBase& templat, Vector& flag)
164  {
165    size_t vsize(templat.size());
166    flag = Vector(vsize, 1.0);
167    bool nan=false;
168    for (size_t i=0; i<vsize; i++)
169      if (std::isnan(templat(i))) {
170        flag(i)=0;
171        nan=true;
172      }
173    return nan;
174  }
175
176
177  void sort_index(std::vector<size_t>& sort_index, const VectorBase& invec)
178  {
179    assert(invec.gsl_vector_p());
180    gsl_permutation* p = gsl_permutation_alloc(invec.size());
181    int status=gsl_sort_vector_index(p,invec.gsl_vector_p());
182    if (status) {
183      gsl_permutation_free(p);
184      throw utility::GSL_error(std::string("sort_index(vector&,const VectorBase&)",status));     
185    }
186    sort_index=std::vector<size_t>(p->data,p->data+p->size);
187    gsl_permutation_free(p);
188  }
189
190
191  void sort_smallest_index(std::vector<size_t>& sort_index, size_t k, 
192                            const VectorBase& invec)
193  {
194    assert(invec.gsl_vector_p());
195    assert(k<=invec.size());
196    sort_index.resize(k);
197    gsl_sort_vector_smallest_index(&sort_index[0],k,invec.gsl_vector_p());
198  }
199 
200  void sort_largest_index(std::vector<size_t>& sort_index, size_t k, 
201                            const VectorBase& invec)
202  {
203    assert(invec.gsl_vector_p());
204    assert(k<=invec.size());
205    sort_index.resize(k);
206    gsl_sort_vector_largest_index(&sort_index[0],k,invec.gsl_vector_p());
207  }
208
209
210  double sum(const VectorBase& v)
211  {
212    double sum = 0;
213    size_t vsize=v.size();
214    for (size_t i=0; i<vsize; ++i)
215      sum += v(i);
216    return sum;
217  }
218
219
220  std::ostream& operator<<(std::ostream& s, const VectorBase& a)
221  {
222    s.setf(std::ios::dec);
223    s.precision(12);
224    for (size_t j = 0; j < a.size(); ++j) {
225      s << a(j);
226      if ( (j+1)<a.size() )
227        s << s.fill();
228    }
229    return s;
230  }
231
232}}} // of namespace utility, yat, and thep
Note: See TracBrowser for help on using the repository browser.