// -*- c++ -*-
///////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2001 Oh-Wook Kwon, all rights reserved. ohwook@yahoo.com
//
//                          Easy Matrix Template Library
// 
// This Easy Matrix Template Library is provided "as is" without any express 
// or implied warranty of any kind with respect to this software. 
// In particular the authors shall not be liable for any direct, 
// indirect, special, incidental or consequential damages arising 
// in any way from use of the software.
// 
// Everyone is granted permission to copy, modify and redistribute this
// Easy Matrix Template Library, provided:
//  1.  All copies contain this copyright notice.
//  2.  All modified copies shall carry a notice stating who
//      made the last modification and the date of such modification.
//  3.  No charge is made for this software or works derived from it.  
//      This clause shall not be construed as constraining other software
//      distributed on the same medium as this software, nor is a
//      distribution fee considered a charge.
//
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
// Filename: Vector.h// Implementation:
//   Template class for vector algebra
//   Vector index starts from 1.
//   Access by A(i) rather than A[i].
// Note: 
//   Access by A(i) loses nothing but may gain something at least in some cases.
//   All return values are of the double type to prevent overflow or underflow.
// If you find any bugs, please send me an e-mail.
///////////////////////////////////////////////////////////////////////////////
#ifndef _VECTOR_H_#define _VECTOR_H_#include "./MatrixConfig.h"
#include "./MathUtil.h"
#include "./Random.h"

typedef complex<float> FComplex;
typedef complex<double> DComplex;
template <typename T> class Vector;template <typename T> class Matrix;enum VectorType { COL_VECTOR, ROW_VECTOR };


// Template functions with default arguments
#ifdef __GNUC__
template <typename T> T Max(const Vector<T>& A,int* idx=0);
template <typename T> T Min(const Vector<T>& A,int* idx=0);
template <typename T> Matrix<T> ToMatrix(const Vector<T>& A, int m=1);
template <typename T> double Norm(const Vector<T>& A, int p=2);
#elsetemplate <typename T> T Max(const Vector<T>& A,int* idx);
template <typename T> T Min(const Vector<T>& A,int* idx);
template <typename T> Matrix<T> ToMatrix(const Vector<T>& A, int m);
template <typename T> double Norm(const Vector<T>& A, int p);
#endif
// friend template functions
template <typename T> Vector<T> mtl_plus(const Vector<T>& A, const Vector<T>& B);
template <typename T> Vector<T> mtl_minus(const Vector<T>& A, const Vector<T>& B);
template <typename T> Vector<T> mtl_plus(const Vector<T>& A, const double b);
template <typename T> Vector<T> mtl_minus(const Vector<T>& A, const double b);
template <typename T> Vector<T> mtl_times(const Vector<T>& A, const double b);
template <typename T> Vector<T> mtl_divide(const Vector<T>& A, const double b);
template <typename T> Vector<T> mtl_plus(const double b, const Vector<T>& A);
template <typename T> Vector<T> mtl_minus(const double b, const Vector<T>& A);
template <typename T> Vector<T> mtl_times(const double b, const Vector<T>& A);
template <typename T> Vector<T> mtl_divide(const double b, const Vector<T>& A);
template <typename T> Vector<T> mtl_times(const Vector<T>& A, const Vector<T>& B);
template <typename T> Vector<T> mtl_divide(const Vector<T>& A, const Vector<T>& B);
template <typename T> ostream& mtl_ostream(ostream& s, const Vector<T>& A);template <typename T> istream& mtl_istream(istream& s, Vector<T>& A);

