codebase / stimlib

  #ifndef RTS_VECTOR_H

  #define RTS_VECTOR_H

  

  #include <iostream>

  #include <cmath>

  #include <sstream>

  #include <vector>

  #include "../cuda/callable.h"

  namespace stim

  {

  

  

  

  template <class T>

  struct vec : public std::vector<T>

  {

  	using std::vector<T>::size;

  	using std::vector<T>::at;

  	using std::vector<T>::resize;

  	using std::vector<T>::push_back;

  	vec(){

  	}

  

  //	//efficiency constructors, makes construction easier for 1D-4D vectors

  	vec(T rhs)

  	{

  		resize(1, rhs);

  //		cerr << " Created a vector " << endl;

  //		//for(int i=0; i<N; i++)

  //		//	v[i] = rhs;

  	}

  

  //	vec(int s)

  //	{

  //		resize(s,0);

  //	}	

  

  	vec(T x, T y)

  	{

  		push_back(x);

  		push_back(y);

  	}

  	vec(T x, T y, T z)

  	{

  		push_back(x);

  		push_back(y);

  		push_back(z);

  	}

  	vec(T x, T y, T z, T w)

  	{

  		push_back(x);

  		push_back(y);

  		push_back(z);

  		push_back(w);

  	}

  

  	//copy constructor

  	vec( const vec<T>& other){

  		unsigned int N = other.size();

  		for(int i=0; i<N; i++)

  			push_back(other[i]);

  	}

  

  	vec<T> push(T x)

  	{

  		push_back(x);

  		return *this;

  	}

  	vec<T> push(T x, T y)

  	{

  		push_back(x);

  		push_back(y);

  		return *this;

  	}

  	vec<T> push(T x, T y, T z)

  	{

  		push_back(x);

  		push_back(y);

  		push_back(z);

  		return *this;

  	}

  	vec<T> push(T x, T y, T z, T w)

  	{

  		push_back(x);

  		push_back(y);

  		push_back(z);

  		push_back(w);

  		return *this;

  	}

  	template< typename U >

  	operator vec<U>(){

  		unsigned int N = size();

  		vec<U> result;

  		for(int i=0; i<N; i++)

  			result.push_back(at(i));

  

  		return result;

  	}

  

  	//template<class U>

  	//friend vec<U, N>::operator vec<T, N>();

  

  	T len() const

  	{

  		unsigned int N = size();

  

          //compute and return the vector length

          T sum_sq = (T)0;

          for(int i=0; i<N; i++)

          {

              sum_sq += pow( at(i), 2 );

          }

          return sqrt(sum_sq);

  

  	}

  

  	vec<T> cart2sph() const

  	{

  		//convert the vector from cartesian to spherical coordinates

  		//x, y, z -> r, theta, phi (where theta = 0 to 2*pi)

  

  		vec<T> sph;

  		sph.push_back(std::sqrt(at(0)*at(0) + at(1)*at(1) + at(2)*at(2)));

  		sph.push_back(std::atan2(at(1), at(0)));

  

  		if(sph[0] == 0)

  			sph.push_back(0);

  		else

  			sph.push_back(std::acos(at(2) / sph[0]));

  

  		return sph;

  	}

  

  	vec<T> sph2cart() const

  	{

  		//convert the vector from cartesian to spherical coordinates

  		//r, theta, phi -> x, y, z (where theta = 0 to 2*pi)

  

  		vec<T> cart;

  		cart.push_back(at(0) * std::cos(at(1)) * std::sin(at(2)));

  		cart.push_back(at(0) * std::sin(at(1)) * std::sin(at(2)));

  		cart.push_back(at(0) * std::cos(at(2)));

  

  		return cart;

  	}

  

  	vec<T> norm() const

  	{

  		unsigned int N = size();

  

          //compute and return the unit vector

          vec<T> result;

  

          //compute the vector length

          T l = len();

  

          //normalize

          for(int i=0; i<N; i++)

          {

              result.push_back(at(i) / l);

          }

  

          return result;

  	}

  

  	vec<T> cross(const vec<T> rhs) const

  	{

  		vec<T> result;

  

  		//compute the cross product (only valid for 3D vectors)

  		result.push_back(at(1) * rhs[2] - at(2) * rhs[1]);

  		result.push_back(at(2) * rhs[0] - at(0) * rhs[2]);

  		result.push_back(at(0) * rhs[1] - at(1) * rhs[0]);

  

  		return result;

  	}

  

      T dot(vec<T> rhs) const

      {

          T result = (T)0;

  	unsigned int N = size();

          for(int i=0; i<N; i++)

              result += at(i) * rhs[i];

  

          return result;

  

      }

  

  	//arithmetic

  	vec<T> operator+(vec<T> rhs) const

  	{

  		vec<T> result;

  		unsigned int N = size();

  		for(int i=0; i<N; i++)

  		    result.push_back(at(i) + rhs[i]);

  		return result;

  	}

  	vec<T> operator+(T rhs) const

  	{

  		vec<T> result;

  		unsigned int N = size();

  		for(int i=0; i<N; i++)

  		    result.push_back(at(i) + rhs);

  

  		return result;

  	}

  	vec<T> operator-(vec<T> rhs) const

  	{

          vec<T> result;

  		unsigned int N = size();

          for(int i=0; i<N; i++)

              result.push_back(at(i) - rhs[i]);

  

          return result;

  	}

  	vec<T> operator*(T rhs) const

  	{

          vec<T> result;

  		unsigned int N = size();

          for(int i=0; i<N; i++)

              result.push_back(at(i) * rhs);

  

          return result;

  	}

  	vec<T> operator/(T rhs) const

  	{

          vec<T> result;

  		unsigned int N = size();

          for(int i=0; i<N; i++)

              result.push_back(at(i) / rhs);

  

          return result;

  	}

  	vec<T> operator*=(T rhs){

  

  		unsigned int N = size();

  		for(int i=0; i<N; i++)

  			at(i) = at(i) * rhs;

  		return *this;

  	}

  	vec<T> operator+=(vec<T> rhs){

  		unsigned int N = size();

  		for(int i=0; i<N; i++)

  			at(i) += rhs[i];

  		return *this;

  	}

  	vec<T> & operator=(T rhs){

  

  		unsigned int N = size();

  		for(int i=0; i<N; i++)

  			at(i) = rhs;

  		return *this;

  	}

  

  	template<typename Y>

  	vec<T> & operator=(vec<Y> rhs){

  		unsigned int N = rhs.size();

  		resize(N);

  

  		for(int i=0; i<N; i++)

  			at(i) = rhs[i];

  		return *this;

  	}

  	//unary minus

  	vec<T> operator-() const{

  		vec<T> r;

  		//negate the vector

  		unsigned int N = size();

  		for(int i=0; i<N; i++)

  		    r.push_back(-at(i));

  		return r;

  	}

  	std::string str() const

  	{

  		std::stringstream ss;

  

  		unsigned int N = size();

  

  		ss<<"[";

  		for(int i=0; i<N; i++)

  		{

  			ss<<at(i);

  			if(i != N-1)

  				ss<<", ";

  		}

  		ss<<"]";

  

  		return ss.str();

  	}

  

  };

  

  

  }	//end namespace rts

  

  template <typename T>

  std::ostream& operator<<(std::ostream& os, stim::vec<T> v)

  {

      os<<v.str();

      return os;

  }

  

  template <typename T>

  stim::vec<T> operator*(T lhs, stim::vec<T> rhs)

  {

  	stim::vec<T> r;

  

      return rhs * lhs;

  }

  

  //#if __GNUC__ > 3 && __GNUC_MINOR__ > 7

  //template<class T, int N> using rtsVector = rts::vector<T, N>;

  //#endif

  

  #endif