source: trunk/c++_tools/gslapi/vector.cc @ 593

Last change on this file since 593 was 593, checked in by Markus Ringnér, 15 years ago

Fixed std includes to compile with g++ 4.1.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.6 KB
Line 
1// $Id: vector.cc 593 2006-08-25 12:59:21Z markus $
2
3/*
4  Copyright (C) 2003 Daniel Dalevi, Peter Johansson
5  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
6  Copyright (C) 2005 Jari Häkkinen, Peter Johansson, Markus Ringnér
7  Copyright (C) 2006 Jari Häkkinen, Peter Johansson
8
9  This file is part of the thep c++ tools library,
10                                http://lev.thep.lu.se/trac/c++_tools
11
12  The c++ tools library is free software; you can redistribute it
13  and/or modify it under the terms of the GNU General Public License
14  as published by the Free Software Foundation; either version 2 of
15  the License, or (at your option) any later version.
16
17  The c++ tools library is distributed in the hope that it will be
18  useful, but WITHOUT ANY WARRANTY; without even the implied warranty
19  of 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 this program; if not, write to the Free Software
24  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
25  02111-1307, USA.
26*/
27
28#include <c++_tools/gslapi/vector.h>
29#include <c++_tools/gslapi/matrix.h>
30#include <c++_tools/utility/stl_utility.h>
31#include <c++_tools/utility/utility.h>
32
33
34#include <iostream>
35#include <sstream>
36#include <vector>
37#include <utility>
38
39
40namespace theplu {
41namespace gslapi {
42
43
44
45  vector::vector(vector& v, size_t offset, size_t n, size_t stride)
46    : const_view_(NULL)
47  {
48    // Jari, exception handling needed here. Failure in setting up a
49    // proper gsl_vector_view is communicated by NULL pointer in the
50    // view structure (cf. GSL manual). How about GSL error state?
51    view_ = new gsl_vector_view(gsl_vector_subvector_with_stride(v.v_,offset,
52                                                                 stride,n));
53    v_ = &(view_->vector);
54  }
55
56
57
58  vector::vector(const vector& v, size_t offset, size_t n, size_t stride)
59    : view_(NULL)
60  {
61    // Jari, exception handling needed here. Failure in setting up a
62    // proper gsl_vector_view is communicated by NULL pointer in the
63    // view structure (cf. GSL manual). How about GSL error state?
64    const_view_ = new gsl_vector_const_view(
65                   gsl_vector_const_subvector_with_stride(v.v_,offset,stride,n));
66    v_ = const_cast<gsl_vector*>(&(const_view_->vector));
67  }
68
69
70
71  vector::vector(matrix& m, size_t i, bool row)
72    : const_view_(NULL)
73  {
74    view_=new gsl_vector_view(row ?
75                              gsl_matrix_row   (m.gsl_matrix_p(),i) :
76                              gsl_matrix_column(m.gsl_matrix_p(),i)  );
77    v_ = &(view_->vector);
78  }
79
80
81
82  vector::vector(const matrix& m, size_t i, bool row)
83    : view_(NULL)
84  {
85    const_view_ = new gsl_vector_const_view(row ?
86                                   gsl_matrix_const_row   (m.gsl_matrix_p(),i) :
87                                   gsl_matrix_const_column(m.gsl_matrix_p(),i) );
88    v_ = const_cast<gsl_vector*>(&(const_view_->vector));
89  }
90
91
92
93  vector::vector(std::istream& is, char sep) 
94    throw (utility::IO_error,std::exception)
95    : view_(NULL), const_view_(NULL)
96  {
97    // read the data file and store in stl vectors (dynamically
98    // expandable)
99    std::vector<std::vector<double> > data_matrix;
100    u_int nof_columns=0;
101    u_int nof_rows=0;
102    std::string line;
103    while(getline(is, line, '\n')) {
104      // Empty lines
105      if (!line.size())
106          continue;
107      nof_rows++;
108     
109      std::vector<double> v;
110      std::string element;
111      std::stringstream ss(line);
112      bool ok=true;
113      while(ok) {
114        if(sep=='\0')
115          ok=(ss>>element);
116        else 
117          ok=getline(ss, element, sep);
118        if(!ok)
119          break;       
120
121        if(utility::is_double(element)) {
122          v.push_back(atof(element.c_str()));
123        }
124        else if (!element.size() || utility::is_nan(element)) {
125          v.push_back(std::numeric_limits<double>::quiet_NaN());
126        }
127        else {
128          // Jari, this should be communicated with as an exception.
129          // std::cerr << "Warning: '" << element
130          //           << "' is not an integer." << std::endl;
131        }
132      } 
133      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
134          v.push_back(std::numeric_limits<double>::quiet_NaN());
135      if (!nof_columns)
136        nof_columns=v.size();
137      else if (nof_rows && (nof_columns>1)) {
138        std::ostringstream s;
139        s << "vector::vector(std::istream&) data file error:\n"
140          << "    File has inconsistent number of rows (" << nof_rows
141          << ") and columns (" << nof_columns
142          << ").\n    Expected a row or a column vector.";
143        throw utility::IO_error(s.str());
144      }
145      else if (v.size()!=nof_columns) {
146        std::ostringstream s;
147        s << "vector::vector(std::istream&) data file error:\n"
148          << "    Line " << nof_rows << " has " << v.size()
149          << " columns; expected " << nof_columns << " column.";
150        throw utility::IO_error(s.str());
151      }
152      data_matrix.push_back(v);
153    }
154
155    // manipulate the state of the stream to be good
156    is.clear(std::ios::goodbit);
157    // convert the data to a gsl vector
158    v_ = gsl_vector_alloc(nof_rows*nof_columns);
159    size_t n=0;
160    for (size_t i=0; i<nof_rows; i++)
161      for (size_t j=0; j<nof_columns; j++) 
162        gsl_vector_set( v_, n++, data_matrix[i][j] );
163  }
164
165
166
167  vector::~vector(void)
168  {
169    if (view_)
170      delete view_;
171    else if (const_view_)
172      delete const_view_;
173    else if (v_)
174      gsl_vector_free(v_);
175    v_=NULL;
176  }
177
178
179
180  gsl_vector* vector::create_gsl_vector_copy(void) const
181  {
182    gsl_vector* vec = gsl_vector_alloc(size());
183    gsl_vector_memcpy(vec,v_);
184    return vec;
185  }
186
187
188
189  std::pair<double,double> vector::minmax(void) const
190  {
191    double min, max;
192    gsl_vector_minmax(v_,&min,&max);
193    return std::pair<double,double>(min,max);
194  }
195
196
197
198  std::pair<size_t,size_t> vector::minmax_index(void) const
199  {
200    size_t min_index, max_index;
201    gsl_vector_minmax_index(v_,&min_index,&max_index);
202    return std::pair<size_t,size_t>(min_index,max_index);
203  }
204
205
206
207  double vector::sum(void) const
208  {
209    double sum = 0;
210    for (size_t i=0; i<size(); i++)
211      sum += (*this)(i);
212    return( sum );
213  } 
214
215  bool vector::operator==(const vector& a) const
216  {
217    if (size()!=a.size())
218      return false;
219    for (size_t i=0; i<size(); ++i)
220      if (gsl_vector_get(v_,i)!=a(i))
221        return false;
222    return true;
223  }
224
225
226
227  double vector::operator*( const vector &other ) const
228  {
229    double res = 0.0;;
230    for ( size_t i = 0; i < size(); ++i ) 
231      res += other(i) * (*this)(i);
232    return res;
233  }
234
235
236
237  const vector& vector::operator=( const vector& other )
238  {
239    if( this != &other ) {
240      if (view_)
241        delete view_;
242      else if (const_view_)
243        delete const_view_;
244      else if ( v_ )
245        gsl_vector_free( v_ );
246      v_ = other.create_gsl_vector_copy();
247    }
248    return *this;
249  } 
250
251
252
253  std::ostream& operator<<(std::ostream& s, const vector& a)
254  {
255    s.setf(std::ios::dec);
256    s.precision(12);
257    for (size_t j = 0; j < a.size(); ++j) {
258      s << a[j];
259      if ( (j+1)<a.size() )
260        s << s.fill();
261    }
262
263    return s;
264  }
265
266
267}} // of namespace gslapi and namespace thep
Note: See TracBrowser for help on using the repository browser.