#ifndef DISABLE_COMPLEX
template <typename T> Vector<T> mtl_plus(const Vector<T>& A, const complex<T> b);
template <typename T> Vector<T> mtl_minus(const Vector<T>& A, const complex<T> b);
template <typename T> Vector<T> mtl_times(const Vector<T>& A, const complex<T> b);
template <typename T> Vector<T> mtl_divide(const Vector<T>& A, const complex<T> b);
template <typename T> Vector<T> mtl_plus(const complex<T> b, const Vector<T>& A);
template <typename T> Vector<T> mtl_minus(const complex<T> b, const Vector<T>& A);
template <typename T> Vector<T> mtl_times(const complex<T> b, const Vector<T>& A);
template <typename T> Vector<T> mtl_divide(const complex<T> b, const Vector<T>& A);
#endif
template <typename T> int Length(const Vector<T>& A);template <typename T> T Sum(const Vector<T>& A);
template <typename T> Vector<T> Transpose(const Vector<T>& A);
template <typename T> Vector<T> ArrayTranspose(const Vector<T>& A);template <typename T> Matrix<T> OuterProd(const Vector<T>& A, const Vector<T>& B);template <typename T> T InnerProd(const Vector<T>& A, const Vector<T>& B);
template <typename T> Vector<T> Cross(const Vector<T>& A, const Vector<T>& B);template <typename T> T Dot(const Vector<T>& A, const Vector<T>& B);
template <typename T> Vector<T> ArrayMultiply(const Vector<T>& A, const Vector<T>& B);template <typename T> Vector<T> ArrayDivide(const Vector<T>& A, const Vector<T>& B);template <typename T> Vector<int> Int(const Vector<T>& A);
template <typename T> Vector<float> Float(const Vector<T>& A);
template <typename T> Vector<double> Double(const Vector<T>& A);
template <typename T> Vector<T> Abs(const Vector<T>& A);
template <typename T> Vector<T> Sign(const Vector<T>& A);
template <typename T> Vector<T> Pow(const Vector<T>& A, double b);
template <typename T> Vector<T> Pow(const Vector<T>& A, int b);
template <typename T> Vector<T> Exp(const Vector<T>& A);template <typename T> Vector<T> Log(const Vector<T>& A);template <typename T> Vector<T> Digamma(const Vector<T>& A);template <typename T> Vector<T> Gammaln(const Vector<T>& A);template <typename T> Matrix<T> Diag(const Vector<T> A);template <typename T> T Mean(const Vector<T>& A);
template <typename T> Vector<T> Cumsum(const Vector<T>& A);
template <typename T> Vector<T> Cumprod(const Vector<T>& A);template <typename T> string Num2str(const Vector<T>& A);template <typename T> int Sort(const Vector<T>& A, Vector<T>& v, Vector<int>& orderA, int order);
template <typename T> Vector<T> FFT(const Vector<T>& A,int N);
template <typename T> Vector<T> IFFT(const Vector<T>& A,int N);
template <typename T> Vector<T> Merge(const Vector<T>& A,const Vector<T>& B);


// ------------------------------------------------------------------------// Vector declaration// ------------------------------------------------------------------------template <typename T>class Vector {private:	// I chose "vector" in STL as a storage because the "vector" preserves
	// data after resizing. I waste one element but it made implementation easy.	int n;
	VectorType type;
	vector<T> vec;
public:
	friend class Matrix<T>;	Vector(int n_=0, VectorType type_=(VectorType)COL_VECTOR, const T& v_=T());
	Vector(int n_, VectorType type_, const char* v_);
	Vector(const Vector<T>& v_, int n_, VectorType type_);	//Vector(const Vector<T>& A);	virtual ~Vector();

	/////////////////////////////
	// member operator overloading
	/////////////////////////////	T&			operator()(int i);
	const T&	operator()(int i) const;
	Vector<T>	operator+();	Vector<T>	operator-();	Vector<T>& operator+=(const Vector<T>& B);	Vector<T>& operator-=(const Vector<T>& B);	Vector<T>& operator+=(const double b);	Vector<T>& operator-=(const double b);	Vector<T>& operator*=(const double b);	Vector<T>& operator/=(const double b);
#ifndef DISABLE_COMPLEX	Vector<T>& operator+=(const complex<T> b);
	Vector<T>& operator-=(const complex<T> b);
	Vector<T>& operator*=(const complex<T> b);
	Vector<T>& operator/=(const complex<T> b);
#endif
	//Vector<T>& operator=(const Vector<T>& A);	Vector<T>& operator=(const T& a);
#ifdef ALLOW_ACCESS_BY_BRACKET
	// The [] operator does not check if the arguments are in the valid range.
	// I recommend to use it only in implementation of the Vector class.
	T&			operator[](int i);
	const T&	operator[](int i) const;
