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

Last change on this file since 1487 was 1487, checked in by Jari Häkkinen, 13 years ago

Addresses #436. GPL license copy reference should also be updated.

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