rtcVec.h

Go to the documentation of this file.

/*
 * Copyright (C) 2009
 * Robert Bosch LLC
 * Research and Technology Center North America
 * Palo Alto, California
 *
 * All rights reserved.
 *
 *------------------------------------------------------------------------------
 * project ....: Autonomous Technologies
 * file .......: rtcVec.h
 * authors ....: Benjamin Pitzer
 * organization: Robert Bosch LLC
 * creation ...: 08/16/2006
 * modified ...: $Date: 2009-03-13 16:48:41 -0700 (Fri, 13 Mar 2009) $
 * changed by .: $Author: wg75pal $
 * revision ...: $Revision: 99 $
 */
#ifndef RTC_VEC_H
#define RTC_VEC_H

//== INCLUDES ==================================================================
#include <vector>
#include "rtc/rtcMath.h"
#include "rtc/rtcMat.h"

#ifdef WIN32
#undef min
#undef max
#endif
//== NAMESPACES ================================================================
namespace rtc {

// Forward declarations
template <class T, int M> class Vec;        // M-d vector
template <class T, int M, int N> class Mat; // MxN Matrix
template <class T, int M> class SMat;       // MxM Square Matrix

template <class T, int M>
class Vec {
public:
  // Constructors
  Vec();
  Vec(const T* d);
  Vec(const T a);
  Vec(const Vec<T,M>& v);

  // Cast Operation
  template <class U> Vec(const Vec<U,M>& v);

  // Mutators
  void set(const T* d);
  void set(const T a);
  void set(const Vec<T,M>& v);
  Vec<T,M>& operator = (const T a);
  Vec<T,M>& operator = (const T* d);
  Vec<T,M>& operator = (const Vec<T,M>& v);
  template <int N> void setSubVec(const int i, const Vec<T,N>& v);

  // Casting Mutators
  template <class U> void set(const Vec<U,M>& v);
  template <class U> Vec<T,M>& operator = (const Vec<U,M>& v);

  // Accessors
  T operator () (const int i) const;
  T& operator () (const int i);
  T operator [] (const int i) const;
  T& operator [] (const int i);
  T& at(const int i);
  const T& at(const int i) const;

  template <int N> Vec<T,N> getSubVec(const int i) const;

  // Norms and Normalize
  T norm() const;
  T normSqr() const;
  T pnorm(float p) const;
  T normalize();
  Vec<T,M> normalized() const;

  // Reductions: Max/Min, Sum/Product
  T max() const;
  T min() const;
  T sum() const;
  T prod() const;
  Vec<T,M> cumsum() const;

  // General elementwise operations
  void perform(T (*mathFun)(T));
  void perform(T (*mathFun)(T,T), const T arg2);
  void perform(T (*mathFun)(T,T), const Vec<T,M>& v);
  Vec<T,M> performed(T (*mathFun)(T));
  Vec<T,M> performed(T (*mathFun)(T,T), const T arg2);
  Vec<T,M> performed(T (*mathFun)(T,T), const Vec<T,M>& v);
  Vec<T,M> minimize(const Vec<T,M>& other);
  Vec<T,M> maximize(const Vec<T,M>& other);

  // Dot and Outer Products
  T dot(const Vec<T,M>& v) const;
  SMat<T,M> outer(const Vec<T,M>& v) const;
  SMat<T,M> outer() const;

  // Random vectors and sorting
  static Vec<T,M> uniformRand(const T a = T(0), const T b = T(1));
  static Vec<T,M> normalRand(const T mean = T(0), const T stdev = T(1));
  static Vec<T,M> multivariateGauss(const Vec<T,M>& mean, const SMat<T,M>& cov);

  void sort(bool ascending = true);
  Vec<T,M> sorted(bool ascending = true);

  // Addition and subtraction
  Vec<T,M>& add(const Vec<T,M>& v);
  Vec<T,M>& subtract(const Vec<T,M>& v);

