plugins/base1/se.lu.thep.wenni/trunk/lib/c++_tools/gslapi/vector.cc

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

Generated on Tue Nov 21 03:42:21 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)