source: trunk/yat/utility/Vector.cc @ 1797

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

updating copyright statements

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.4 KB
Line 
1// $Id: Vector.cc 1797 2009-02-12 18:07:10Z 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 Jari Häkkinen, Peter Johansson
8  Copyright (C) 2009 Peter Johansson
9
10  This file is part of the yat library, http://dev.thep.lu.se/yat
11
12  The yat library is free software; you can redistribute it and/or
13  modify it under the terms of the GNU General Public License as
14  published by the Free Software Foundation; either version 3 of the
15  License, or (at your option) any later version.
16
17  The yat library is distributed in the hope that it will be useful,
18  but WITHOUT ANY WARRANTY; without even the implied warranty of
19  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20  General Public License for more details.
21
22  You should have received a copy of the GNU General Public License
23  along with yat. If not, see <http://www.gnu.org/licenses/>.
24*/
25
26#include "Vector.h"
27#include "utility.h"
28
29#include <algorithm>
30#include <cassert>
31#include <cmath>
32#include <iostream>
33#include <string>
34#include <sstream>
35#include <utility>
36#include <vector>
37
38namespace theplu {
39namespace yat {
40namespace utility {
41
42
43  Vector::Vector(void)
44    : VectorMutable()
45  {
46  }
47
48
49  Vector::Vector(size_t n, double init_value)
50    : VectorMutable(create_gsl_vector(n))
51  {
52    if (n)
53      all(init_value);
54    assert(n==0 || vec_);
55    assert(n==0 || const_vec_);
56  }
57
58
59  Vector::Vector(const Vector& other)
60    : VectorMutable(create_gsl_vector_copy(other))
61  {
62  }
63
64
65  Vector::Vector(const VectorBase& other)
66    : VectorMutable(create_gsl_vector_copy(other))
67  {
68  }
69
70
71  Vector::Vector(std::istream& is, char sep)
72    throw (utility::IO_error, std::exception)
73    : VectorMutable()
74  {
75    if (!is.good())
76      throw utility::IO_error("Vector: istream is not good");
77
78    // read the data file and store in stl vectors (dynamically
79    // expandable)
80    std::vector<std::vector<double> > data_matrix;
81    unsigned int nof_columns=0;
82    unsigned int nof_rows=0;
83    std::string line;
84    while(getline(is, line, '\n')) {
85      // Empty lines
86      if (!line.size())
87          continue;
88      nof_rows++;
89
90      data_matrix.push_back(std::vector<double>());
91      std::vector<double>& v=data_matrix.back();
92      std::string element;
93      std::stringstream ss(line);
94      bool ok=true;
95      while(ok) {
96        if(sep=='\0')
97          ok=(ss>>element);
98        else 
99          ok=getline(ss, element, sep);
100        if(!ok)
101          break;       
102
103        if (!element.size())
104          v.push_back(std::numeric_limits<double>::quiet_NaN());
105        else {
106          try {
107            v.push_back(convert<double>(element));
108          }
109          catch (std::runtime_error& e) {
110            std::stringstream ss(e.what());
111            ss << "\nVector.cc: " << element
112               << " is not accepted as a Vector element\n";
113            throw IO_error(ss.str());
114          }
115        }
116      } 
117      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
118          v.push_back(std::numeric_limits<double>::quiet_NaN());
119      if (!nof_columns)
120        nof_columns=v.size();
121      else if (nof_rows && (nof_columns>1)) {
122        std::ostringstream s;
123        s << "Vector::Vector(std::istream&) data file error:\n"
124          << "    File has inconsistent number of rows (" << nof_rows
125          << ") and columns (" << nof_columns
126          << ").\n    Expected a row or a column Vector.";
127        throw utility::IO_error(s.str());
128      }
129      else if (v.size()!=nof_columns) {
130        std::ostringstream s;
131        s << "Vector::Vector(std::istream&) data file error:\n"
132          << "    Line " << nof_rows << " has " << v.size()
133          << " columns; expected " << nof_columns << " column.";
134        throw utility::IO_error(s.str());
135      }
136    }
137
138    // manipulate the state of the stream to be good
139    is.clear(std::ios::goodbit);
140
141    if (!nof_columns || !nof_rows)
142      return;
143
144    // convert the data to a gsl vector
145    vec_ = gsl_vector_alloc(nof_rows*nof_columns);
146    if (!vec_)
147      throw utility::GSL_error("Vector::Vector failed to allocate memory");
148    size_t n=0;
149    // if gsl error handler disabled, out of bounds index will not
150    // abort the program.
151    for (size_t i=0; i<nof_rows; i++)
152      for (size_t j=0; j<nof_columns; j++) 
153        gsl_vector_set( vec_, n++, data_matrix[i][j] );
154    const_vec_ = vec_;
155  }
156
157
158  Vector::~Vector(void)
159  {
160    delete_allocated_memory();
161  }
162
163
164  const Vector& Vector::assign(const VectorBase& other)
165  {
166    // avoid self assignment
167    if (this == &other)
168      return *this;
169    if (!other.size()){
170      delete_allocated_memory();
171      return *this;
172    }
173    // indirect self assignment if begin <= other.begin() < end
174    if (begin() <= other.begin() && other.begin()<end()) {
175      if (size() != other.size()) {
176        gsl_vector* tmp = create_gsl_vector_copy(other);
177        delete_allocated_memory();
178        const_vec_ = vec_ = tmp;
179      }
180      return *this;
181    } 
182    if (size()==other.size()){
183      if (gsl_vector_memcpy(vec_, other.gsl_vector_p()))
184        throw utility::GSL_error("Vector::assign memcpy failed");
185    }
186    else {
187      delete_allocated_memory();
188      vec_ = create_gsl_vector_copy(other);
189    }
190    const_vec_ = vec_;
191    return *this;
192  } 
193
194
195  gsl_vector* Vector::create_gsl_vector_copy(const VectorBase& other) const
196  {
197    if (!other.size())
198      return NULL;
199    gsl_vector* vec = gsl_vector_alloc(other.size());
200    if (!vec)
201      throw utility::GSL_error("Vector::create_gsl_vector_copy failed to allocate memory");
202    if (gsl_vector_memcpy(vec, other.gsl_vector_p()))
203      throw utility::GSL_error("Vector::create_gsl_vector_copy memcpy failed");
204    return vec;
205  }
206
207
208  gsl_vector* Vector::create_gsl_vector(size_t n) const
209  {
210    if (!n)
211      return NULL;
212    gsl_vector* vec = gsl_vector_alloc(n);
213    if (!vec)
214      throw utility::GSL_error("Vector::create_gsl_vector failed to allocate memory");
215    return vec;
216  }
217
218
219  void Vector::delete_allocated_memory(void) 
220  {
221    if (vec_)
222      gsl_vector_free(vec_);
223    const_vec_ = vec_ = NULL;
224  }
225
226
227  bool Vector::isview(void) const
228  {
229    return false;
230  }
231
232
233  void Vector::resize(size_t n, double init_value)
234  {
235    delete_allocated_memory();
236    if (!n)
237      return;
238    const_vec_ = vec_ = gsl_vector_alloc(n);
239    if (!vec_)
240      throw utility::GSL_error("Vector::resize failed to allocate memory");
241    all(init_value);
242  }
243
244
245  void swap(Vector& v, Vector& w)
246  {
247    assert(v.gsl_vector_p()); assert(w.gsl_vector_p());
248    int status=gsl_vector_swap(v.gsl_vector_p(),w.gsl_vector_p());
249    if (status)
250      throw utility::GSL_error(std::string("swap(Vector&,Vector&)",status));
251  }
252
253
254  const Vector& Vector::operator=( const VectorBase& other )
255  {
256    return assign(other);
257  } 
258
259
260  const Vector& Vector::operator=( const Vector& other )
261  {
262    return assign(other);
263  }
264
265
266}}} // of namespace utility, yat, and thep
Note: See TracBrowser for help on using the repository browser.