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

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

going back to previous design in which view and const_view are in different classes. Having them in same didnt work as expected. There is a problem converting vector::iterator to vector::const_iterator. I'll open a separate ticket for this issue.

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