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

Last change on this file since 1172 was 1172, checked in by Peter, 15 years ago

allowing resize to zero size. In Matrix if resize is called with only one of columns or rows set to zero a zero by zero matrix is created. I think we can live with that. Though an assertion is called because it is strange.

  • 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 1172 2008-02-27 14:32:15Z 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/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 "utility.h"
29#include "yat/random/random.h"
30
31#include <algorithm>
32#include <cassert>
33#include <cmath>
34#include <iostream>
35#include <sstream>
36#include <utility>
37#include <vector>
38
39namespace theplu {
40namespace yat {
41namespace utility {
42
43
44  Vector::Vector(void)
45    : VectorMutable()
46  {
47  }
48
49
50  Vector::Vector(size_t n, double init_value)
51    : VectorMutable(gsl_vector_alloc(n))
52  {
53    if (!vec_)
54      throw utility::GSL_error("Vector::Vector failed to allocate memory");
55    assert(const_vec_);
56    all(init_value);
57  }
58
59
60  Vector::Vector(const Vector& other)
61    : VectorMutable(create_gsl_vector_copy(other))
62  {
63  }
64
65
66  Vector::Vector(const VectorBase& other)
67    : VectorMutable(create_gsl_vector_copy(other))
68  {
69  }
70
71
72  Vector::Vector(std::istream& is, char sep)
73    throw (utility::IO_error, std::exception)
74    : VectorMutable()
75  {
76    if (!is.good())
77      throw utility::IO_error("Vector: istream is not good");
78
79    // read the data file and store in stl vectors (dynamically
80    // expandable)
81    std::vector<std::vector<double> > data_matrix;
82    u_int nof_columns=0;
83    u_int nof_rows=0;
84    std::string line;
85    while(getline(is, line, '\n')) {
86      // Empty lines
87      if (!line.size())
88          continue;
89      nof_rows++;
90
91      std::vector<double> v;
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(utility::is_double(element)) {
104          v.push_back(atof(element.c_str()));
105        }
106        else if (!element.size() || utility::is_nan(element)) {
107          v.push_back(std::numeric_limits<double>::quiet_NaN());
108        }
109        else {
110          std::stringstream ss("Warning: '");
111          ss << element << "' is not a double.";
112          throw IO_error(ss.str());
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    gsl_vector* vec = gsl_vector_alloc(other.size());
197    if (!vec)
198      throw utility::GSL_error("Vector::create_gsl_vector_copy failed to allocate memory");
199    if (gsl_vector_memcpy(vec, other.gsl_vector_p()))
200      throw utility::GSL_error("Vector::create_gsl_vector_copy memcpy failed");
201    return vec;
202  }
203
204
205  void Vector::delete_allocated_memory(void) 
206  {
207    if (vec_)
208      gsl_vector_free(vec_);
209    const_vec_ = vec_ = NULL;
210  }
211
212
213  bool Vector::isview(void) const
214  {
215    return false;
216  }
217
218
219  void Vector::resize(size_t n, double init_value)
220  {
221    delete_allocated_memory();
222    if (!n)
223      return;
224    const_vec_ = vec_ = gsl_vector_alloc(n);
225    if (!vec_)
226      throw utility::GSL_error("Vector::resize failed to allocate memory");
227    all(init_value);
228  }
229
230
231  void set_basis(Vector& v, size_t i)
232  {
233    assert(v.gsl_vector_p());
234    gsl_vector_set_basis(v.gsl_vector_p(),i);
235  }
236
237
238  void sort(Vector& v)
239  {
240    assert(v.gsl_vector_p());
241    std::sort(v.begin(), v.end());
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.