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

// $Id: Vector.cc 1275 2008-04-11 06:10:12Z jari $

/*
  Copyright (C) 2003 Daniel Dalevi, Peter Johansson
  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
  Copyright (C) 2005, 2006, 2007 Jari Häkkinen, Peter Johansson, Markus Ringnér
  Copyright (C) 2008 Jari Häkkinen, Peter Johansson

  This file is part of the yat library, http://trac.thep.lu.se/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 2 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  02111-1307, USA.
*/

#include "Vector.h"
#include "utility.h"

#include <algorithm>
#include <cassert>
#include <cmath>
#include <iostream>
#include <sstream>
#include <utility>
#include <vector>

namespace theplu {
namespace yat {
namespace utility {


  Vector::Vector(void)
    : VectorMutable()
  {
  }


  Vector::Vector(size_t n, double init_value)
    : VectorMutable(gsl_vector_alloc(n))
  {
    if (!vec_)
      throw utility::GSL_error("Vector::Vector failed to allocate memory");
    assert(const_vec_);
    all(init_value);
  }


  Vector::Vector(const Vector& other)
    : VectorMutable(create_gsl_vector_copy(other))
  {
  }


  Vector::Vector(const VectorBase& other)
    : VectorMutable(create_gsl_vector_copy(other))
  {
  }


  Vector::Vector(std::istream& is, char sep)
    throw (utility::IO_error, std::exception)
    : VectorMutable()
  {
    if (!is.good())
      throw utility::IO_error("Vector: istream is not good");

    // read the data file and store in stl vectors (dynamically
    // expandable)
    std::vector<std::vector<double> > data_matrix;
    unsigned int nof_columns=0;
    unsigned int nof_rows=0;
    std::string line;
    while(getline(is, line, '\n')) {
      // Empty lines
      if (!line.size())
          continue;
      nof_rows++;

      std::vector<double> v;
      std::string element;
      std::stringstream ss(line);
      bool ok=true;
      while(ok) {
        if(sep=='\0')
          ok=(ss>>element);
        else 
          ok=getline(ss, element, sep);
        if(!ok)
          break;        

        if(utility::is_double(element)) {
          v.push_back(atof(element.c_str()));
        }
        else if (!element.size() || utility::is_nan(element)) {
          v.push_back(std::numeric_limits<double>::quiet_NaN());
        }
        else {
          std::stringstream ss("Warning: '");
          ss << element << "' is not a double.";
          throw IO_error(ss.str());
        }
      } 
      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
          v.push_back(std::numeric_limits<double>::quiet_NaN());
      if (!nof_columns)
        nof_columns=v.size();
      else if (nof_rows && (nof_columns>1)) {
        std::ostringstream s;
        s << "Vector::Vector(std::istream&) data file error:\n"
          << "    File has inconsistent number of rows (" << nof_rows
          << ") and columns (" << nof_columns
          << ").\n    Expected a row or a column Vector.";
        throw utility::IO_error(s.str());
      }
      else if (v.size()!=nof_columns) {
        std::ostringstream s;
        s << "Vector::Vector(std::istream&) data file error:\n"
          << "    Line " << nof_rows << " has " << v.size()
          << " columns; expected " << nof_columns << " column.";
        throw utility::IO_error(s.str());
      }
      data_matrix.push_back(v);
    }

    // manipulate the state of the stream to be good
    is.clear(std::ios::goodbit);

    if (!nof_columns || !nof_rows)
      return;

    // convert the data to a gsl vector
    vec_ = gsl_vector_alloc(nof_rows*nof_columns);
    if (!vec_)
      throw utility::GSL_error("Vector::Vector failed to allocate memory");
    size_t n=0;
    // if gsl error handler disabled, out of bounds index will not
    // abort the program.
    for (size_t i=0; i<nof_rows; i++)
      for (size_t j=0; j<nof_columns; j++) 
        gsl_vector_set( vec_, n++, data_matrix[i][j] );
    const_vec_ = vec_;
  }


  Vector::~Vector(void)
  {
    delete_allocated_memory();
  }


  const Vector& Vector::assign(const VectorBase& other)
  {
    // avoid self assignment
    if (this == &other)
      return *this;
    if (!other.size()){
      delete_allocated_memory();
      return *this;
    }
    // indirect self assignment if begin <= other.begin() < end 
    if (begin() <= other.begin() && other.begin()<end()) {
      if (size() != other.size()) {
        gsl_vector* tmp = create_gsl_vector_copy(other);
        delete_allocated_memory();
        const_vec_ = vec_ = tmp;
      }
      return *this;
    } 
    if (size()==other.size()){
      if (gsl_vector_memcpy(vec_, other.gsl_vector_p()))
        throw utility::GSL_error("Vector::assign memcpy failed");
    }
    else {
      delete_allocated_memory();
      vec_ = create_gsl_vector_copy(other);
    }
    const_vec_ = vec_;
    return *this;
  } 


  gsl_vector* Vector::create_gsl_vector_copy(const VectorBase& other) const
  {
    if (!other.size())
      return NULL;
    gsl_vector* vec = gsl_vector_alloc(other.size());
    if (!vec)
      throw utility::GSL_error("Vector::create_gsl_vector_copy failed to allocate memory");
    if (gsl_vector_memcpy(vec, other.gsl_vector_p()))
      throw utility::GSL_error("Vector::create_gsl_vector_copy memcpy failed");
    return vec;
  }


  void Vector::delete_allocated_memory(void) 
  {
    if (vec_)
      gsl_vector_free(vec_);
    const_vec_ = vec_ = NULL;
  }


  bool Vector::isview(void) const
  {
    return false;
  }


  void Vector::resize(size_t n, double init_value)
  {
    delete_allocated_memory();
    if (!n)
      return;
    const_vec_ = vec_ = gsl_vector_alloc(n);
    if (!vec_)
      throw utility::GSL_error("Vector::resize failed to allocate memory");
    all(init_value);
  }


  void swap(Vector& v, Vector& w)
  {
    assert(v.gsl_vector_p()); assert(w.gsl_vector_p());
    int status=gsl_vector_swap(v.gsl_vector_p(),w.gsl_vector_p());
    if (status)
      throw utility::GSL_error(std::string("swap(Vector&,Vector&)",status));
  }


  const Vector& Vector::operator=( const VectorBase& other )
  {
    return assign(other);
  } 


  const Vector& Vector::operator=( const Vector& other )
  {
    return assign(other);
  }


}}} // of namespace utility, yat, and thep