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

Last change on this file since 1706 was 1682, checked in by Peter, 13 years ago

prefer std::abs and re-use code

  • 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 1682 2008-12-29 21:27:32Z 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, 2007 Jari Häkkinen, Peter Johansson, Markus Ringnér
7  Copyright (C) 2008 Peter Johansson
8
9  This file is part of the yat library, http://dev.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 3 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 yat. If not, see <http://www.gnu.org/licenses/>.
23*/
24
25#include "VectorBase.h"
26#include "Vector.h"
27#include "utility.h"
28
29#include <gsl/gsl_sort_vector.h>
30
31#include <algorithm>
32#include <cassert>
33#include <cmath>
34#include <ostream>
35#include <string>
36#include <utility>
37#include <vector>
38
39namespace theplu {
40namespace yat {
41namespace utility {
42
43
44  VectorBase::VectorBase(const gsl_vector* v)
45    : const_vec_(v)
46  {
47  }
48
49
50  VectorBase::~VectorBase(void)
51  {
52  }
53
54
55  VectorBase::const_iterator VectorBase::begin(void) const
56  {
57    if (const_vec_)
58      return const_iterator(&(this->operator()(0)), const_vec_->stride);
59    return const_iterator(NULL, 1);
60  }
61
62
63  VectorBase::const_iterator VectorBase::end(void) const
64  {
65    if (const_vec_)
66      return const_iterator(&(this->operator()(0))+const_vec_->stride*size(), 
67                            const_vec_->stride);
68    return const_iterator(NULL, 1);
69  }
70
71
72  bool VectorBase::equal(const VectorBase& other, const double d) const
73  {
74    if (this==&other)
75      return true;
76    if (size()!=other.size())
77      return false;
78    // if gsl error handler disabled, out of bounds index will not
79    // abort the program.
80    for (size_t i=0; i<size(); ++i)
81      if (std::abs( (*this)(i)-other(i) ) > d ||
82          std::isnan((*this)(i)) || std::isnan(other(i)) )
83        return false;
84    return true;
85  }
86
87
88  const gsl_vector* VectorBase::gsl_vector_p(void) const
89  {
90    return const_vec_;
91  }
92
93
94  size_t VectorBase::size(void) const
95  {
96    if (!const_vec_)
97      return 0;
98    return const_vec_->size;
99  }
100
101
102  const double& VectorBase::operator()(size_t i) const
103  {
104    const double* d=gsl_vector_const_ptr(const_vec_, i);
105    if (!d)
106      throw utility::GSL_error("VectorBase::operator()",GSL_EINVAL);
107    return *d;
108  }
109
110
111  bool VectorBase::operator==(const VectorBase& other) const
112  {
113    return equal(other);
114  }
115
116
117  bool VectorBase::operator!=(const VectorBase& other) const
118  {
119    return !equal(other);
120  }
121
122
123  double VectorBase::operator*( const VectorBase &other ) const
124  {
125    double res = 0.0;;
126    for ( size_t i = 0; i < size(); ++i ) 
127      res += other(i) * (*this)(i);
128    return res;
129  }
130
131
132  bool isnull(const VectorBase& v)
133  {
134    return gsl_vector_isnull(v.gsl_vector_p());
135  }
136
137
138  double max(const VectorBase& v)
139  {
140    return gsl_vector_max(v.gsl_vector_p());
141  }
142
143
144  size_t max_index(const VectorBase& v)
145  {
146    return gsl_vector_max_index(v.gsl_vector_p());
147  }
148
149
150  double min(const VectorBase& v)
151  {
152    return gsl_vector_min(v.gsl_vector_p());
153  }
154
155
156  size_t min_index(const VectorBase& v)
157  {
158    return gsl_vector_min_index(v.gsl_vector_p());
159  }
160
161
162  bool nan(const VectorBase& templat, Vector& flag)
163  {
164    size_t vsize(templat.size());
165    flag = Vector(vsize, 1.0);
166    bool nan=false;
167    for (size_t i=0; i<vsize; i++)
168      if (std::isnan(templat(i))) {
169        flag(i)=0;
170        nan=true;
171      }
172    return nan;
173  }
174
175
176  void sort_index(std::vector<size_t>& sort_index, const VectorBase& invec)
177  {
178    assert(invec.gsl_vector_p());
179    gsl_permutation* p = gsl_permutation_alloc(invec.size());
180    int status=gsl_sort_vector_index(p,invec.gsl_vector_p());
181    if (status) {
182      gsl_permutation_free(p);
183      throw utility::GSL_error(std::string("sort_index(vector&,const VectorBase&)",status));     
184    }
185    std::vector<size_t> tmp(p->data,p->data+p->size);
186    sort_index.swap(tmp);
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.