  // Addition and subtraction operator
  Vec<T,M> operator + (const Vec<T,M>& v) const;
  void operator += (const Vec<T,M>& v);
  Vec<T,M> operator - (const Vec<T,M>& v) const;
  void operator -= (const Vec<T,M>& v);
  Vec<T,M> operator - () const;
  Vec<T,M> add(const T a);
  Vec<T,M> operator + (const T a) const;
  void operator += (const T a);
  Vec<T,M> subtract(const T a);
  Vec<T,M> operator - (const T a) const;
  void operator -= (const T a);

  // Multiplication and division operator
  Vec<T,M> operator * (const T a) const;
  Vec<T,M> operator / (const T a) const;
  void operator *= (const T a);
  void operator /= (const T a);

  // Equality and inequality tests
  int compareTo(const Vec<T,M>& v) const;
  bool equalTo(const Vec<T,M>& v, const T tol = T(0)) const;

  // Equality and inequality tests operator
  bool operator == (const Vec<T,M>& v) const;
  bool operator != (const Vec<T,M>& v) const;
  Vec<bool,M> operator >= (const Vec<T,M>& v) const;
  Vec<bool,M> operator <= (const Vec<T,M>& v) const;
  Vec<bool,M> operator > (const Vec<T,M>& v) const;
  Vec<bool,M> operator < (const Vec<T,M>& v) const;

  // Element-wise operations
  Vec<T,M> operator * (const Vec<T,M>& v) const;
  Vec<T,M> operator / (const Vec<T,M>& v) const;
  void operator *= (const Vec<T,M>& v);
  void operator /= (const Vec<T,M>& v);

  // Serialization
  bool write(OutputHandler& oh) const;
  bool read(InputHandler& ih);

