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

Last change on this file since 1680 was 1680, checked in by Peter, 12 years ago

correcting includes

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