#endif

	/////////////////////////////
	// friend operator overloading
	/////////////////////////////	// GCC-2.95-2 did not allow friend operator overloading definition outside class declaration.
	friend Vector<T> operator+(const Vector<T>& A, const Vector<T>& B) {		return mtl_plus(A,B);	}	friend Vector<T> operator-(const Vector<T>& A, const Vector<T>& B) {		return mtl_minus(A,B);	}	friend Vector<T> operator+(const Vector<T>& A, const double b) {		return mtl_plus(A,b);	}	friend Vector<T> operator-(const Vector<T>& A, const double b) {		return mtl_minus(A,b);	}	friend Vector<T> operator*(const Vector<T>& A, const double b) {		return mtl_times(A,b);	}	friend Vector<T> operator/(const Vector<T>& A, const double b) {		return mtl_divide(A,b);	}	friend Vector<T> operator+(const double b, const Vector<T>& A) {		return mtl_plus(b,A);	}	friend Vector<T> operator-(const double b, const Vector<T>& A) {		return mtl_minus(b,A);	}	friend Vector<T> operator*(const double b, const Vector<T>& A) {		return mtl_times(b,A);	}	friend Vector<T> operator/(const double b, const Vector<T>& A) {		return mtl_divide(b,A);	}	friend ostream& operator<<(ostream& s, const Vector<T>& A) {		return mtl_ostream(s,A);	}
	friend ostream& operator<<(ostream& s, Vector<T>& A) {
		return mtl_ostream(s,A);
	}	friend istream& operator>>(istream& s, Vector<T>& A) {		return mtl_istream(s,A);	}

#ifndef DISABLE_COMPLEX
	friend Vector<T> operator+(const Vector<T>& A, const complex<T> b) {
		return mtl_plus(A,b);
	}
	friend Vector<T> operator-(const Vector<T>& A, const complex<T> b) {
		return mtl_minus(A,b);
	}
	friend Vector<T> operator*(const Vector<T>& A, const complex<T> b) {
		return mtl_times(A,b);
	}
	friend Vector<T> operator/(const Vector<T>& A, const complex<T> b) {
		return mtl_divide(A,b);
	}
	friend Vector<T> operator+(const complex<T> b, const Vector<T>& A) {
		return mtl_plus(b,A);
	}
	friend Vector<T> operator-(const complex<T> b, const Vector<T>& A) {
		return mtl_minus(b,A);
	}
	friend Vector<T> operator*(const complex<T> b, const Vector<T>& A) {
		return mtl_times(b,A);
	}
	friend Vector<T> operator/(const complex<T> b, const Vector<T>& A) {
		return mtl_divide(b,A);
	}
#endif
	/////////////////////////////
	// friend functions
	/////////////////////////////
#ifdef __GNUC__
	friend Vector<T> mtl_plus<T>(const Vector<T>& A, const Vector<T>& B);
	friend Vector<T> mtl_minus<T>(const Vector<T>& A, const Vector<T>& B);
	friend Vector<T> mtl_plus<T>(const Vector<T>& A, const double b);
	friend Vector<T> mtl_minus<T>(const Vector<T>& A, const double b);
	friend Vector<T> mtl_times<T>(const Vector<T>& A, const double b);
	friend Vector<T> mtl_divide<T>(const Vector<T>& A, const double b);
	friend Vector<T> mtl_plus<T>(const double b, const Vector<T>& A);
	friend Vector<T> mtl_minus<T>(const double b, const Vector<T>& A);
	friend Vector<T> mtl_times<T>(const double b, const Vector<T>& A);
	friend Vector<T> mtl_divide<T>(const double b, const Vector<T>& A);
	friend Vector<T> mtl_times<T>(const Vector<T>& A, const Vector<T>& B);
	friend Vector<T> mtl_divide<T>(const Vector<T>& A, const Vector<T>& B);
	friend ostream& mtl_ostream<T>(ostream& s, const Vector<T>& A);	friend istream& mtl_istream<T>(istream& s, Vector<T>& A);

#ifndef DISABLE_COMPLEX
	friend Vector<T> mtl_plus<T>(const Vector<T>& A, const complex<T> b);
	friend Vector<T> mtl_minus<T>(const Vector<T>& A, const complex<T> b);
	friend Vector<T> mtl_times<T>(const Vector<T>& A, const complex<T> b);
	friend Vector<T> mtl_divide<T>(const Vector<T>& A, const complex<T> b);
	friend Vector<T> mtl_plus<T>(const complex<T> b, const Vector<T>& A);
	friend Vector<T> mtl_minus<T>(const complex<T> b, const Vector<T>& A);
	friend Vector<T> mtl_times<T>(const complex<T> b, const Vector<T>& A);
	friend Vector<T> mtl_divide<T>(const complex<T> b, const Vector<T>& A);
