stim/math/vector.h

#ifndef STIM_VECTOR_H
#define STIM_VECTOR_H

#include <iostream>
#include <cmath>
#include <sstream>
#include <vector>

#include <stim/cuda/cudatools/callable.h>
#include <stim/math/vec3.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(){

	}

	/// Create a vector with a set dimension d
	vec(size_t d)
	{
		resize(d,0);
	}


//	//efficiency constructors, makes construction easier for 1D-4D vectors
	vec(T x, T y)
	{
		resize(2, 0);
		at(0) = x;
		at(1) = y;
	}
	vec(T x, T y, T z)
	{
		resize(3, 0);
		at(0) = x;
		at(1) = y;
		at(2) = z;
	}
	vec(T x, T y, T z, T w)
	{
		resize(4, 0);
		at(0) = x;
		at(1) = y;
		at(2) = z;
		at(3) = w;
	}

	vec(std::string str){
		std::stringstream ss(str);

		T c;
		while(ss >> c){
			push_back(c);
		}

	}

	
	//copy constructor
	vec( const vec<T>& other){
		size_t N = other.size();
		resize(N);							//resize the current vector to match the copy
		for(size_t i=0; i<N; i++){	//copy each element
			at(i) = other[i];
		}
	}

	//I'm not sure what these were doing here.
	//Keep them now, we'll worry about it later.
	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;
	}

	/// Casting operator. Creates a new vector with a new type U.
	template< typename U >
	operator vec<U>(){
		size_t N = size();
		vec<U> result;
		for(int i=0; i<N; i++)
			result.push_back(at(i));

		return result;
	}


	/// computes the Euclidean length of the vector
	T len() const
	{
		size_t N = size();

        //compute and return the vector length
        T sum_sq = (T)0;
        for(size_t i=0; i<N; i++)
        {
            sum_sq += pow( at(i), 2 );
        }
        return sqrt(sum_sq);

	}

	
	vec<T> cyl2cart() const
	{
		vec<T> cyl;
		cyl.push_back(at(0)*std::sin(at(1)));
		cyl.push_back(at(0)*std::cos(at(1)));
		cyl.push_back(at(2));
		return(cyl);
		
	}
	/// Convert the vector from cartesian to spherical coordinates (x, y, z -> r, theta, phi where theta = [0, 2*pi])
	vec<T> cart2sph() const
	{


		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;
	}

	/// Convert the vector from cartesian to spherical coordinates (r, theta, phi -> x, y, z where theta = [0, 2*pi])
	vec<T> sph2cart() const
	{
		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;
	}

	/// Computes the normalized vector (where each coordinate is divided by the L2 norm)
	vec<T> norm() const
	{
		size_t N = size();

        //compute and return the unit vector
        vec<T> result;

        //compute the vector length
        T l = len();

        //normalize
        for(size_t i=0; i<N; i++)
        {
            result.push_back(at(i) / l);
        }

        return result;
	}

	/// Computes the cross product of a 3-dimensional vector
	vec<T> cross(const vec<T> rhs) const
	{

		vec<T> result(3);

		//compute the cross product (only valid for 3D vectors)
		result[0] = (at(1) * rhs[2] - at(2) * rhs[1]);
		result[1] = (at(2) * rhs[0] - at(0) * rhs[2]);
		result[2] = (at(0) * rhs[1] - at(1) * rhs[0]);

		return result;
	}

	/// Compute the Euclidean inner (dot) product
    T dot(vec<T> rhs) const
    {
        T result = (T)0;
        size_t N = size();
        for(int i=0; i<N; i++)
            result += at(i) * rhs[i];

        return result;

    }

	/// Arithmetic addition operator

    /// @param rhs is the right-hand-side operator for the addition
	vec<T> operator+(vec<T> rhs) const
	{
		size_t N = size();
		vec<T> result(N);

		for(int i=0; i<N; i++)
		    result[i] = at(i) + rhs[i];

		return result;
	}

	/// Arithmetic addition to a scalar

	/// @param rhs is the right-hand-side operator for the addition
	vec<T> operator+(T rhs) const
	{
		size_t N = size();

		vec<T> result(N);
		for(int i=0; i<N; i++)
		    result[i] = at(i) + rhs;

		return result;
	}

	/// Arithmetic subtraction operator

	/// @param rhs is the right-hand-side operator for the subtraction
	vec<T> operator-(vec<T> rhs) const
	{
		size_t N = size();

        vec<T> result(N);

        for(size_t i=0; i<N; i++)
            result[i] = at(i) - rhs[i];

        return result;
	}
	/// Arithmetic subtraction to a scalar

	/// @param rhs is the right-hand-side operator for the addition
	vec<T> operator-(T rhs) const
	{
		size_t N = size();

		vec<T> result(N);
		for(size_t i=0; i<N; i++)
		    result[i] = at(i) - rhs;

		return result;
	}

	/// Arithmetic scalar multiplication operator

	/// @param rhs is the right-hand-side operator for the subtraction
	vec<T> operator*(T rhs) const
	{
		size_t N = size();

        vec<T> result(N);

        for(size_t i=0; i<N; i++)
            result[i] = at(i) * rhs;

        return result;
	}

	/// Arithmetic scalar division operator

	/// @param rhs is the right-hand-side operator for the subtraction
	vec<T> operator/(T rhs) const
	{
		size_t N = size();

        vec<T> result(N);

        for(size_t i=0; i<N; i++)
            result[i] = at(i) / rhs;

        return result;
	}

	/// Multiplication by a scalar, followed by assignment
	vec<T> operator*=(T rhs){

		size_t N = size();
		for(size_t i=0; i<N; i++)
			at(i) = at(i) * rhs;
		return *this;
	}

	/// Addition and assignment
	vec<T> operator+=(vec<T> rhs){
		size_t N = size();
		for(size_t i=0; i<N; i++)
			at(i) += rhs[i];
		return *this;
	}

	/// Assign a scalar to all values
	vec<T> & operator=(T rhs){

		size_t N = size();
		for(size_t i=0; i<N; i++)
			at(i) = rhs;
		return *this;
	}

	/// Cast to a vec3
	operator stim::vec3<T>(){
		stim::vec3<T> r;
		size_t N = std::min<size_t>(size(), 3);
		for(size_t i = 0; i < N; i++)
			r[i] = at(i);
		return r;
	}

	/// Casting and assignment
	template<typename Y>
	vec<T> & operator=(vec<Y> rhs){
		size_t N = rhs.size();
		resize(N);

		for(size_t i=0; i<N; i++)
			at(i) = rhs[i];
		return *this;
	}

	/// Unary minus (returns the negative of the vector)
	vec<T> operator-() const{

		size_t N = size();

		vec<T> r(N);

		//negate the vector
		for(size_t i=0; i<N; i++)
		    r[i] = -at(i);

		return r;
	}


	/// Outputs the vector as a string
	std::string str() const
	{
		std::stringstream ss;

		size_t N = size();

		ss<<"[";
		for(size_t 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;
}


/// Multiply a vector by a constant when the vector is on the right hand side
template <typename T>
stim::vec<T> operator*(T lhs, stim::vec<T> rhs)
{
	stim::vec<T> r;

    return rhs * lhs;
}

#endif