  //data
  T x[M]; 
};

// Global operators for cases where Vec<T,M>
// is the second argument in a binary operator
template <class T, int M> Vec<T,M> operator + (const T a, const Vec<T,M>& v);
template <class T, int M> Vec<T,M> operator - (const T a, const Vec<T,M>& v);
template <class T, int M> Vec<T,M> operator * (const T a, const Vec<T,M>& v);

// ASCII stream IO
template <class T, int M> std::ostream& operator<<(std::ostream& os, const Vec<T,M>& vec);
template <class T, int M> std::istream& operator>>(std::istream& is, Vec<T,M>& vec);

//==============================================================================
// Vec<T,M>
//==============================================================================

// Constructors

template <class T, int M>
inline Vec<T,M>::Vec() {}

template <class T, int M>
inline Vec<T,M>::Vec(const T* d) {
  set(d);
}

template <class T, int M>
inline Vec<T,M>::Vec(const T a) {
  set(a);
}

template <class T, int M>
inline Vec<T,M>::Vec(const Vec<T,M>& v) {
  set(v);
}

// Casting Operation

template <class T, int M> template <class U>
inline Vec<T,M>::Vec(const Vec<U,M>& v) {
  set<U>(v);
}

// Mutators

template <class T, int M>
inline void Vec<T,M>::set(const T a) {
  if (a == T(0))
    memset(x,0,M*sizeof(T));
  else
    for (int i=0;i<M;++i)
      x[i] = a;
}

template <class T, int M>
inline void Vec<T,M>::set(const T* d) {
  for (int i=0;i<M;i++) x[i] = d[i];
}

template <class T, int M>
inline void Vec<T,M>::set(const Vec<T,M>& v) {
  memcpy((void*)x,(void*)v.x,M*sizeof(T));
}

template <class T, int M>
inline Vec<T,M>& Vec<T,M>::operator = (const T a) {
  set(a);
  return *this;
}

template <class T, int M>
inline Vec<T,M>& Vec<T,M>::operator = (const T* d) {
  set(d);
  return *this;
}

template <class T, int M>
inline Vec<T,M>& Vec<T,M>::operator = (const Vec<T,M>& v) {
  set(v);
  return *this;
}

template <class T, int M> template <int N>
inline void Vec<T,M>::setSubVec(const int i, const Vec<T,N>& v) {
  for (int j=i,k=0;j<M && k<N;j++,k++) x[j] = v.x[k];
}

// Casting Mutators

template <class T, int M> template <class U>
inline void Vec<T,M>::set(const Vec<U,M>& v) {
  for (int i=0;i<M;i++) x[i] = T(v.x[i]);
}

template <class T, int M> template <class U>
inline Vec<T,M>& Vec<T,M>::operator = (const Vec<U,M>& v) {
  set<U>(v);
  return *this;
}

// Accessors

template <class T, int M>
inline T Vec<T,M>::operator () (const int i) const {
#if MATMATH_CHECK_BOUNDS
  if (i < 0 || i > M) {
    std::cerr << "Vec<" << M << ">::(" << i << "): index out of range\n";
    std::cerr << std::flush;
    exit(1);
}
#endif
  return x[i];
}

template <class T, int M>
inline T& Vec<T,M>::operator () (const int i) {
#if MATMATH_CHECK_BOUNDS
  if (i < 0 || i > M) {
    std::cerr << "&Vec<" << M << ">::(" << i << "): index out of range\n";
    std::cerr << std::flush;
    exit(1);
  }
#endif
  return x[i];
}

template <class T, int M>
inline T Vec<T,M>::operator [] (const int i) const {
#if MATMATH_CHECK_BOUNDS
  if (i < 0 || i > M) {
    std::cerr << "Vec<" << M << ">::(" << i << "): index out of range\n";
    std::cerr << std::flush;
    exit(1);
}
#endif
  return x[i];
}

template <class T, int M>
inline T& Vec<T,M>::operator [] (const int i) {
#if MATMATH_CHECK_BOUNDS
  if (i < 0 || i > M) {
    std::cerr << "&Vec<" << M << ">::(" << i << "): index out of range\n";
    std::cerr << std::flush;
    exit(1);
  }
#endif
  return x[i];
}

template <class T, int M>
inline T& Vec<T,M>::at(const int i) {
#if MATMATH_CHECK_BOUNDS
  if (i < 0 || i > M) {
    std::cerr << "&Vec<" << M << ">::(" << i << "): index out of range\n";
    std::cerr << std::flush;
    exit(1);
  }
#endif
  return x[i];
}

template <class T, int M>
inline const T& Vec<T,M>::at(const int i) const {
#if MATMATH_CHECK_BOUNDS
  if (i < 0 || i > M) {
    std::cerr << "&Vec<" << M << ">::(" << i << "): index out of range\n";
    std::cerr << std::flush;
    exit(1);
  }
#endif
  return x[i];
}

template <class T, int M> template <int N>
inline Vec<T,N> Vec<T,M>::getSubVec(const int i) const {
  Vec<T,N> v(T(0));
  for (int j=i,k=0;j<M && k<N;j++,k++) v.x[k] = x[j];
  return v;
}

// Norms and Normalize

template <class T, int M>
inline T Vec<T,M>::norm() const {
  return T(sqrt(double(normSqr())));
}

template <class T, int M>
inline T Vec<T,M>::normSqr() const {
  T sum = T(0);
  for (int i=0;i<M;i++) sum += x[i]*x[i];
  return sum;
}

template <class T, int M>
inline T Vec<T,M>::pnorm(float p) const {
  T sum = T(0);
  for (int i=0;i<M;i++) sum += pow(fabs((double)x[i]), (double)p);
  return pow((double)sum, (double)1.0/p);
}

template <class T, int M>
inline T Vec<T,M>::normalize() {
  T l = norm();
  if (l > T(0)) for (int i=0;i<M;i++) x[i] /= l;
  return l;
}

template <class T, int M>
inline Vec<T,M> Vec<T,M>::normalized() const {
  Vec<T,M> v(x);
  v.normalize();
  return v;
}

// Maximum/Minimum value and Sum

template <class T, int M>
inline T Vec<T,M>::max() const {
  T m = x[0];
  for (int i=1;i<M;i++) if (x[i]>m) m = x[i];
  return m;
}

template <class T, int M>
inline T Vec<T,M>::min() const {
  T m = x[0];
  for (int i=1;i<M;i++) if (x[i]<m) m = x[i];
  return m;
}

template <class T, int M>
inline T Vec<T,M>::sum() const {
  T s = T(0);
  for (int i=0;i<M;i++) s += x[i];
  return s;
}

template <class T, int M>
inline Vec<T,M> Vec<T,M>::cumsum() const {
  Vec<T,M> v;
  T s = T(0);
  for (int i=0;i<M;i++) {
    s += x[i];
    v.x[i]=s;
  }
  return v;
}

template <class T, int M>
inline T Vec<T,M>::prod() const {
  T p = T(1);
  for (int i=0;i<M;i++) p *= x[i];
  return p;
}

// General elementwise operations

template <class T, int M>
inline void Vec<T,M>::perform(T (*mathFun)(T)) {
  for (int i=0;i<M;i++) {
    x[i] = (*mathFun)(x[i]);
  }
}

template <class T, int M>
inline void Vec<T,M>::perform(T (*mathFun)(T,T), const T arg2) {
  for (int i=0;i<M;i++) {
    x[i] = (*mathFun)(x[i], arg2);
  }
}

template <class T, int M>
inline void Vec<T,M>::perform(T (*mathFun)(T,T), const Vec<T,M>& v) {
  for (int i=0;i<M;i++) {
    x[i] = (*mathFun)(x[i], v.x[i]);
  }
}

template <class T, int M>
inline Vec<T,M> Vec<T,M>::performed(T (*mathFun)(T)) {
  Vec<T,M> v;
  for (int i=0;i<M;i++) {
    v.x[i] = (*mathFun)(x[i]);
  }
  return v;
}

template <class T, int M>
inline Vec<T,M> Vec<T,M>::performed(T (*mathFun)(T,T), const T arg2) {
  Vec<T,M> v;
  for (int i=0;i<M;i++) {
    v.x[i] = (*mathFun)(x[i], arg2);
  }
  return v;
}

template <class T, int M>
inline Vec<T,M> Vec<T,M>::performed(T (*mathFun)(T,T), const Vec<T,M>& vp) {
  Vec<T,M> v;
  for (int i=0;i<M;i++) {
    v.x[i] = (*mathFun)(x[i], vp.x[i]);
  }
  return v;
}

template <class T, int M>
inline Vec<T,M> Vec<T,M>::minimize(const Vec<T,M>& other) {
  for (int i=0;i<M;i++) {
    if (other.x[i] < x[i]) x[i] = other.x[i];
  }
  return *this;
}

template <class T, int M>
inline Vec<T,M> Vec<T,M>::maximize(const Vec<T,M>& other) {
  for (int i=0;i<M;i++) {
    if (other.x[i] > x[i]) x[i] = other.x[i];
  }
  return *this;
}

// Dot, Cross, and Outer Products

template <class T, int M>
inline T Vec<T,M>::dot(const Vec<T,M>& v) const {
  T sum = T(0);
  for (int i=0;i<M;i++) sum += x[i]*v.x[i];
  return sum;
}

template <class T, int M>
inline SMat<T,M> Vec<T,M>::outer(const Vec<T,M>& v) const {
  SMat<T,M> m;
  for (int i=0,k=0;i<M;i++) for (int j=0;j<M;j++,k++) m.x[k] = x[i]*v.x[j];
  return m;
}

template <class T, int M>
inline SMat<T,M> Vec<T,M>::outer() const {
  SMat<T,M> m;
  for (int i=0,k=0;i<M;i++) for (int j=0;j<M;j++,k++) m.x[k] = x[i]*x[j];
  return m;
}

// Random vectors and sorting

template <class T, int M>
inline Vec<T,M> Vec<T,M>::uniformRand(const T a, const T b) {
  Vec<T,M> v;
  for (int i=0;i<M;i++) v.x[i] = rtc_uniform_rand<T>(a,b);
  return v;
}

template <class T, int M>
inline Vec<T,M> Vec<T,M>::normalRand(const T mean, const T stdev) {
  Vec<T,M> v;
  for (int i=0;i<M;i++) v.x[i] = rtc_normal_rand<T>(mean, stdev);
  return v;
}

template <class T, int M>
inline Vec<T,M> Vec<T,M>::multivariateGauss(const Vec<T,M>& mean, const SMat<T,M>& cov) {
  Vec<T,M> v;
  SMat<T,M> S(cov);
  int n=S.choleskyDecomp();
  S.transpose();
  Vec<T,M> X = normalRand();
  v = mean + S*X;
  return v;
}


template <class T, int M>
inline void Vec<T,M>::sort(bool ascending) {
  int count;
  do {
    count = 0;
    for (int i=0;i<(M-1);i++) {
if (ascending) {
  if (x[i] > x[i+1]) {
    rtc_swap(x[i],x[i+1]); count++;
  }
} else {
  if (x[i] < x[i+1]) {
    rtc_swap(x[i+1],x[i]); count++;
  }
}
    }
  } while (count > 0);
}

template <class T, int M>
inline Vec<T,M> Vec<T,M>::sorted(bool ascending) {
  Vec<T,M> v(x);
  v.sort(ascending);
  return v;
}

// Addition and subtraction

template <class T, int M>
inline Vec<T,M>& Vec<T,M>::add(const Vec<T,M>& v) {
  for (int i=0;i<M;i++) x[i] = x[i] + v.x[i];
  return *this;
}

template <class T, int M>
inline Vec<T,M>& Vec<T,M>::subtract(const Vec<T,M>& v) {
  for (int i=0;i<M;i++) x[i] = x[i] - v.x[i];
  return *this;
}

// Addition and subtraction operator

template <class T, int M>
inline Vec<T,M> Vec<T,M>::operator + (const Vec<T,M>& v) const {
  Vec<T,M> vp;
  for (int i=0;i<M;i++) vp.x[i] = x[i] + v.x[i];
  return vp;
}

template <class T, int M>
inline void Vec<T,M>::operator += (const Vec<T,M>& v) {
  for (int i=0;i<M;i++) x[i] += v.x[i];
}

template <class T, int M>
inline Vec<T,M> Vec<T,M>::operator - (const Vec<T,M>& v) const {
  Vec<T,M> vp;
  for (int i=0;i<M;i++) vp.x[i] = x[i] - v.x[i];
  return vp;
}

template <class T, int M>
inline void Vec<T,M>::operator -= (const Vec<T,M>& v) {
  for (int i=0;i<M;i++) x[i] -= v.x[i];
}

template <class T, int M>
inline Vec<T,M> Vec<T,M>::operator - () const {
  Vec<T,M> vp;
  for (int i=0;i<M;i++) vp.x[i] = -x[i];
  return vp;
}

template <class T, int M>
inline Vec<T,M> Vec<T,M>::add(const T a) {
  for (int i=0;i<M;i++) x[i] = x[i] - a;
  return *this;
}

template <class T, int M>
inline Vec<T,M> Vec<T,M>::operator + (const T a) const {
  Vec<T,M> vp;
  for (int i=0;i<M;i++) vp.x[i] = x[i] + a;
  return vp;
}

template <class T, int M>
inline void Vec<T,M>::operator += (const T a) {
  for (int i=0;i<M;i++) x[i] += a;
}

template <class T, int M>
inline Vec<T,M> Vec<T,M>::subtract(const T a) {
  for (int i=0;i<M;i++) x[i] = x[i] - a;
  return *this;
}

template <class T, int M>
inline Vec<T,M> Vec<T,M>::operator - (const T a) const {
  Vec<T,M> vp;
  for (int i=0;i<M;i++) vp.x[i] = x[i] - a;
  return vp;
}

template <class T, int M>
inline void Vec<T,M>::operator -= (const T a) {
  for (int i=0;i<M;i++) x[i] -= a;
}

template <class T, int M>
inline Vec<T,M> operator + (const T a, const Vec<T,M>& v) {
  Vec<T,M> vp;
  for (int i=0;i<M;i++) vp.x[i] = a + v.x[i];
  return vp;
}

template <class T, int M>
inline Vec<T,M> operator - (const T a, const Vec<T,M>& v) {
  Vec<T,M> vp;
  for (int i=0;i<M;i++) vp.x[i] = a - v.x[i];
  return vp;
}

// Multiplication and division

template <class T, int M>
inline Vec<T,M> Vec<T,M>::operator * (const T a) const {
  Vec<T,M> vp;
  for (int i=0;i<M;i++) vp.x[i] = x[i]*a;
  return vp;
}

template <class T, int M>
inline Vec<T,M> Vec<T,M>::operator / (const T a) const {
  Vec<T,M> vp;
  for (int i=0;i<M;i++) vp.x[i] = x[i]/a;
  return vp;
}

template <class T, int M>
inline void Vec<T,M>::operator *= (const T a) {
  for (int i=0;i<M;i++) x[i] *= a;
}

template <class T, int M>
inline void Vec<T,M>::operator /= (const T a) {
  for (int i=0;i<M;i++) x[i] /= a;
}

template <class T, int M>
inline Vec<T,M> operator * (const T a, const Vec<T,M>& v) {
  Vec<T,M> vp;
  for (int i=0;i<M;i++) vp.x[i] = a*v.x[i];
  return vp;
}

// Equality and inequality tests.

template <class T, int M>
inline bool Vec<T,M>::operator == (const Vec<T,M>& v) const {
  for (int i=0;i<M;i++)
    if (x[i] != v.x[i])
      return(false);
  return true;
}

template <class T, int M>
inline bool Vec<T,M>::operator != (const Vec<T,M>& v) const {
  for (int i=0;i<M;i++)
    if (x[i] != v.x[i])
      return(true);
  return false;
}

//  /** Element-wise test for equality.
//   * @return vector of boolean results
//   */
//  template <class T, int M>
//  inline Vec<bool,M> Vec<T,M>::operator == (const Vec<T,M>& v) const {
//    Vec<bool,M> b(false);
//    for (int i=0;i<M;i++) if (x[i] == v.x[i]) b.x[i] = true;
//    return b;
//  }

//  /** Element-wise test for inequality.
//   * @return vector of boolean results
//   */
//  template <class T, int M>
//  inline Vec<bool,M> Vec<T,M>::operator != (const Vec<T,M>& v) const {
//    Vec<bool,M> b(false);
//    for (int i=0;i<M;i++) if (x[i] != v.x[i]) b.x[i] = true;
//    return b;
//  }

template <class T, int M>
inline Vec<bool,M> Vec<T,M>::operator >= (const Vec<T,M>& v) const {
  Vec<bool,M> b(false);
  for (int i=0;i<M;i++) if (x[i] >= v.x[i]) b.x[i] = true;
  return b;
}

template <class T, int M>
inline Vec<bool,M> Vec<T,M>::operator <= (const Vec<T,M>& v) const {
  Vec<bool,M> b(false);
  for (int i=0;i<M;i++) if (x[i] <= v.x[i]) b.x[i] = true;
  return b;
}

template <class T, int M>
inline Vec<bool,M> Vec<T,M>::operator > (const Vec<T,M>& v) const {
  Vec<bool,M> b(false);
  for (int i=0;i<M;i++) if (x[i] > v.x[i]) b.x[i] = true;
  return b;
}

template <class T, int M>
inline Vec<bool,M> Vec<T,M>::operator < (const Vec<T,M>& v) const {
  Vec<bool,M> b(false);
  for (int i=0;i<M;i++) if (x[i] < v.x[i]) b.x[i] = true;
  return b;
}

template <class T, int M>
inline int Vec<T,M>::compareTo(const Vec<T,M>& v) const {
  int g=0, l=0;
  for (int i=0;i<M;i++) {
    if (x[i] < v.x[i]) l++;
    if (x[i] > v.x[i]) g++;
  }
  if (l==M) return -1;
  else if (g==M) return 1;
  else return 0;
}

template <class T, int M>
inline bool Vec<T,M>::equalTo(const Vec<T,M>& v, const T tol) const {
  bool t = true;
  for (int i=0;i<M;i++) if (rtc_abs(x[i] - v.x[i]) > tol) t = false;
  return t;
}

// Element-wise multiplication and division

template <class T, int M>
inline Vec<T,M> Vec<T,M>::operator * (const Vec<T,M>& v) const {
  Vec<T,M> vp;
  for (int i=0;i<M;i++) vp.x[i] = x[i]*v.x[i];
  return vp;
}

template <class T, int M>
inline Vec<T,M> Vec<T,M>::operator / (const Vec<T,M>& v) const {
  Vec<T,M> vp;
  for (int i=0;i<M;i++) vp.x[i] = x[i]/v.x[i];
  return vp;
}

template <class T, int M>
inline void Vec<T,M>::operator *= (const Vec<T,M>& v) {
  for (int i=0;i<M;i++) x[i] *= v.x[i];
}

template <class T, int M>
inline void Vec<T,M>::operator /= (const Vec<T,M>& v) {
  for (int i=0;i<M;i++) x[i] /= v.x[i];
}

// Serialization routines

template <class T, int M>
inline bool Vec<T,M>::write(OutputHandler& oh) const {
  return oh.write((char *)(x),M*sizeof(T));
}

template <class T, int M>
inline bool Vec<T,M>::read(InputHandler& ih) {
  return ih.read((char *)(x),M*sizeof(T));
}

template <class T, int M>
std::ostream& operator<<(std::ostream& os, const Vec<T,M>& vec) {
int minFieldWidth = os.precision()+2;

os << "[";
for (int i=0; i<M; ++i)
  os << std::setw(minFieldWidth) << vec.x[i] << " ";
os << "]";

return os;
}

template <class T, int M>
std::istream& operator>>(std::istream& is, Vec<T,M>& vec) {
  using namespace std;
  vector<T> data;
  string vecString;
  stringstream vecStringStream;

  getline(is, vecString, ']');
  int sPos = (int)vecString.find('[');
  if (sPos == (int)string::npos)
    throw Exception("format error: expecting formated vector to start with '['");

  //erase the starting '['
  //note the ending ']' was removed by the getline function as the delim
  vecString.erase(0,sPos+1);
  trim(vecString);
  vecStringStream.str(vecString);

  //determine num of rows and cols
  int colCount = 0;
  T tmpVal;
  while(vecStringStream.good()){
    vecStringStream >> tmpVal;
    data.push_back(tmpVal);
    ++colCount;
  }

  //check that dimensions agree
  if (colCount != M){
    std::stringstream ss;
    ss << "format error: formated text has " << colCount << " columns";
    ss << " while destination vector has " << M  << " columns" << endl;
    throw Exception(ss.str());
  }

  //copy extracted data
  for (int i=0;i<M;i++){
    vec.x[i] = data[i];
  }

  return is;
}

template <class T, int M>
bool rtc_write(OutputHandler& oh, const Vec<T,M>& data)
{
  return data.write(oh);
};

template <class T, int M>
bool rtc_read(InputHandler& ih, Vec<T,M>& data)
{
  return data.read(ih);
};

//==============================================================================
} // NAMESPACE rtc
//==============================================================================
#endif // RTC_VEC_H defined
//==============================================================================