#endif
	friend int Length<T>(const Vector<T>& A);	friend T Sum<T>(const Vector<T>& A);	friend Vector<T> Transpose<T>(const Vector<T>& A);	friend Vector<T> ArrayTranspose<T>(const Vector<T>& A);
	friend Matrix<T> OuterProd<T>(const Vector<T>& A, const Vector<T>& B);	friend T InnerProd<T>(const Vector<T>& A, const Vector<T>& B);
	friend Vector<T> Cross<T>(const Vector<T>& A, const Vector<T>& B);
	friend T Dot<T>(const Vector<T>& A, const Vector<T>& B);	friend Vector<T> ArrayMultiply<T>(const Vector<T>& A, const Vector<T>& B);	friend Vector<T> ArrayDivide<T>(const Vector<T>& A, const Vector<T>& B);
	friend Vector<int> Int<T>(const Vector<T>& A);	friend Vector<float> Float<T>(const Vector<T>& A);
	friend Vector<double> Double<T>(const Vector<T>& A);
	friend Vector<T> Abs<T>(const Vector<T>& A);
	friend Vector<T> Sign<T>(const Vector<T>& A);
	friend Vector<T> Pow<T>(const Vector<T>& A, double b);
	friend Vector<T> Pow<T>(const Vector<T>& A, int b);
	friend Vector<T> Exp<T>(const Vector<T>& A);	friend Vector<T> Log<T>(const Vector<T>& A);	friend Vector<T> Digamma<T>(const Vector<T>& A);	friend Vector<T> Gammaln<T>(const Vector<T>& A);	friend Matrix<T> Diag<T>(const Vector<T> A);	friend T Mean<T>(const Vector<T>& A);
	friend Vector<T> Cumsum<T>(const Vector<T>& A);
	friend Vector<T> Cumprod<T>(const Vector<T>& A);	friend string Num2str<T>(const Vector<T>& A);	friend int Sort<T>(const Vector<T>& A, Vector<T>& v, Vector<int>& orderA, int order);
	friend Vector<T> FFT<T>(const Vector<T>& A,int N);
	friend Vector<T> IFFT<T>(const Vector<T>& A,int N);
	friend Vector<T> Merge<T>(const Vector<T>& A,const Vector<T>& B);#else
	friend Vector<T> mtl_plus(const Vector<T>& A, const Vector<T>& B);
	friend Vector<T> mtl_minus(const Vector<T>& A, const Vector<T>& B);
	friend Vector<T> mtl_plus(const Vector<T>& A, const double b);
	friend Vector<T> mtl_minus(const Vector<T>& A, const double b);
	friend Vector<T> mtl_times(const Vector<T>& A, const double b);
	friend Vector<T> mtl_divide(const Vector<T>& A, const double b);
	friend Vector<T> mtl_plus(const double b, const Vector<T>& A);
	friend Vector<T> mtl_minus(const double b, const Vector<T>& A);
	friend Vector<T> mtl_times(const double b, const Vector<T>& A);
	friend Vector<T> mtl_divide(const double b, const Vector<T>& A);
	friend Vector<T> mtl_times(const Vector<T>& A, const Vector<T>& B);
	friend Vector<T> mtl_divide(const Vector<T>& A, const Vector<T>& B);
	friend ostream& mtl_ostream(ostream& s, const Vector<T>& A);	friend istream& mtl_istream(istream& s, Vector<T>& A);

#ifndef DISABLE_COMPLEX
	friend Vector<T> mtl_plus(const Vector<T>& A, const complex<T> b);
	friend Vector<T> mtl_minus(const Vector<T>& A, const complex<T> b);
	friend Vector<T> mtl_times(const Vector<T>& A, const complex<T> b);
	friend Vector<T> mtl_divide(const Vector<T>& A, const complex<T> b);
	friend Vector<T> mtl_plus(const complex<T> b, const Vector<T>& A);
	friend Vector<T> mtl_minus(const complex<T> b, const Vector<T>& A);
	friend Vector<T> mtl_times(const complex<T> b, const Vector<T>& A);
	friend Vector<T> mtl_divide(const complex<T> b, const Vector<T>& A);
