// FILE ex_vector2.cc

#include <iostream>
#include <fstream>
#include <cmath>
#include <cassert>
#include <valarray>

using namespace std;



class Vector
{
private:
  valarray<double> vec_;
  
public:
  Vector (int size)
    : vec_(size)
  { }
  
  // assignment
  Vector& operator= (const Vector& v)
  {
    if (vec_.size () != v.vec_.size ())
      {
	vec_.resize (v.vec_.size ());
      }
    vec_ = v.vec_;
  }
  
  // set all elements to value
  Vector& operator= (double value)
  {
    vec_ = value;
  }
  
  // get size
  int size () const
  {
    return (int)vec_.size ();
  }
  
  // resize and initialize with 0
  void resize (int newsize)
  {
    vec_.resize ((size_t)newsize);
  }
  
  // element access
  double operator[] (int i) const { assert (0<=i && i<size()); return vec_[i]; }
  double& operator[] (int i) { assert (0<=i && i<size()); return vec_[i]; }
  
  void print (ostream& os) const
  {
    os << '[';
    if (vec_.size () > 0) os << vec_[0];
    for (size_t i=1; i<vec_.size(); ++i) os << ", " << vec_[i];
    os << ']';
  }
  
  Vector& operator+= (const Vector& v)
  {
    assert (v.size() == size());
    vec_ += v.vec_;
    return *this;
  }
  
  Vector& operator-= (const Vector& v)
  {
    assert (v.size() == size());
    vec_ -= v.vec_;
    return *this;
  }
  
  Vector& operator*= (double factor)
  {
    vec_ *= factor;
    return *this;
  }
  
  double l2norm () const
  {
    double sum = 0;
    for (int i=0; i<size(); ++i)
      {
	sum += vec_[i]*vec_[i];
      }
    return sqrt(sum);
  }
  
}; // class Vector


inline
ostream& operator<< (ostream& os, const Vector& v)
{
  v.print (os);
  return os;
}


double inner_product (const Vector& v, const Vector& w)
{
  assert (v.size () > 0 && v.size() == w.size());
  double sum = 0;
  for (int i=0; i<v.size(); ++i)
    {
      sum += v[i] * w[i];
    }
  return sum;
}



int main ()
{
  Vector x(5);
  Vector y(5);
  
  x = 0; // initialize all (5) elements with 0
  x[0] = 1;
  x[2] = 4;
  cout << "x = " << x << endl;
  
  y = 1; // initialize all (5) elements with 1
  y[1] = -3;
  y[2] = 3;
  cout << "y = " << y << endl;
  
  x += y;
  cout << "x += y -> x = " << x << endl;
  
  cout << "||x|| = " << x.l2norm () << endl;
  
  cout << "(x, y) = " << inner_product (x, y) << endl;

  for (int i=0; i<5; ++i)
    {
      x[i] = i;
    }
  cout << "x = " << x << endl;
  
  cout << "writing x and y to file 'data.gpd'..." << endl;
  // file output

  ofstream ofs;
  ofs.open("data.gpd");
  ofs << "#x   y" << endl;
  for (int i=0; i<5; ++i)
    {
      ofs << x[i] << " " << y[i] << endl;
    }
  ofs.close ();
  cout << "done" << endl;
  cout << "in Linux, these data can be visualized using gnuplot" << endl;
  cout << "plot 'data.gpd' using 1:2 w lp" << endl;
  
}