source: trunk/yat/utility/vector.cc @ 899

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

changed euclidean to fulfill design criteria, removes sum_squared_deviation from AveragerPairWeighted?, also changed vector::sort to use stl rather than gsl. refs #248 and fixes #253

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.1 KB
Line 
1// $Id: vector.cc 899 2007-09-26 21:26:25Z 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 "vector.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  vector::vector(void)
46    : v_(NULL), view_(NULL), view_const_(NULL), proxy_v_(NULL)
47  {
48  }
49
50
51  vector::vector(size_t n, double init_value)
52    : v_(gsl_vector_alloc(n)), view_(NULL), view_const_(NULL),
53      proxy_v_(v_)
54  {
55    if (!v_)
56      throw utility::GSL_error("vector::vector failed to allocate memory");
57
58    all(init_value);
59  }
60
61
62  vector::vector(const vector& other)
63    : v_(other.create_gsl_vector_copy()), view_(NULL),
64      view_const_(NULL), proxy_v_(v_)
65  {
66  }
67
68
69  vector::vector(vector& v, size_t offset, size_t n, size_t stride)
70    : view_const_(NULL)
71  {
72    view_ = new gsl_vector_view(gsl_vector_subvector_with_stride(v.v_,offset,
73                                                                 stride,n));
74    if (!view_)
75      throw utility::GSL_error("vector::vector failed to setup view");
76    proxy_v_ = v_ = &(view_->vector);
77  }
78
79
80  vector::vector(const vector& v, size_t offset, size_t n, size_t stride)
81    : v_(NULL), view_(NULL)
82  {
83    view_const_ = new gsl_vector_const_view(
84                   gsl_vector_const_subvector_with_stride(v.v_,offset,stride,n));
85    if (!view_const_)
86      throw utility::GSL_error("vector::vector failed to setup view");
87    proxy_v_ = &(view_const_->vector);
88  }
89
90
91  vector::vector(matrix& m, size_t i, bool row)
92    : view_const_(NULL)
93  {
94    view_=new gsl_vector_view(row ?
95                              gsl_matrix_row   (m.gsl_matrix_p(),i) :
96                              gsl_matrix_column(m.gsl_matrix_p(),i)  );
97    if (!view_)
98      throw utility::GSL_error("vector::vector failed to setup view");
99    proxy_v_ = v_ = &(view_->vector);
100  }
101
102
103  vector::vector(const matrix& m, size_t i, bool row)
104    : v_(NULL), view_(NULL)
105  {
106    view_const_ = new gsl_vector_const_view(row ?
107                                   gsl_matrix_const_row   (m.gsl_matrix_p(),i) :
108                                   gsl_matrix_const_column(m.gsl_matrix_p(),i) );
109    if (!view_const_)
110      throw utility::GSL_error("vector::vector failed to setup view");
111    proxy_v_ = &(view_const_->vector);
112  }
113
114
115  vector::vector(std::istream& is, char sep)
116    throw (utility::IO_error, std::exception)
117    : view_(NULL), view_const_(NULL)
118  {
119    // read the data file and store in stl vectors (dynamically
120    // expandable)
121    std::vector<std::vector<double> > data_matrix;
122    u_int nof_columns=0;
123    u_int nof_rows=0;
124    std::string line;
125    while(getline(is, line, '\n')) {
126      // Empty lines
127      if (!line.size())
128          continue;
129      nof_rows++;
130
131      std::vector<double> v;
132      std::string element;
133      std::stringstream ss(line);
134      bool ok=true;
135      while(ok) {
136        if(sep=='\0')
137          ok=(ss>>element);
138        else 
139          ok=getline(ss, element, sep);
140        if(!ok)
141          break;       
142
143        if(utility::is_double(element)) {
144          v.push_back(atof(element.c_str()));
145        }
146        else if (!element.size() || utility::is_nan(element)) {
147          v.push_back(std::numeric_limits<double>::quiet_NaN());
148        }
149        else {
150          std::stringstream ss("Warning: '");
151          ss << element << "' is not a double.";
152          throw IO_error(ss.str());
153        }
154      } 
155      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
156          v.push_back(std::numeric_limits<double>::quiet_NaN());
157      if (!nof_columns)
158        nof_columns=v.size();
159      else if (nof_rows && (nof_columns>1)) {
160        std::ostringstream s;
161        s << "vector::vector(std::istream&) data file error:\n"
162          << "    File has inconsistent number of rows (" << nof_rows
163          << ") and columns (" << nof_columns
164          << ").\n    Expected a row or a column vector.";
165        throw utility::IO_error(s.str());
166      }
167      else if (v.size()!=nof_columns) {
168        std::ostringstream s;
169        s << "vector::vector(std::istream&) data file error:\n"
170          << "    Line " << nof_rows << " has " << v.size()
171          << " columns; expected " << nof_columns << " column.";
172        throw utility::IO_error(s.str());
173      }
174      data_matrix.push_back(v);
175    }
176
177    // manipulate the state of the stream to be good
178    is.clear(std::ios::goodbit);
179    // convert the data to a gsl vector
180    proxy_v_ = v_ = gsl_vector_alloc(nof_rows*nof_columns);
181    if (!v_)
182      throw utility::GSL_error("vector::vector failed to allocate memory");
183    size_t n=0;
184    // if gsl error handler disabled, out of bounds index will not
185    // abort the program.
186    for (size_t i=0; i<nof_rows; i++)
187      for (size_t j=0; j<nof_columns; j++) 
188        gsl_vector_set( v_, n++, data_matrix[i][j] );
189  }
190
191
192  vector::~vector(void)
193  {
194    delete_allocated_memory();
195  }
196
197
198  vector::iterator vector::begin(void)
199  {
200    return vector::iterator(*this, 0);
201  }
202
203
204  vector::const_iterator vector::begin(void) const
205  {
206    return vector::const_iterator(*this, 0);
207  }
208
209
210  const vector& vector::clone(const vector& other)
211  {
212    if (this!=&other) {
213
214      delete_allocated_memory();
215
216      if (other.view_) {
217        view_ = new gsl_vector_view(*other.view_);
218        proxy_v_ = v_ = &(view_->vector);
219      }
220      else if (other.view_const_) {
221        view_const_ = new gsl_vector_const_view(*other.view_const_);
222        proxy_v_ = &(view_const_->vector);
223        v_=NULL;
224      }
225      else if (other.v_)
226        proxy_v_ = v_ = other.create_gsl_vector_copy();
227
228    }
229    return *this;
230  } 
231
232
233  gsl_vector* vector::create_gsl_vector_copy(void) const
234  {
235    gsl_vector* vec = gsl_vector_alloc(size());
236    if (!vec)
237      throw utility::GSL_error("vector::create_gsl_vector_copy failed to allocate memory");
238    if (gsl_vector_memcpy(vec, proxy_v_))
239      throw utility::GSL_error("vector::create_gsl_matrix_copy dimension mis-match");
240    return vec;
241  }
242
243
244  void vector::delete_allocated_memory(void)
245  {
246    if (view_)
247      delete view_;
248    else if (view_const_)
249      delete view_const_;
250    else if (v_)
251      gsl_vector_free(v_);
252    proxy_v_=v_=NULL;
253  }
254
255
256  void vector::div(const vector& other)
257  {
258    assert(v_);
259    int status=gsl_vector_div(v_,other.gsl_vector_p());
260    if (status)
261      throw utility::GSL_error(std::string("vector::div",status));
262  }
263
264
265  vector::iterator vector::end(void)
266  {
267    return vector::iterator(*this, size());
268  }
269
270
271  vector::const_iterator vector::end(void) const
272  {
273    return vector::const_iterator(*this, size());
274  }
275
276
277  bool vector::equal(const vector& other, const double d) const
278  {
279    if (this==&other)
280      return true;
281    if (size()!=other.size())
282      return false;
283    // if gsl error handler disabled, out of bounds index will not
284    // abort the program.
285    for (size_t i=0; i<size(); ++i)
286      // The two last condition checks are needed for NaN detection
287      if (fabs( (*this)(i)-other(i) ) > d ||
288          (*this)(i)!=(*this)(i) || other(i)!=other(i))
289        return false;
290    return true;
291  }
292
293
294  const gsl_vector* vector::gsl_vector_p(void) const
295  {
296    return proxy_v_;
297  }
298
299
300  gsl_vector* vector::gsl_vector_p(void)
301  {
302    return v_;
303  }
304
305
306  bool vector::isview(void) const
307  {
308    return view_ || view_const_;
309  }
310
311
312  void vector::mul(const vector& other)
313  {
314    assert(v_);
315    int status=gsl_vector_mul(v_,other.gsl_vector_p());
316    if (status)
317      throw utility::GSL_error(std::string("vector::mul",status));
318  }
319
320
321  void vector::resize(size_t n, double init_value)
322  {
323    delete_allocated_memory();
324    proxy_v_ = v_ = gsl_vector_alloc(n);
325    if (!v_)
326      throw utility::GSL_error("vector::vector failed to allocate memory");
327    all(init_value);
328  } 
329
330
331  void vector::reverse(void)
332  {
333    assert(v_);
334    gsl_vector_reverse(v_);
335  }
336
337
338  void vector::all(const double& value)
339  {
340    assert(v_);
341    gsl_vector_set_all(v_,value);
342  }
343
344
345  size_t vector::size(void) const
346  {
347    if (!proxy_v_)
348      return 0;
349    return proxy_v_->size;
350  }
351
352
353  void vector::swap(size_t i, size_t j)
354  {
355    assert(v_);
356    int status=gsl_vector_swap_elements(v_, i, j);
357    if (status)
358      throw utility::GSL_error(std::string("vector::swap_elements",status));
359  }
360
361
362  double& vector::operator()(size_t i)
363  {
364    assert(v_);
365    double* d=gsl_vector_ptr(v_, i);
366    if (!d)
367      throw utility::GSL_error("vector::operator()",GSL_EINVAL);
368    return *d;
369  }
370
371
372  const double& vector::operator()(size_t i) const
373  {
374    const double* d=gsl_vector_const_ptr(proxy_v_, i);
375    if (!d)
376      throw utility::GSL_error("vector::operator()",GSL_EINVAL);
377    return *d;
378  }
379
380
381  double& vector::operator[](size_t i)
382  {
383    return this->operator()(i);
384  }
385
386
387  const double& vector::operator[](size_t i) const
388  {
389    return this->operator()(i);
390  }
391
392
393  bool vector::operator==(const vector& other) const
394  {
395    return equal(other);
396  }
397
398
399  bool vector::operator!=(const vector& other) const
400  {
401    return !equal(other);
402  }
403
404
405  double vector::operator*( const vector &other ) const
406  {
407    double res = 0.0;;
408    for ( size_t i = 0; i < size(); ++i ) 
409      res += other(i) * (*this)(i);
410    return res;
411  }
412
413
414  const vector& vector::operator=( const vector& other )
415  {
416    assert(v_);
417    if (gsl_vector_memcpy(v_,other.gsl_vector_p()))
418      throw utility::GSL_error("vector::set dimension mis-match");
419    return *this;
420  } 
421
422
423  const vector& vector::operator+=(const vector& other)
424  {
425    assert(v_);
426    int status=gsl_vector_add(v_, other.gsl_vector_p());
427    if (status)
428      throw utility::GSL_error(std::string("vector::add", status));
429    return *this;
430  }
431
432
433  const vector& vector::operator+=(double d)
434  {
435    assert(v_);
436    gsl_vector_add_constant(v_, d);
437    return *this;
438  }
439
440
441  const vector& vector::operator-=(const vector& other)
442  {
443    assert(v_);
444    int status=gsl_vector_sub(v_, other.gsl_vector_p());
445    if (status)
446      throw utility::GSL_error(std::string("vector::sub", status));
447    return *this;
448  }
449
450
451  const vector& vector::operator-=(const double d)
452  {
453    assert(v_);
454    gsl_vector_add_constant(v_, -d);
455    return *this;
456  }
457
458
459  const vector& vector::operator*=(const double d)
460  {
461    assert(v_);
462    gsl_vector_scale(v_, d);
463    return *this;
464  }
465
466
467  bool isnull(const vector& v)
468  {
469    return gsl_vector_isnull(v.gsl_vector_p());
470  }
471
472
473  double max(const vector& v)
474  {
475    return gsl_vector_max(v.gsl_vector_p());
476  }
477
478
479  size_t max_index(const vector& v)
480  {
481    return gsl_vector_max_index(v.gsl_vector_p());
482  }
483
484
485  double min(const vector& v)
486  {
487    return gsl_vector_min(v.gsl_vector_p());
488  }
489
490
491  size_t min_index(const vector& v)
492  {
493    return gsl_vector_min_index(v.gsl_vector_p());
494  }
495
496
497  bool nan(const vector& templat, vector& flag)
498  {
499    size_t vsize=templat.size();
500    if (vsize!=flag.size())
501      flag.clone(vector(vsize,1.0));
502    else
503      flag.all(1.0);
504    bool nan=false;
505    for (size_t i=0; i<vsize; i++)
506      if (std::isnan(templat(i))) {
507        flag(i)=0;
508        nan=true;
509      }
510    return nan;
511  }
512
513
514  void set_basis(vector& v, size_t i)
515  {
516    assert(v.gsl_vector_p());
517    gsl_vector_set_basis(v.gsl_vector_p(),i);
518  }
519
520
521  void shuffle(vector& invec)
522  {
523    vector vec(invec);
524    random::DiscreteUniform rnd;
525    for (size_t i=0; i<vec.size(); i++){
526      size_t index = rnd(vec.size()-i);
527      invec[i]=vec(index);
528      vec(index)=vec(vec.size()-i-1);
529    }
530  }
531
532
533  void sort(vector& v)
534  {
535    assert(v.gsl_vector_p());
536    std::sort(v.begin(), v.end());
537  }
538
539
540  void sort_index(std::vector<size_t>& sort_index, const vector& invec)
541  {
542    assert(invec.gsl_vector_p());
543    gsl_permutation* p = gsl_permutation_alloc(invec.size());
544    int status=gsl_sort_vector_index(p,invec.gsl_vector_p());
545    if (status) {
546      gsl_permutation_free(p);
547      throw utility::GSL_error(std::string("sort_index(vector&,const vector&)",status));     
548    }
549    sort_index=std::vector<size_t>(p->data,p->data+p->size);
550    gsl_permutation_free(p);
551  }
552
553
554  void sort_smallest_index(std::vector<size_t>& sort_index, size_t k, 
555                            const vector& invec)
556  {
557    assert(invec.gsl_vector_p());
558    assert(k<=invec.size());
559    sort_index.resize(k);
560    gsl_sort_vector_smallest_index(&sort_index[0],k,invec.gsl_vector_p());
561  }
562 
563  void sort_largest_index(std::vector<size_t>& sort_index, size_t k, 
564                            const vector& invec)
565  {
566    assert(invec.gsl_vector_p());
567    assert(k<=invec.size());
568    sort_index.resize(k);
569    gsl_sort_vector_largest_index(&sort_index[0],k,invec.gsl_vector_p());
570  }
571
572
573  double sum(const vector& v)
574  {
575    double sum = 0;
576    size_t vsize=v.size();
577    for (size_t i=0; i<vsize; ++i)
578      sum += v[i];
579    return sum;
580  }
581
582
583  void swap(vector& v, vector& w)
584  {
585    assert(v.gsl_vector_p()); assert(w.gsl_vector_p());
586    int status=gsl_vector_swap(v.gsl_vector_p(),w.gsl_vector_p());
587    if (status)
588      throw utility::GSL_error(std::string("swap(vector&,vector&)",status));
589  }
590
591
592  std::ostream& operator<<(std::ostream& s, const vector& a)
593  {
594    s.setf(std::ios::dec);
595    s.precision(12);
596    for (size_t j = 0; j < a.size(); ++j) {
597      s << a[j];
598      if ( (j+1)<a.size() )
599        s << s.fill();
600    }
601    return s;
602  }
603
604}}} // of namespace utility, yat, and thep
Note: See TracBrowser for help on using the repository browser.