#endif
	friend int Length(const Vector<T>& A);	friend T Sum(const Vector<T>& A);
	friend Vector<T> Transpose(const Vector<T>& A);	friend Vector<T> ArrayTranspose(const Vector<T>& A);
	friend Matrix<T> OuterProd(const Vector<T>& A, const Vector<T>& B);	friend T InnerProd(const Vector<T>& A, const Vector<T>& B);
	friend Vector<T> Cross(const Vector<T>& A, const Vector<T>& B);
	friend T Dot(const Vector<T>& A, const Vector<T>& B);	friend Vector<T> ArrayMultiply(const Vector<T>& A, const Vector<T>& B);	friend Vector<T> ArrayDivide(const Vector<T>& A, const Vector<T>& B);	friend Vector<int> Int(const Vector<T>& A);
	friend Vector<float> Float(const Vector<T>& A);
	friend Vector<double> Double(const Vector<T>& A);
	friend Vector<T> Abs(const Vector<T>& A);
	friend Vector<T> Sign(const Vector<T>& A);	friend Vector<T> Pow(const Vector<T>& A, double b);
	friend Vector<T> Pow(const Vector<T>& A, int b);
	friend Vector<T> Exp(const Vector<T>& A);	friend Vector<T> Log(const Vector<T>& A);	friend Vector<T> Digamma(const Vector<T>& A);	friend Vector<T> Gammaln(const Vector<T>& A);	friend Matrix<T> Diag(const Vector<T> A);	friend T Mean(const Vector<T>& A);
	friend Vector<T> Cumsum(const Vector<T>& A);
	friend Vector<T> Cumprod(const Vector<T>& A);	friend string Num2str(const Vector<T>& A);	friend int Sort(const Vector<T>& A, Vector<T>& v, Vector<int>& orderA, int order);
	friend Vector<T> FFT(const Vector<T>& A,int N);
	friend Vector<T> IFFT(const Vector<T>& A,int N);
	friend Vector<T> Merge(const Vector<T>& A,const Vector<T>& B);

#endif	
	/////////////////////////////
	// member functions
	/////////////////////////////	int Size() const;
	VectorType Type() const;
	int SetType(VectorType type_);
	Vector<T> t() const;	Vector<T> ct() const;
	int Resize(int new_n);	int Add(const T& a);	int Erase(int c);	void Clear();
	bool Empty() const;
	T Sum_() const;	T Max_(int* idx=0) const;	T Min_(int* idx=0) const;
	int Print(ostream& s, const char *name=0);	int Print(ostream& s, string& name);
	int Read(istream& s, string& name);
	int Read(istream& s, const char *name=0);	Vector<T> SubVector(int n1, int n2);	int Save(const string& fname);
	int SaveMatlab(const string& fname);
	void SetPermIdentity();

	int IsOne(double tolerance=1e-6) const;
	int IsZero(double tolerance=1e-6) const;
	int IsUnit(double tolerance=1e-6) const;

	void FromInt(const Vector<int>& A);
	void FromFloat(const Vector<float>& A);
	void FromDouble(const Vector<double>& A);
	Vector<int> ToInt();
	Vector<float> ToFloat();
	Vector<double> ToDouble();

public:
	/////////////////////////////
	// exported static functions
	/////////////////////////////
	static Vector<T> Rand(int m);
	static Vector<T> Randn(int m);
	static Vector<T> CRandn(int m);

	/////////////////////////////
	// local static functions
	/////////////////////////////
	static double normp(const Vector<T>& A,int p);
	static int FFTBase(const Vector<T>& A,int N,int isign,Vector<T>& C);
	static int FFTCore(Vector<T>& data,int nn,int isign);
};typedef Vector<float> FVector;typedef Vector<double> DVector;typedef Vector<FComplex> CFVector;
typedef Vector<DComplex> CDVector;
// ------------------------------------------------------------------------// Vector implementation// ------------------------------------------------------------------------#ifndef local_max
#define local_max(a, b)  (((a) > (b)) ? (a) : (b))
#endif
#ifndef local_min
#define local_min(a, b)  (((a) < (b)) ? (a) : (b)) 
#endif


/////////////////////////////
// constructor and destructor
/////////////////////////////template <typename T> Vector<T>::Vector(int n_,VectorType type_,const T& v_) : n(n_),type(type_) {
	vec.resize(n_+1);	for(int i=1;i<=n_;i++) vec[i]=v_;}


template <typename T> Vector<T>::Vector(const Vector<T>& v_,int n_,VectorType type_) : n(n_),type(type_) {
	vec.resize(n_+1);
	int i;
	for(i=1;i<=n_;i++) {
		if(i<=v_.Size()) vec[i]=v_[i];
		else vec[i]=T();
	}
}


template <typename T> Vector<T>::Vector(int n_,VectorType type_,const char* v_) : n(n_),type(type_) {
	vec.resize(n_+1);
	Vector<string> num=mtl_split(v_);
	for(int i=1;i<=n_;i++) {
		if(i<=num.Size()) vec[i]=T(atof(num[i].c_str()));