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

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

fixes #434

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