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

Last change on this file since 1566 was 1519, checked in by Peter, 13 years ago

speeding up QuantileNormalizer? by avoiding some copying and by only sorting each column once and not twice as before.

  • 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 1519 2008-09-21 04:35:39Z 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#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    std::vector<size_t> tmp(p->data,p->data+p->size);
185    sort_index.swap(tmp);
186    gsl_permutation_free(p);
187  }
188
189
190  void sort_smallest_index(std::vector<size_t>& sort_index, size_t k, 
191                            const VectorBase& invec)
192  {
193    assert(invec.gsl_vector_p());
194    assert(k<=invec.size());
195    sort_index.resize(k);
196    gsl_sort_vector_smallest_index(&sort_index[0],k,invec.gsl_vector_p());
197  }
198 
199  void sort_largest_index(std::vector<size_t>& sort_index, size_t k, 
200                            const VectorBase& invec)
201  {
202    assert(invec.gsl_vector_p());
203    assert(k<=invec.size());
204    sort_index.resize(k);
205    gsl_sort_vector_largest_index(&sort_index[0],k,invec.gsl_vector_p());
206  }
207
208
209  double sum(const VectorBase& v)
210  {
211    double sum = 0;
212    size_t vsize=v.size();
213    for (size_t i=0; i<vsize; ++i)
214      sum += v(i);
215    return sum;
216  }
217
218
219  std::ostream& operator<<(std::ostream& s, const VectorBase& a)
220  {
221    s.setf(std::ios::dec);
222    s.precision(12);
223    for (size_t j = 0; j < a.size(); ++j) {
224      s << a(j);
225      if ( (j+1)<a.size() )
226        s << s.fill();
227    }
228    return s;
229  }
230
231}}} // of namespace utility, yat, and thep
Note: See TracBrowser for help on using the repository browser.