source: trunk/src/vector.cc @ 125

Last change on this file since 125 was 125, checked in by Peter, 17 years ago

shuffle function added

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.2 KB
Line 
1// $Id: vector.cc 125 2004-08-02 15:15:34Z peter $
2
3#include "vector.h"
4#include "random_singleton.h"
5
6#include <iostream>
7#include <sstream>
8#include <string>
9#include <vector>
10#include <utility>
11
12
13
14
15
16namespace theplu {
17namespace gslapi {
18
19
20
21  vector::vector(void)
22    : v_(NULL)
23  {
24  }
25
26
27
28  vector::vector(const size_t n, const double init_value)
29  {
30    v_ = gsl_vector_alloc(n);
31    set_all(init_value);
32  }
33
34
35
36  vector::vector(const vector& other, const std::vector<size_t>& index)
37  {
38    if (index.size()){
39      v_ = gsl_vector_alloc(index.size());
40      for (size_t i=0; i<index.size(); i++) 
41        gsl_vector_set(v_, i, other(index[i]));
42    }
43    else
44      v_ = other.gsl_vector_copy();
45  }
46
47
48
49  vector::vector(gsl_vector* vec)
50    : v_(vec)
51  {
52  }
53
54
55
56  vector::vector(std::istream& is)
57  {
58    using namespace std;
59 
60    // read the data file and store in stl vectors (dynamically expandable)
61    std::vector<std::vector<double> > data_matrix;
62    u_int nof_columns=0;
63    u_int nof_rows = 0;
64    string s;
65    for (nof_rows = 0; getline(is, s, '\n'); nof_rows++) {
66      istringstream line(s);
67      std::vector<double> v;
68      string tmp_string;
69      while (line >> tmp_string) {
70        if(!is.good()) {
71          // Jari. change to throw exception, remove unnecessary includes
72          cerr << "vector::vector(std::istream&): "
73               << "error reading data file!" << endl;
74          exit(1);
75        }
76        double t=atof(tmp_string.c_str());
77        // Here we should check that it actually was a double!!!!!
78        v.push_back(t);
79      }
80      if(nof_columns==0)
81        nof_columns=v.size();
82      else if(v.size()!=nof_columns) {
83        // Jari. change to throw exception, remove unnecessary includes
84        cerr << "vector::vector(std::istream&) data file error: "
85             << "line" << nof_rows+1 << " has " << v.size() 
86             << " columns; expected " << nof_columns
87             << " columns"
88             << endl; 
89        exit(1);
90      }
91      data_matrix.push_back(v);
92    }
93 
94    // manipulate the state of the stream to be good
95    is.clear(std::ios::goodbit);
96 
97    // convert the data to a gsl vector and check that data file is a
98    // column vector or a row vector
99    if(nof_columns==1) {
100      v_ = gsl_vector_alloc ( nof_rows );
101      for(u_int i=0;i<nof_rows;i++) 
102        gsl_vector_set( v_, i, data_matrix[i][0] );
103    }
104    else if(nof_rows==1){
105      v_ = gsl_vector_alloc ( nof_columns );
106      for(u_int i=0;i<nof_columns;i++) 
107        gsl_vector_set( v_, i, data_matrix[0][i] );
108    }
109    else {
110      // Jari. change to throw exception, remove unnecessary includes
111      cerr << "vector::vector(std::istream&) data file error: "
112           << "file has " << nof_rows << " rows and " << nof_columns
113           << " columns; expected a row vector or a column vector"
114           << endl;
115      exit(1);
116    }
117  }
118
119
120
121  vector::~vector(void)
122  {
123    if (v_) {
124      gsl_vector_free(v_);
125      v_=NULL;
126    }
127  }
128
129
130
131  gsl_vector* vector::gsl_vector_copy(void) const
132  {
133    gsl_vector* vec = gsl_vector_alloc(size());
134    gsl_vector_memcpy(vec,v_);
135    return vec;
136  }
137
138
139
140
141
142  std::pair<size_t,size_t> vector::minmax_index(void) const
143  {
144    size_t min_index=0;
145    size_t max_index=0;
146    void gsl_vector_minmax_index (const gsl_vector * v_, size_t * 
147                                  min_index, size_t * max_index);
148    return std::pair<size_t,size_t>(min_index, max_index);
149  }
150
151
152
153
154  std::pair<size_t,size_t> vector::minmax_index(const std::vector<size_t>& subset ) const
155  {
156    size_t min_index = subset[0];
157    size_t max_index = subset[0];
158    for (unsigned int i=1; i<subset.size(); i++) {
159      if (gsl_vector_get( v_, subset[i]) < gsl_vector_get( v_, min_index))
160        min_index = subset[i];
161      else if (gsl_vector_get( v_, subset[i]) > gsl_vector_get( v_, max_index))
162        max_index = subset[i];
163    }
164    return std::pair<size_t,size_t>(min_index, max_index);
165  }
166
167
168
169  vector vector::mul_elements( const vector& other ) const
170  {
171    vector res( *this );
172    gsl_vector_mul( res.v_, other.v_ );
173    return res;
174  }
175
176  vector vector::shuffle(vector& vec) const
177  {
178    vector tmp_vec(vec);
179    theplu::cpptools::random_singleton* 
180      rnd=theplu::cpptools::random_singleton::get_instance();
181    for (size_t i=0; i<vec.size(); i++){
182      size_t index = rnd->get_uniform_int(vec.size()-i);
183      vec(i)=tmp_vec(index);
184      tmp_vec(index)=tmp_vec(vec.size()-i-1);
185    }
186   
187    return vec;
188  }
189
190  double vector::sum(void) const
191  {
192    double sum = 0;
193    for (u_int i=0; i<size(); i++)
194      sum += gsl_vector_get( v_, i );
195    return( sum );
196  } 
197
198
199
200  vector vector::operator-(void) const
201  {
202    vector res( *this );
203    gsl_vector_scale (res.v_,-1.);
204    return res;
205  }
206
207
208
209  double vector::operator*( const vector &other ) const
210  {
211    // Jari, check for gsl_support
212    double res = 0.0;;
213    for ( size_t i = 0; i < size(); ++i ) 
214      res += other[i] * (*this)[i];
215    return res;
216  }
217
218
219
220  vector vector::operator+(const vector &other) const
221  {
222    vector res(*this);
223    gsl_vector_add(res.v_,other.v_);
224    return res;
225  }
226
227
228
229  vector vector::operator-( const vector &other ) const
230  {
231    vector res( *this );
232    gsl_vector_sub(res.v_,other.v_);
233    return res;
234  }
235
236
237
238  vector& vector::operator=( const vector& other )
239  {
240    if( this != &other ) {
241      if ( v_ )
242        gsl_vector_free( v_ );
243      v_ = other.gsl_vector_copy();
244    }
245    return *this;
246  } 
247
248
249
250  std::ostream& theplu::gslapi::operator<<(std::ostream& s, const vector& a)
251  {
252    s.setf(std::ios::fixed);
253
254    for (size_t j = 0; j < a.size(); ++j) {
255      s << a[j];
256      if ( (j+1)<a.size() )
257        s << " ";
258    }
259
260    return s;
261  }
262
263
264}} // of namespace gslapi and namespace thep
Note: See TracBrowser for help on using the repository browser.